- Email: [email protected]
Bilinear Continuous-Time Systems Identification via Hartley-Based Modulating Functions
Automatica, Vol. 34, No. 4, pp. 499—503, 1998 ( 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0005-1098/98 $19.00#0.00
PII: S0005–1098(97)00216–1
Technical Communique
Bilinear Continuous-Time Systems Identification via Hartley-Based Modulating Functions* S. DANIEL-BERHE-‡ and H. UNBEHAUENKey Words—Bilinear continuous-time systems; identification; Hartley modulating functions method.
Hsu, 1982), etc., have been used to approximate the I/O-sequences of a bilinear system in order to identify the unknown parameters. The present study is focused towards the bilinear continuous-time systems identification using Hartley modulating functions and is organized as follows. In the next part, the structure of bilinear continuous-time systems is discussed. Then identification of the proposed method is presented followed by a typical illustrative example. Finally, the paper finishes with some discussions and conclusions as well as future research plans.
Abstract—This paper highlights the relevance and merits of the Hartley modulating functions (HMF) method for the identification of bilinear continuous-time (BCT) systems from recorded input and noise-contaminated output data and it provides an insight into parameter estimation of a wider range of nonlinear systems in practice. The methodology replaces the I/O-differential equation representing the dynamic system behavior by the Hartley spectrum equation. As a result it involves the known derivatives of the modulating function instead of the derivatives of the input and noisy output data by applying integral transformation to signals. A frequency weighted least-squares algorithm is also applied in the identification and a normalized root mean square criterion is used to investigate some computational considerations and the bias of the estimates. Results of the simulation studies demonstrate the appropriateness of the approach and its efficiency. ( 1998 Elsevier Science Ltd. All rights reserved.
2. Structure of bilinear continuous-time systems The mathematical description of a continuous-time singleinput—single-output (SISO) bilinear system can be given in the following state-space representation: x5 (t)"Ax(t)#bu(t)#u(t)Nx(t), y(t)"cTx(t),
1. Introduction Clearly, most real systems are nonlinear and continuous in time and a linear representation of their actual behavior may lead to limitations in applications. Correspondingly, nonlinear parametric systems may be modeled using special or general nonlinear models. Bilinear systems are one of the special subclass and well-structured nonlinear systems. In the past few decades attentions have been focused on the class of bilinear systems from the theoretical and the practical point of view (Rink and Mohler, 1968; Mohler, 1973, 1991; Bruni et al., 1974; Haber and Keviczky, 1976, etc.). They arise naturally in a wide variety of applications and are found to be quite suitable for fitting some important actual processes, for instance, nuclear fission, hydraulic drive, turbo-pump dynamics, convection heat transfer, etc., have been modeled approximately by bilinear systems (Mohler, 1973, 1991; Goodhart et al., 1991; Reuter, 1994, etc.). In fact, it appears that many important processes, not only in engineering, but also in biology, socio-economics, ecology, etc., may be modeled by bilinear systems, refer Mohler (1973, 1991), etc. Recently, a bilinear approach for system modeling and control of an industrial plant as well as the characteristics of practical bilinear model structures have been presented by Dunoyer et al. (1996a, b), respectively. Furthermore, the Walsh functions (Rao et al., 1976; Karanam et al., 1978), the block-pulse functions (Cheng and
(1)
where u(t) and y(t) are the measured scalar control input and output, respectively, x(t) is an n-dimensional state vector, A and N are the n]n, matrices of real constants, b and c are the n]1 vectors of real constants. Bilinear systems can be formulated in different forms, however, for the detail mathematical analysis of this presentation it is assumed to be represented in the phase variable canonical form of the following matrices of A, N, b and c: 0 1 . . . 0 . . . . . . . . A" , . . . . 0 . . . 0 1 !a . . . !a !a 1 n~1 n 0 . . N" . 0 n 1
* Received 13 May 1997; received in final form 10 November 1997. This paper was recommended for publication in revised form by Editor P. Dorato. Corresponding author Prof. Dr.-Ing. H. Unbehauen. Tel. #49/234 700 4071; Fax #49/ 234 7094 101; E-mail [email protected]. - Control Engineering Laboratory, Faculty of Electrical Engineering, Ruhr-University Bochum, Building IC 3/131, D-44780 Bochum, Germany. ‡ On German Academic Exchange Service Research Fellowship and on leave from the Department of Electrical Engineering, Addis Ababa University, Addis Ababa, Ethiopia.
. . .
0 . . . . . , . . . . 0 . . . n n
b 1 . . b" . b n~1 b n
(2) 1 . . and C" . . 0 0
Representation of equations (1) and (2) can be rewritten in the following simplified form: xR (t)"x (t)#b u(t), j j`1 j
j"1, 2, 2 , n!1,
n xR (t)" + [!a #n u(t)] x (t)#b u(t), n i i i n i/1 y(t)"x (t). 1 499
(3)
500
Technical Communiques
Now, consider a third-order BCT observability canonical form given by equation (4) and rewrite in measured I/O-variables, i.e. xR (t) 0 1 0 1 xR (t) " 0 0 1 2 xR (t) !a !a !a 3 1 2 3
x (t) 1 x (t) 2 x (t) 3
b 0 0 0 1 # b u(t)#u(t) 0 0 0 2 b n n n 3 1 2 3
G
/ (t), 04t4¹ / (t)" m , m 0 otherwise x (t) 1 x (t) , 2 x (t) 3
/(l) (t) exists for all l"0, 1, 2 , k!1, and m
2
!a [y¨ (t)!b uR (t)!b u(t)] 3 1 2
][yR (t)!b u(t)]#n u(t) 1 3 (5)
Rearranging the terms leads to the preceding equation of the form 3 3 3 j~1 (3) y(t)"! + a #(i~1) y(t) # + b (3~i) u(t) # + + a b (j~1~i) u(t) i i j i i/1 i/1 j/2 i/1 (6) 3 3 j~1 # + n u(t) (i~1) y(t) ! + + n b u(t) (j~1~i) u(t) , i j i i/1 j/2 i/1 (0) where l(t)"(d0/dt0) l(t)"l(t). The above derivation can be generalized for an nth order of the form given in equation (2). The main purpose is to estimate the BCT system parameters a , i b and n , i"1, 2, 3 from the input—output data records. Generi i ally, equation (6) is a standard form of a nonlinear input—output differential operator model (Pearson, 1992) of the form (7)
where # (h) are given functions of the parameter vector h with j # "1, ! (u, y) and ) (u, y) are specified functions of the I/O0 jk k signals, P (p) are fixed polynomials of degree n in the differential jk operator p"d/dt, j"0, 1, 2 , n , and k"0, 1, 2 , n . 1 2 3. Identification In general, the modulating functions method starts with multiplying both sides of equation (6) by / (t) and integrating over m the length of the modulating function which results in
P
P
T
P
3 T 3 T (3) / (t) y(t) dt"!+ a / (t) (i~1) y(t) dt# + b / (t) (3~i) u(t) dt m i m i m 0 0 0 i/1 i/1
P
3 j~1 T #+ + ab / (t) (j~1~i) u(t) dt j i m 0 j/2 i/1
P
where
P
AB
P
T
0
m(t) /(i) (t) dt m
AB AB
in :T m(t) cas((!1)i(k#m!j) u (t) dt 0 0hgggggigggggj 2 i Hm((~1) (k`m~j)u0) k k in " + (!1)j (k#m!j)iui cas@ 0 j 2 j/0
AB
]H ((!1)i (k#m!j)) u , m 0 "HM (i) (mu ) 0 m
(11)
where HM (i) (mu ) is the mth HMF spectral component of the 0 m ith derivative of a continuous signal m(t), and H (mu ) is the m 0 Hartley transform (HT) (Bracewll, 1986) of m(t) defined by
P
H (u)" m
=
m(t) cas ut dt.
(12)
~=
Its corresponding transform of discrete sequence Mm(k¹/(N!1))N is given by
A B A B
1 N~1 k¹ 2nlk HK (l)" + m cas m N N!1 N k/0 for k, l"0, 1, 2 , N!1,
(13)
where N is the number of samples within the fixed interval [0, ¹ ]. Moreover, the integral for the product of two signals in the form m (t) m(i) (t) is given by 1 2 HM (0,i) (mu )"E (mu )* HM (i) (mu )#O (mu )* HM (i) (!mu ) 0 m1 0 0 m1 0 0 m1 ,m2 m2 m2 "H (mu ) ? HM (i) (mu ), (14) 0 0 m1 m2 where the operator ? symbol represents the two convolutions in short form for simplicity, E (mu )"1 [H (mu )# m1 0 0 2 m1 H (mu )] and O (mu )"1 [H (mu )!H (mu )]. m1 0 m2 0 0 m1 0 2 m1 Now, applying equations (11) and (14) into equation (8) leads to the HMF-model of the BCT system example 3 3 HM (3) (mu )"! + a HM (i~1) (mu )# + b HM (3~i) (mu ) 0 i y 0 i u 0 y i/1 i/1
3 T (i~1) #+ n / (t) u(t) y(t) dt i m 0 i/1 3 j~1 T !+ + n b / (t) u(t) (j~1~i) u(t) dt, j i m 0 j/2 i/1
0
/ (t)m(i)(t) dt"(!1)i m
]cas@
#b u(t)#n u(t) y(t)#n u(t) 3 1 2
n1 n2 + + # (h)! (u,y) P (p) ) (u, y)"0, j jk jk k j/0 k/1
T
AB
1
][y¨ (t)!b uR (t)!b u(t)]. 1 2
P
k k "(!1)i + (!1)j (k#m!j)iui (!1)i 0 j j/0
... y (t)!b u¨ (t)!b uR (t)"!a y(t)!a [yR (t)!b u(t)]
1
(10)
Now, the integrals of equation (8) may be calculated analytically by repeatedly applying integration-by-parts until all the derivatives of the I/O-signals shift to known / (t). For instance, m
Replacing the state variables x (t), i"1, 2, 3 by their corresi ponding I/O-signals and writing in the I/O-differential equation form implies 2
/(m) (t)"0, for t"0 and t"¹. m
(4)
y(t)"[1 0 0] x(t).
1
is a kth-order Hartley modulating function (Patra and Unbehauen, 1995) relative to fixed time interval [0, ¹ ], cas(ut)"cos(ut)#sin(ut), and it has the following three properties:
3 j~1 # + + a b HM (j~1~i) (mu ) j i u 0 j/2 i/1 (8)
k k / (t)" + (!1)j cas((k#m!j)u t), m"0,$1, 2 (9) m 0 j j/0
3 # + n H (mu ) ? HM (i~1) (mu ) i u 0 y 0 i/1 3 j~1 1 !+ + n b HM (j~1~i) (mu ). 0 ( j!1!i) j i u2 j/2 i/1
(15)
Technical Communiques Let z(mu )"HM (3) (mu ), assuming an additional error e(mu ) 0 0 0 y and rearranging terms of equation (15), it can be rewritten as a linear regression, z(mu )"uT(mu )h#e(mu ), 0 0 0
(16)
Letting z(mu )"HM (2) (mu ) and rewriting equation (21) in the y 0 0 form of equation (16) with uT(mu )"[!HM (mu )!HM (1) (mu ) F HM (1) (mu ) HM (mu ) F 0 u 0 u y 0 y 0 0 H ? HM (mu )H ? HM (1) (mu ) F !HM 2 (mu )] u y 0 u y 0 u 0
where and
uT(mu )"[!HM !HM (1)!HM (2) F HM (2) HM (1) HM F 0 y y y u u u
hT"[a a F b b #a b F n n F n b ] 1 2 1 2 2 1 1 2 2 1
H ? HM H ? HM (1) H ? HM (2) F !HM 2!1 HM (1) ] u y u y u y u 2 u2
and
hT"[a a a F b (b #a b ) (b #a b #a b ) F n n n F 1 2 3 1 2 3 1 3 2 1 3 2 1 2 3 (n b #n b )n b ]. 2 1 3 2 3 1 Let a sequence of observations be made for m"0, $1, 2 , $M. Then, (2M#1) regression equations of the form (16) can be represented as a vector equation and minimizing a cost function given by J(h)"1 eT(Mu , h)We(Mu )"1 [z(Mu ) 0 0 2 0 2 !W(Mu )h]T W[z(Mu )!W(Mu )h] 0 0 0
(17)
with respect to the parameter vector h leads to a weighted least squares estimate hK ,hK (Mu )"[(T(Mu )WW(Mu )]~1(T(Mu ) Wz(Mu ), w W 0 0 0 0 0 (18) where W is a positive-definite symmetric frequency dependent weighting matrix, WT(Mu )"[u(!Mu ) 2 u(!u )u(0)u(u ) 2 u(Mu )] 0 0 0 0 0 and zT(Mu )"[z(!Mu ) 2 z(!u )z(0)z(u ) 2 z(Mu )]. 0 0 0 0 0 4. Illustrative example To illustrate the approach and its effectiveness a simulation study for purposes of comparison is performed for a 2nd-order BCT system given by in its state space representation as
C D C
DC D C D C DC D
xR (t) 0 1 " xR (t) !a !a 2 1 2 #u(t)
0 0 n n 1 2
501
(22)
leads to a linear regression. To investigate the methodology u(t)"0.5#0.25Mcos 0.5nt#sin ntN over [0, ¹]"[0, 9 s] is considered. For simulation purposes the I/O-data are generated by solving the BCT system differential equations using Runge—Kutta formulas, and the output signal is also mixed with noise for a Monte-Carlo simulation depending on a chosen value of the noise-to-signal ratio (NSR) (Fig. 1). Having the I/O-data, first compute the HT of the I/O-data using the standard parabolic rule and then the Hartley spectra of the I/O-signals, their derivatives and convolutions. Next, estimate the parameters of the HMF-model applying a weighted LS algorithm (18) with a weighted factor inversely proportional to frequency. Finally, repeat the indicated computations for a required number of Monte Carlo runs (50) and calculate the average value of the parameters. The root mean square (RMS) normalized error criterion defined by E*hE"M1 + nh [(hK !h )/h ]2N1@2 is used to measure the nh i/1 i *536% *536% bias of the estimates (Table 1) and to study the computational considerations (Table 2) associated with the method, where h4 i and h are the estimated and true parameters, respectively, *536% and n is the number of parameters. In the parameter estimation h the record length is discretized into 210 subintervals and only (2M#1)"11 spectral data points are taken. The estimated parameters of the BCT model are shown in Fig. 2a and b. In addition, the values of the estimated BCT model parameters are given in Table 1 for different levels of the NSR. In addition, the RMS normalized error E*hE criterion is applied for different values of N of the HT according to equation (13) for further investigation and to study some of the computational considerations. The results of the simulation study reported in Table 2 show that as N increased from N"256 to 1024 in the noise free or noisy case the accuracy of the estimated parameters also increased which is directly due to the corresponding smaller sampling time selection (t "¹/N). 4!.1
x (t) b 1 # 1 u(t) x (t) b 2 2 x (t) 1 , x (t) 2
(19)
y(t)"[1 0] x(t), where Ma , a , b , b , n , n N"M2, 3, 0.7, 1.2,!1.9, 2.7N. 1 2 1 2 1 2 The corresponding I/O-differential equation results in y¨ (t)"!a y(t)!a yR (t)#b uR (t)#b u(t)#a b u(t) 1 2 1 2 2 1 #n u(t)y(t)#n u(t)yR (t)!n b u2(t). 1 2 2 1
(20)
Now, modulating equation (20) using n"2 as well as applying their corresponding HMF spectrum notations discussed in Section 3 leads to the HMF-model of the form HM (2) (mu )"!a HM (mu )!a HM (1) (mu )#b HM (1) (mu ) y 0 1 y 0 2 y 0 1 u 0 #(b #a b ) HM (mu )#n H ? HM (mu ) 2 2 1 u 0 1 u y 0 #n H ? HM (1) (mu )!n b HM 2 (mu ). 2 u y 0 2 1 u 0
(21)
Fig. 1. I/O-response plots of BCT system example (19) for various levels of NSR of y(t).
502
Technical Communiques
Table 1. Estimated parameter values of the BCT model HMF-model parameters (mean) NSR in % 0 1 2 3 5 10 20 30 40 50 60 70
aL :2 1
aL :3 2
bK :0.7 1
nL :!1.9 nL :2.7 1 2
2.000 2.002 1.997 1.996 2.001 1.973 1.960 2.012 2.081 2.101 1.952 1.958
3.000 3.002 2.999 2.996 3.003 2.978 2.980 3.033 3.142 3.153 3.015 3.096
0.700 0.700 0.700 0.700 0.699 0.697 0.695 0.705 0.700 0.717 0.692 0.705
!1.900 !1.895 !1.908 !1.902 !1.905 !1.959 !2.005 !1.920 !1.900 !1.743 !2.195 !2.173
HMF-models parameters, continued NSR b #Y a b 2 2 1 in % :3.3
n Yb 2 1 :1.89
0 1 2 3 5 10 20 30 40 50 60 70
1.890 1.893 1.885 1.886 1.892 1.853 1.841 1.898 1.960 2.071 1.782 1.793
3.300 3.301 3.298 3.296 3.304 3.282 3.282 3.322 3.385 3.397 3.307 3.305
2.700 2.703 2.697 2.693 2.702 2.667 2.695 2.720 2.969 3.079 2.792 2.906
RMS E*hE
Estimated bK "b #Y a b 2 2 2 1 !aL bK 2 1 :1.2
0.00000 0.00132 0.00202 0.00167 0.00138 0.01595 0.02465 0.00789 0.04757 0.07774 0.06466 0.06614
1.200 1.200 1.199 1.199 1.205 1.206 1.211 1.184 1.186 1.136 1.221 1.122
Table 2. RMS normalized error E*hE for different order N of HT (13) and NSR cases E*hE NSR in %
N"256
N"512
0% 10% 70%
1.375]10~3 2.392]10~2 7.616]10~2
3.627]10~4 2.199]10~2 7.043]10~2
N"1024 0.000000 1.595]10~2 6.614]10~2
Fig. 2. Bilinear state space (19) or its I/O-differential equation (20) parameters estimate, aL "2, aL "3, nL "!1.9, nL "2.7, 1 2 1 2 bK "0.7, (b #Y a b )"3.3, (n Y b )"1.89 and bK "1.2. 2 1 2 1 2 1 2
Furthermore, the plot comparing the actual BCT model output y(t) with its estimate generated from its corresponding highly noise contaminated signal (70% NSR case) is presented in Fig. 3. It is evident from the results presented in Figs. 2 and 3 as well as Tables 1 and 2 that the identification scheme shows its excellent performance and effectiveness for estimation of physically-based parameters of BCT system. 5. Discussions and conclusions Identification of BCT system from recorded I/O-data via a frequency dependent weighted LS algorithm and HMF-method has been presented. Referring the I/O-differential equation of the BCT system example and its corresponding HMF-model it shows that the structure of the dynamics is not influenced by applying HMF-method, however, the I/O-measured signals are replaced by their corresponding spectra. This property of the approach results in its appropriateness towards identifying the system parameters that relate directly the estimated to their physicallybased dynamic representation. The developed procedure and the obtained results may form a basis for identification of a wider class of nonlinear systems. Special attention is also paid to measure the bias of the estimates and to investigate the
Fig. 3. Actual model output y(t) and its estimate from 70% NSR case.
computational considerations associated with the technique using root mean square criterion. A batch scheme recursive identification of the methodology is under study.
Technical Communiques Acknowledgements—The first author thanks the German academic exchange service for providing the financial support for his stay in Ruhr-University, Bochum. He is also thankful to the Department of Electrical Engineering, Faculty of Technology, Addis Ababa University, Addis Ababa, Ethiopia, for sanctioning leave of absence. References Bracewll, R. N. (1986). ¹he Hartley ¹ransform. Oxford University Press, New York. Bruni, C., G. DiPillo and G. Koch (1974). Bilinear systems: an appealing class of nearly linear systems in theory and applications. IEEE ¹rans. Autmat. Control, AC 19(4), 334—348. Cheng, B. and N. S. Hsu (1982). Analysis and parameter estimation systems via block-pulse functions. IEEE Int. J. Control, 36(1), 53—65. Dunoyer, A., L. Balmer, K. J. Burnham and D. J. G. James (1996a). On the characteristics of practical bilinear model structures. Systems Sci., 22 (2), 43—58. Dunoyer, A., L. Balmer, K. J. Burnham and D. J. G. James (1996b). Systems modeling and control for industrial plant: a bilinear approach. In 13th IFAC ¹riennial ¼orld Congress, San Francisco, USA, pp. 347—352. Goodhart, S. G., K. J. Bur Burnham and D. J. G. James (1991).
503
A bilinear self-tuning controller for industrial heating plant. In Proc. IEE Int. Conf. Control, pp. 779—783. Haber, R. and L. Keviczky (1976). Identification of nonlinear dynamic systems. In 4th IFAC Symp. on Ident. and Syst.Par. Est., USSR (1), pp. 62—112. Karanam, V. R., P. A. Frick and R. R. Mohler (1978). Bilinear system identification by Walsh functions. IEEE ¹rans. Automat. Control, AC 23 (4), 709—713. Mohler, R. R. (1973). Bilinear Control Process. Academic Press, New York. Mohler, R. R. (1991). Nonlinear Systems: »ol. II Application to Bilinear Control. Prentice Hall, Englewood Cliffs, NJ. Patra, A. and H. Unbehauen (1995). Identification of a class of nonlinear continuous-time systems using Hartley modulating functions. Int. J. Control, 62(2), 1431—1451. Pearson, A. E. (1992). Explicit parameter identification for a class of nonlinear input/output differential operator models. In Proc. of 31st IEEE CDC, Arizona, pp. 3656—3660. Rao, K. V., P. A. Frick and R. R. Mohler (1976). On bilinear system identification by Walsh functions. In 4th IFAC Symp. on Ident. and Syst. Par. Est., USSR (3), 350—359. Reuter, H. (1994). State space identification of bilinear canonical forms. In Proc. IEE Int. Conf. Control, pp. 833—838. Rink, R. E. and R. R. Mohler (1968). Completely controllable bilinear systems. NSIAM J. Control, 6(3), 477—486.
PII: S0005–1098(97)00216–1
Technical Communique
Bilinear Continuous-Time Systems Identification via Hartley-Based Modulating Functions* S. DANIEL-BERHE-‡ and H. UNBEHAUENKey Words—Bilinear continuous-time systems; identification; Hartley modulating functions method.
Hsu, 1982), etc., have been used to approximate the I/O-sequences of a bilinear system in order to identify the unknown parameters. The present study is focused towards the bilinear continuous-time systems identification using Hartley modulating functions and is organized as follows. In the next part, the structure of bilinear continuous-time systems is discussed. Then identification of the proposed method is presented followed by a typical illustrative example. Finally, the paper finishes with some discussions and conclusions as well as future research plans.
Abstract—This paper highlights the relevance and merits of the Hartley modulating functions (HMF) method for the identification of bilinear continuous-time (BCT) systems from recorded input and noise-contaminated output data and it provides an insight into parameter estimation of a wider range of nonlinear systems in practice. The methodology replaces the I/O-differential equation representing the dynamic system behavior by the Hartley spectrum equation. As a result it involves the known derivatives of the modulating function instead of the derivatives of the input and noisy output data by applying integral transformation to signals. A frequency weighted least-squares algorithm is also applied in the identification and a normalized root mean square criterion is used to investigate some computational considerations and the bias of the estimates. Results of the simulation studies demonstrate the appropriateness of the approach and its efficiency. ( 1998 Elsevier Science Ltd. All rights reserved.
2. Structure of bilinear continuous-time systems The mathematical description of a continuous-time singleinput—single-output (SISO) bilinear system can be given in the following state-space representation: x5 (t)"Ax(t)#bu(t)#u(t)Nx(t), y(t)"cTx(t),
1. Introduction Clearly, most real systems are nonlinear and continuous in time and a linear representation of their actual behavior may lead to limitations in applications. Correspondingly, nonlinear parametric systems may be modeled using special or general nonlinear models. Bilinear systems are one of the special subclass and well-structured nonlinear systems. In the past few decades attentions have been focused on the class of bilinear systems from the theoretical and the practical point of view (Rink and Mohler, 1968; Mohler, 1973, 1991; Bruni et al., 1974; Haber and Keviczky, 1976, etc.). They arise naturally in a wide variety of applications and are found to be quite suitable for fitting some important actual processes, for instance, nuclear fission, hydraulic drive, turbo-pump dynamics, convection heat transfer, etc., have been modeled approximately by bilinear systems (Mohler, 1973, 1991; Goodhart et al., 1991; Reuter, 1994, etc.). In fact, it appears that many important processes, not only in engineering, but also in biology, socio-economics, ecology, etc., may be modeled by bilinear systems, refer Mohler (1973, 1991), etc. Recently, a bilinear approach for system modeling and control of an industrial plant as well as the characteristics of practical bilinear model structures have been presented by Dunoyer et al. (1996a, b), respectively. Furthermore, the Walsh functions (Rao et al., 1976; Karanam et al., 1978), the block-pulse functions (Cheng and
(1)
where u(t) and y(t) are the measured scalar control input and output, respectively, x(t) is an n-dimensional state vector, A and N are the n]n, matrices of real constants, b and c are the n]1 vectors of real constants. Bilinear systems can be formulated in different forms, however, for the detail mathematical analysis of this presentation it is assumed to be represented in the phase variable canonical form of the following matrices of A, N, b and c: 0 1 . . . 0 . . . . . . . . A" , . . . . 0 . . . 0 1 !a . . . !a !a 1 n~1 n 0 . . N" . 0 n 1
* Received 13 May 1997; received in final form 10 November 1997. This paper was recommended for publication in revised form by Editor P. Dorato. Corresponding author Prof. Dr.-Ing. H. Unbehauen. Tel. #49/234 700 4071; Fax #49/ 234 7094 101; E-mail [email protected]. - Control Engineering Laboratory, Faculty of Electrical Engineering, Ruhr-University Bochum, Building IC 3/131, D-44780 Bochum, Germany. ‡ On German Academic Exchange Service Research Fellowship and on leave from the Department of Electrical Engineering, Addis Ababa University, Addis Ababa, Ethiopia.
. . .
0 . . . . . , . . . . 0 . . . n n
b 1 . . b" . b n~1 b n
(2) 1 . . and C" . . 0 0
Representation of equations (1) and (2) can be rewritten in the following simplified form: xR (t)"x (t)#b u(t), j j`1 j
j"1, 2, 2 , n!1,
n xR (t)" + [!a #n u(t)] x (t)#b u(t), n i i i n i/1 y(t)"x (t). 1 499
(3)
500
Technical Communiques
Now, consider a third-order BCT observability canonical form given by equation (4) and rewrite in measured I/O-variables, i.e. xR (t) 0 1 0 1 xR (t) " 0 0 1 2 xR (t) !a !a !a 3 1 2 3
x (t) 1 x (t) 2 x (t) 3
b 0 0 0 1 # b u(t)#u(t) 0 0 0 2 b n n n 3 1 2 3
G
/ (t), 04t4¹ / (t)" m , m 0 otherwise x (t) 1 x (t) , 2 x (t) 3
/(l) (t) exists for all l"0, 1, 2 , k!1, and m
2
!a [y¨ (t)!b uR (t)!b u(t)] 3 1 2
][yR (t)!b u(t)]#n u(t) 1 3 (5)
Rearranging the terms leads to the preceding equation of the form 3 3 3 j~1 (3) y(t)"! + a #(i~1) y(t) # + b (3~i) u(t) # + + a b (j~1~i) u(t) i i j i i/1 i/1 j/2 i/1 (6) 3 3 j~1 # + n u(t) (i~1) y(t) ! + + n b u(t) (j~1~i) u(t) , i j i i/1 j/2 i/1 (0) where l(t)"(d0/dt0) l(t)"l(t). The above derivation can be generalized for an nth order of the form given in equation (2). The main purpose is to estimate the BCT system parameters a , i b and n , i"1, 2, 3 from the input—output data records. Generi i ally, equation (6) is a standard form of a nonlinear input—output differential operator model (Pearson, 1992) of the form (7)
where # (h) are given functions of the parameter vector h with j # "1, ! (u, y) and ) (u, y) are specified functions of the I/O0 jk k signals, P (p) are fixed polynomials of degree n in the differential jk operator p"d/dt, j"0, 1, 2 , n , and k"0, 1, 2 , n . 1 2 3. Identification In general, the modulating functions method starts with multiplying both sides of equation (6) by / (t) and integrating over m the length of the modulating function which results in
P
P
T
P
3 T 3 T (3) / (t) y(t) dt"!+ a / (t) (i~1) y(t) dt# + b / (t) (3~i) u(t) dt m i m i m 0 0 0 i/1 i/1
P
3 j~1 T #+ + ab / (t) (j~1~i) u(t) dt j i m 0 j/2 i/1
P
where
P
AB
P
T
0
m(t) /(i) (t) dt m
AB AB
in :T m(t) cas((!1)i(k#m!j) u (t) dt 0 0hgggggigggggj 2 i Hm((~1) (k`m~j)u0) k k in " + (!1)j (k#m!j)iui cas@ 0 j 2 j/0
AB
]H ((!1)i (k#m!j)) u , m 0 "HM (i) (mu ) 0 m
(11)
where HM (i) (mu ) is the mth HMF spectral component of the 0 m ith derivative of a continuous signal m(t), and H (mu ) is the m 0 Hartley transform (HT) (Bracewll, 1986) of m(t) defined by
P
H (u)" m
=
m(t) cas ut dt.
(12)
~=
Its corresponding transform of discrete sequence Mm(k¹/(N!1))N is given by
A B A B
1 N~1 k¹ 2nlk HK (l)" + m cas m N N!1 N k/0 for k, l"0, 1, 2 , N!1,
(13)
where N is the number of samples within the fixed interval [0, ¹ ]. Moreover, the integral for the product of two signals in the form m (t) m(i) (t) is given by 1 2 HM (0,i) (mu )"E (mu )* HM (i) (mu )#O (mu )* HM (i) (!mu ) 0 m1 0 0 m1 0 0 m1 ,m2 m2 m2 "H (mu ) ? HM (i) (mu ), (14) 0 0 m1 m2 where the operator ? symbol represents the two convolutions in short form for simplicity, E (mu )"1 [H (mu )# m1 0 0 2 m1 H (mu )] and O (mu )"1 [H (mu )!H (mu )]. m1 0 m2 0 0 m1 0 2 m1 Now, applying equations (11) and (14) into equation (8) leads to the HMF-model of the BCT system example 3 3 HM (3) (mu )"! + a HM (i~1) (mu )# + b HM (3~i) (mu ) 0 i y 0 i u 0 y i/1 i/1
3 T (i~1) #+ n / (t) u(t) y(t) dt i m 0 i/1 3 j~1 T !+ + n b / (t) u(t) (j~1~i) u(t) dt, j i m 0 j/2 i/1
0
/ (t)m(i)(t) dt"(!1)i m
]cas@
#b u(t)#n u(t) y(t)#n u(t) 3 1 2
n1 n2 + + # (h)! (u,y) P (p) ) (u, y)"0, j jk jk k j/0 k/1
T
AB
1
][y¨ (t)!b uR (t)!b u(t)]. 1 2
P
k k "(!1)i + (!1)j (k#m!j)iui (!1)i 0 j j/0
... y (t)!b u¨ (t)!b uR (t)"!a y(t)!a [yR (t)!b u(t)]
1
(10)
Now, the integrals of equation (8) may be calculated analytically by repeatedly applying integration-by-parts until all the derivatives of the I/O-signals shift to known / (t). For instance, m
Replacing the state variables x (t), i"1, 2, 3 by their corresi ponding I/O-signals and writing in the I/O-differential equation form implies 2
/(m) (t)"0, for t"0 and t"¹. m
(4)
y(t)"[1 0 0] x(t).
1
is a kth-order Hartley modulating function (Patra and Unbehauen, 1995) relative to fixed time interval [0, ¹ ], cas(ut)"cos(ut)#sin(ut), and it has the following three properties:
3 j~1 # + + a b HM (j~1~i) (mu ) j i u 0 j/2 i/1 (8)
k k / (t)" + (!1)j cas((k#m!j)u t), m"0,$1, 2 (9) m 0 j j/0
3 # + n H (mu ) ? HM (i~1) (mu ) i u 0 y 0 i/1 3 j~1 1 !+ + n b HM (j~1~i) (mu ). 0 ( j!1!i) j i u2 j/2 i/1
(15)
Technical Communiques Let z(mu )"HM (3) (mu ), assuming an additional error e(mu ) 0 0 0 y and rearranging terms of equation (15), it can be rewritten as a linear regression, z(mu )"uT(mu )h#e(mu ), 0 0 0
(16)
Letting z(mu )"HM (2) (mu ) and rewriting equation (21) in the y 0 0 form of equation (16) with uT(mu )"[!HM (mu )!HM (1) (mu ) F HM (1) (mu ) HM (mu ) F 0 u 0 u y 0 y 0 0 H ? HM (mu )H ? HM (1) (mu ) F !HM 2 (mu )] u y 0 u y 0 u 0
where and
uT(mu )"[!HM !HM (1)!HM (2) F HM (2) HM (1) HM F 0 y y y u u u
hT"[a a F b b #a b F n n F n b ] 1 2 1 2 2 1 1 2 2 1
H ? HM H ? HM (1) H ? HM (2) F !HM 2!1 HM (1) ] u y u y u y u 2 u2
and
hT"[a a a F b (b #a b ) (b #a b #a b ) F n n n F 1 2 3 1 2 3 1 3 2 1 3 2 1 2 3 (n b #n b )n b ]. 2 1 3 2 3 1 Let a sequence of observations be made for m"0, $1, 2 , $M. Then, (2M#1) regression equations of the form (16) can be represented as a vector equation and minimizing a cost function given by J(h)"1 eT(Mu , h)We(Mu )"1 [z(Mu ) 0 0 2 0 2 !W(Mu )h]T W[z(Mu )!W(Mu )h] 0 0 0
(17)
with respect to the parameter vector h leads to a weighted least squares estimate hK ,hK (Mu )"[(T(Mu )WW(Mu )]~1(T(Mu ) Wz(Mu ), w W 0 0 0 0 0 (18) where W is a positive-definite symmetric frequency dependent weighting matrix, WT(Mu )"[u(!Mu ) 2 u(!u )u(0)u(u ) 2 u(Mu )] 0 0 0 0 0 and zT(Mu )"[z(!Mu ) 2 z(!u )z(0)z(u ) 2 z(Mu )]. 0 0 0 0 0 4. Illustrative example To illustrate the approach and its effectiveness a simulation study for purposes of comparison is performed for a 2nd-order BCT system given by in its state space representation as
C D C
DC D C D C DC D
xR (t) 0 1 " xR (t) !a !a 2 1 2 #u(t)
0 0 n n 1 2
501
(22)
leads to a linear regression. To investigate the methodology u(t)"0.5#0.25Mcos 0.5nt#sin ntN over [0, ¹]"[0, 9 s] is considered. For simulation purposes the I/O-data are generated by solving the BCT system differential equations using Runge—Kutta formulas, and the output signal is also mixed with noise for a Monte-Carlo simulation depending on a chosen value of the noise-to-signal ratio (NSR) (Fig. 1). Having the I/O-data, first compute the HT of the I/O-data using the standard parabolic rule and then the Hartley spectra of the I/O-signals, their derivatives and convolutions. Next, estimate the parameters of the HMF-model applying a weighted LS algorithm (18) with a weighted factor inversely proportional to frequency. Finally, repeat the indicated computations for a required number of Monte Carlo runs (50) and calculate the average value of the parameters. The root mean square (RMS) normalized error criterion defined by E*hE"M1 + nh [(hK !h )/h ]2N1@2 is used to measure the nh i/1 i *536% *536% bias of the estimates (Table 1) and to study the computational considerations (Table 2) associated with the method, where h4 i and h are the estimated and true parameters, respectively, *536% and n is the number of parameters. In the parameter estimation h the record length is discretized into 210 subintervals and only (2M#1)"11 spectral data points are taken. The estimated parameters of the BCT model are shown in Fig. 2a and b. In addition, the values of the estimated BCT model parameters are given in Table 1 for different levels of the NSR. In addition, the RMS normalized error E*hE criterion is applied for different values of N of the HT according to equation (13) for further investigation and to study some of the computational considerations. The results of the simulation study reported in Table 2 show that as N increased from N"256 to 1024 in the noise free or noisy case the accuracy of the estimated parameters also increased which is directly due to the corresponding smaller sampling time selection (t "¹/N). 4!.1
x (t) b 1 # 1 u(t) x (t) b 2 2 x (t) 1 , x (t) 2
(19)
y(t)"[1 0] x(t), where Ma , a , b , b , n , n N"M2, 3, 0.7, 1.2,!1.9, 2.7N. 1 2 1 2 1 2 The corresponding I/O-differential equation results in y¨ (t)"!a y(t)!a yR (t)#b uR (t)#b u(t)#a b u(t) 1 2 1 2 2 1 #n u(t)y(t)#n u(t)yR (t)!n b u2(t). 1 2 2 1
(20)
Now, modulating equation (20) using n"2 as well as applying their corresponding HMF spectrum notations discussed in Section 3 leads to the HMF-model of the form HM (2) (mu )"!a HM (mu )!a HM (1) (mu )#b HM (1) (mu ) y 0 1 y 0 2 y 0 1 u 0 #(b #a b ) HM (mu )#n H ? HM (mu ) 2 2 1 u 0 1 u y 0 #n H ? HM (1) (mu )!n b HM 2 (mu ). 2 u y 0 2 1 u 0
(21)
Fig. 1. I/O-response plots of BCT system example (19) for various levels of NSR of y(t).
502
Technical Communiques
Table 1. Estimated parameter values of the BCT model HMF-model parameters (mean) NSR in % 0 1 2 3 5 10 20 30 40 50 60 70
aL :2 1
aL :3 2
bK :0.7 1
nL :!1.9 nL :2.7 1 2
2.000 2.002 1.997 1.996 2.001 1.973 1.960 2.012 2.081 2.101 1.952 1.958
3.000 3.002 2.999 2.996 3.003 2.978 2.980 3.033 3.142 3.153 3.015 3.096
0.700 0.700 0.700 0.700 0.699 0.697 0.695 0.705 0.700 0.717 0.692 0.705
!1.900 !1.895 !1.908 !1.902 !1.905 !1.959 !2.005 !1.920 !1.900 !1.743 !2.195 !2.173
HMF-models parameters, continued NSR b #Y a b 2 2 1 in % :3.3
n Yb 2 1 :1.89
0 1 2 3 5 10 20 30 40 50 60 70
1.890 1.893 1.885 1.886 1.892 1.853 1.841 1.898 1.960 2.071 1.782 1.793
3.300 3.301 3.298 3.296 3.304 3.282 3.282 3.322 3.385 3.397 3.307 3.305
2.700 2.703 2.697 2.693 2.702 2.667 2.695 2.720 2.969 3.079 2.792 2.906
RMS E*hE
Estimated bK "b #Y a b 2 2 2 1 !aL bK 2 1 :1.2
0.00000 0.00132 0.00202 0.00167 0.00138 0.01595 0.02465 0.00789 0.04757 0.07774 0.06466 0.06614
1.200 1.200 1.199 1.199 1.205 1.206 1.211 1.184 1.186 1.136 1.221 1.122
Table 2. RMS normalized error E*hE for different order N of HT (13) and NSR cases E*hE NSR in %
N"256
N"512
0% 10% 70%
1.375]10~3 2.392]10~2 7.616]10~2
3.627]10~4 2.199]10~2 7.043]10~2
N"1024 0.000000 1.595]10~2 6.614]10~2
Fig. 2. Bilinear state space (19) or its I/O-differential equation (20) parameters estimate, aL "2, aL "3, nL "!1.9, nL "2.7, 1 2 1 2 bK "0.7, (b #Y a b )"3.3, (n Y b )"1.89 and bK "1.2. 2 1 2 1 2 1 2
Furthermore, the plot comparing the actual BCT model output y(t) with its estimate generated from its corresponding highly noise contaminated signal (70% NSR case) is presented in Fig. 3. It is evident from the results presented in Figs. 2 and 3 as well as Tables 1 and 2 that the identification scheme shows its excellent performance and effectiveness for estimation of physically-based parameters of BCT system. 5. Discussions and conclusions Identification of BCT system from recorded I/O-data via a frequency dependent weighted LS algorithm and HMF-method has been presented. Referring the I/O-differential equation of the BCT system example and its corresponding HMF-model it shows that the structure of the dynamics is not influenced by applying HMF-method, however, the I/O-measured signals are replaced by their corresponding spectra. This property of the approach results in its appropriateness towards identifying the system parameters that relate directly the estimated to their physicallybased dynamic representation. The developed procedure and the obtained results may form a basis for identification of a wider class of nonlinear systems. Special attention is also paid to measure the bias of the estimates and to investigate the
Fig. 3. Actual model output y(t) and its estimate from 70% NSR case.
computational considerations associated with the technique using root mean square criterion. A batch scheme recursive identification of the methodology is under study.
Technical Communiques Acknowledgements—The first author thanks the German academic exchange service for providing the financial support for his stay in Ruhr-University, Bochum. He is also thankful to the Department of Electrical Engineering, Faculty of Technology, Addis Ababa University, Addis Ababa, Ethiopia, for sanctioning leave of absence. References Bracewll, R. N. (1986). ¹he Hartley ¹ransform. Oxford University Press, New York. Bruni, C., G. DiPillo and G. Koch (1974). Bilinear systems: an appealing class of nearly linear systems in theory and applications. IEEE ¹rans. Autmat. Control, AC 19(4), 334—348. Cheng, B. and N. S. Hsu (1982). Analysis and parameter estimation systems via block-pulse functions. IEEE Int. J. Control, 36(1), 53—65. Dunoyer, A., L. Balmer, K. J. Burnham and D. J. G. James (1996a). On the characteristics of practical bilinear model structures. Systems Sci., 22 (2), 43—58. Dunoyer, A., L. Balmer, K. J. Burnham and D. J. G. James (1996b). Systems modeling and control for industrial plant: a bilinear approach. In 13th IFAC ¹riennial ¼orld Congress, San Francisco, USA, pp. 347—352. Goodhart, S. G., K. J. Bur Burnham and D. J. G. James (1991).
503
A bilinear self-tuning controller for industrial heating plant. In Proc. IEE Int. Conf. Control, pp. 779—783. Haber, R. and L. Keviczky (1976). Identification of nonlinear dynamic systems. In 4th IFAC Symp. on Ident. and Syst.Par. Est., USSR (1), pp. 62—112. Karanam, V. R., P. A. Frick and R. R. Mohler (1978). Bilinear system identification by Walsh functions. IEEE ¹rans. Automat. Control, AC 23 (4), 709—713. Mohler, R. R. (1973). Bilinear Control Process. Academic Press, New York. Mohler, R. R. (1991). Nonlinear Systems: »ol. II Application to Bilinear Control. Prentice Hall, Englewood Cliffs, NJ. Patra, A. and H. Unbehauen (1995). Identification of a class of nonlinear continuous-time systems using Hartley modulating functions. Int. J. Control, 62(2), 1431—1451. Pearson, A. E. (1992). Explicit parameter identification for a class of nonlinear input/output differential operator models. In Proc. of 31st IEEE CDC, Arizona, pp. 3656—3660. Rao, K. V., P. A. Frick and R. R. Mohler (1976). On bilinear system identification by Walsh functions. In 4th IFAC Symp. on Ident. and Syst. Par. Est., USSR (3), 350—359. Reuter, H. (1994). State space identification of bilinear canonical forms. In Proc. IEE Int. Conf. Control, pp. 833—838. Rink, R. E. and R. R. Mohler (1968). Completely controllable bilinear systems. NSIAM J. Control, 6(3), 477—486.