- Email: [email protected]
Experimental physical parameter estimation of a thyristor driven DC-motor using the HMF-method
Control Engineering Practice 6 (1998) 615—626
Experimental physical parameter estimation of a thyristor driven DC-motor using the HMF-method1 S. Daniel-Berhe, H. Unbehauen* Control Engineering Laboratory, Faculty of Electrical Engineering, Ruhr-University Bochum, ESR IC 3/131, D-44780 Bochum, Germany Received 28 April 1997; accepted 17 January 1998
Abstract This paper deals with the physical parameter estimation of the nonlinear continuous-time dynamics of a thyristor-driven dc-motor experimental set-up using the Hartley modulating functions (HMF) method. The approach applies a frequency-weighted leastsquares formulation to estimate the nonlinear continuous-time system parameters based on input-output data records over a finite time interval. The laboratory experimental set-up of the plant is facilitated by a computer-aided control system design software system, CADACS. The method has given satisfactory experimental results, and demonstrates its feasibility for use in parameter estimation of physically-based continuous-time systems. Due to its encouraging results, further studies will continue, in order to implement the methodology in the form of batch scheme recursive estimation for nonlinear continuous-time dynamics. ( 1998 Elsevier Science ¸td. All rights reserved. Keywords: Physical parameter estimation; dc motors; modulating functions; nonlinear; continuous-time systems
1. Introduction Most real, dynamic systems are continuous in time and nonlinear, and can be represented by a set of differential equations in which the parameters are physically defined. For numerical treatment, the dynamics of real systems are usually described by discrete-time models, for instance in the form of difference equations. However, a number of practical and theoretical considerations indicate that continuous-time treatment of these systems is often relevant, leading to continuous-time parameterestimation methods. Furthermore, identification of nonlinear continuous-time systems from available sampled input-output data records is desirable, to enhance the interpretation of the nonlinear dynamic system behavior, and to support the controller design.
* All correspondence should be sent to Professor Dr.-Ing. Heinz Unbehauen. Tel.: ##49 234 700 4071; fax ##49 234 709-4101; e-mail: [email protected] 1 The original version of this paper was presented at the 11th IFAC Symposium on System Identification, which was held in Fukuoka, Japan during 8—11 July 1997.
In recent years the identification of continuous-time linear and nonlinear systems is a topic that has received considerable attention. For instance, the relevance and importance of continuous-time system identification, with detailed reviews, have been reported in the literature by Young (1981), Unbehauen and Rao (1987, 1990). Sinha and Rao (1991), etc. These works discussed the relationship of the estimated parameters to their physically-based continuous-time system descriptions, and reviewed the progress of research on parameter estimation for continuous-time models of dynamic systems over the past few decades. Furthermore, the need for models/methods for nonlinear continuous-time systems identification is due to the fact that a conventional continuous-time linear model holds in only a small neighborhood of the actual equilibrium state, and leads to limitations in applications. Today, quite a number of methods are available for the parameter identification of linear continuous-time systems. Unfortunately these methods cannot be applied for the identification of nonlinear continuous-time systems. The use of standard linear techniques has been mostly restricted to only very special classes of nonlinear systems, such as bilinear systems, Hammerstein systems, etc.
0967-0661/98/$ — See front matter ( 1998 Elsevier Science Ltd. All rights reserved PII S 0 9 6 7 - 0 6 6 1 ( 9 8 ) 0 0 0 3 6 - 7
616
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
In general, some of the familiar approaches of continuous-time systems identification include Poisson moment functionals (Saha et al., 1991), block pulse functions (Kung and Shih, 1986), Walsh functions (Karanam et al., 1978), orthogonal polynomials (Wang and Marteau, 1987), general hybrid orthogonal functions (Patra, 1989), Laguerre polynomials (Wang et al., 1984), time moments (Subrahmanyam et al., 1995), delayed state variable filters (Tsang and Billings, 1994), Shinbrot’s method of moment functional (Shinbrot, 1957), Hermite modulating functions (Jalili et al., 1992), trigonometric modulating functions (Pearson and Lee, 1983), Fourier modulating functions (Pearson and Lee, 1985; Pearson, 1992), splinetype modulating functions (Maletinsky, 1979), Hartley modulating functions methods (Patra and Unbehauen, 1995). Most of the aforementioned methods were first developed for linear continuous-time systems, and further developed for special nonlinear systems, and have shown promising results. Generally, they require further studies for nonlinear continuous-time system identification. Relatively, the literature of moduling functions methods shows very promising results for the identification of a large class of nonlinear continuous-time systems (Unbehauen and Rao, 1987). The concept of modulating functions, introduced by Shinbrot in 1957, is commonly known as Shinbrot’s method of moment functionals. Essentially, the use of modulating functions involves converting a differential expression involving input-output signals on a specified time interval, into a sequence of algebraic equations. Shinbrot explored parameter identification of higher-order nonlinear dynamic systems by applying integral transformation techniques. In particular, Fourier modulating functions methods have shown encouraging performance for the identification of linear, integrable as well as convolvable, nonlinear continuoustime models (Pearson and Lee 1985, Pearson 1992). Basically, the Hartley modulating functions (Patra and Unbehauen, 1995) is closely related to the Fourier modulating functions. Compared to the very efficient Fourier modulating functions method, the HMF method has additional important advantages: the Hartley modulating functions are real-valued, and their Hartley spectra can be computed efficiently with the help of fast algorithms. This new methodology is applicable to a large class of nonlinear continuous-time systems by defining a set of HMF for characterizing the continuous process signals. The approach has adopted a class of nonlinear differential operator models, essentially developed for Fourier modulating functions (Pearson, 1992), which has been applied to parameter identification by formulating an efficient frequency-weighted least-squares algorithm. It seems that modulating functions methods provide the most promising alternative for use in parameter estimation of physically-based nonlinear continuous-time
systems (Unbehauen and Rao, 1987). Therefore, in this paper a new, feasible continuous-time parameter estimation method, based on the Hartley modulating functions approach is applied, which directly relates the estimated parameters to their physically-based continuous-time systems representation. Basically, the demonstration of this experimental set-up is intended to show the feasibility of the modulating functions method for physical parameter estimation of nonlinear continuoustime systems. This paper is organized as follows. In Section 2, the background of the Hartley modulating functions, and the Hartley transform and spectra are discussed. The experimental set-up and plant description are outlined in Section 3. A new parameter-estimation algorithm is discussed in Section 4. Experimental results are documented in Section 5. Finally, in Section 6, conclusions and future research plans are presented.
2. Background The general idea of modulating functions was first introduced by Shinbrot (1957). The application of modulating functions to convert a differential equation into an exact algebraic form is well known in the field of system identification, depending on several choices of modulating functions, for instance Shinbrot’s method of moment functionals (Shinbrot, 1957), spline-type modulating functions (Maletinsky, 1979), Hermite modulating functions (Jalili et al., 1992), trigonometric modulating functions (Pearson and Lee, 1983), Fourier modulating functions (Pearson and Lee, 1985; Pearson, 1992), Hartley modulating functions methods (Patra and Unbehauen, 1995), etc. The Hartley modulating functions have recently been introduced as a new member of the modulating functions family. This modulating function is closely related to the Fourier modulating functions. A member of the family of an nth-order HMF / (t), m"0, $1, $2,2 relative to m a fixed time interval [0, ¹ ] is given by
AB
n n / (t)" + (!1)i cas(n#m!i)u t, m 0 i i/0
(1)
where m is the modulating frequency index, cas ut"cos ut#sin ut and u "2n/¹ is essentially the 0 resolving frequency. Furthermore, / (t) is an HMF of m order n on a fixed interval [0, ¹] if it is sufficiently smoother and satisfies for each m"0, $1, $2,2, $M the two-point boundary conditions dk / (t)"0, dtk m
for t"0 and t"¹, k"0, 1,2, n!1
as depicted in Fig. 1.
(2)
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
617
The approach uses Hartley transforms and spectra. The continuous and discrete Hartley transforms (HT) (Hartley, 1942) of a signal x are given by =
P~= x(t)cas ut dt
H (u)"HMx(t)N" x
(3)
and
G A BH
A B A B
k¹ 2nlk 1 N~1 " + x cas (4) N N N k/0 respectively, for l"0, 1,2, N!1, where N is the number of samples within the fixed interval [0, ¹]. The mth HMF spectral component of a signal x(t) is given by k¹ HK (l)"H x x N
T
P0
HM (mu )" x 0
Fig. 1. The family of 2nd-order HMF, / (t), m"0, 1, 2; n"2 and m t3[0, ¹ ].
x(t)/ (t) dt m
AB AB
n :T x(t)cas(n#m!j)u t dt n 0 0 " + (!1)j j hggggiggggj j/0 Hx ((n`m~j)u0) n n H ((n#m!j)u ). (5) " + (!1)j x 0 j j/0 If x(i)(t), i"1, 2,2, n is the ith derivative of the signal x(t), then its corresponding HMF spectrum is given by
ristor-driven dc-motor system and analog input-output display instruments, as shown in Fig. 2a. The armature terminal and field motor circuit input voltages are set from a PC, using computer-aided design and analysis of control systems (CADACS) software (Unbehauen, 1993; Schmid, 1992) through D/A converters. The software has a capability to link the results of a computer-aided analysis and design procedure with the real-time experimental environment. Consequently, the input-output responses can be displayed on the PC screen via converters and interfacing circuits, as well as on the analog instruments. Fig. 2c shows the standard model of a dc motor. The main aim is to apply the measured input-output data to a suitable HMF algorithm of the dc-motor model, and estimate the physically-based system parameters.
AB
n n HM (i) (mu )" + (!1)j cas@(in/2)(n#m!j )i x 0 j j/0 ]ui H ((!1)i(n#m!j)u ), (6) 0 x 0 where cas@(t)"cos(t)!sin(t). Moreover, if the signals x (t) and x (t) are in the form of x (t)x(i) (t), then the 1 2 1 2 HMF spectrum for such a product is given by
G
(l#m#n!j)iH ((!1)i(l#m#n!j)u ) x2 0 n in = n 1 !(!l!m!n#j)iHx2 ((!1)i(!l!m!n#j)u0) HM (0,i) (mu )" + H (lu ) + (!1)j ui cas@ ] x1,x2 0 x1 0 0 j 2 2 #(l!m!n#j)iH ((!1)i(l!m!n#j)u ) l/~= j/0 x2 0 #(!l#m#n!j)iH ((!1)i(!l#m#n!j)u ) x2 0 "H (mu ) ? HM (i)(mu ), x1 0 x2 0
AB
AB
where the operator ? symbol represents the description in short form for simplicity. The proof of (7) is given in the appendix.
3. Experimental set-up and plant description The configuration of the experimental thyristor-driven dc-motor system is shown in Fig. 2(a—c), along with a detailed circuit of the dc motor and load. Primarily, the experimental set-up of the plant consists of a personal computer (PC) with real-time operating system software and interfacing circuits, digital-to-analog (D/A) and analog-to-digital (A/D) converters, a thy-
H
(7)
In general, dc-motor systems have extensive applications, varying from industrial drivers to medical and other instrumentation controls, because of their versatile speed-load characteristics. A dc-motor system can be represented mathematically by the following differential equations: di (t) f u (t)"R i (t)#¸ f ff f dt di (t) a #k i (t)u(t) u (t)"R i (t)#¸ a aa a dt mf J
du(t) "k i (t)i (t)!q (t)!Bu(t) mf a l dt
618
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Fig. 3. Block diagram of the dc motor and load experimental set-up.
Fig. 2(a—c). Laboratory experimental thyristor-driven dc motor system layout. (a) Scheme of the experimental set-up under study. (b) Part of the thyristor-driven dc motor installation. (c) Equivalent circuit representation of the dc motor and load.
q (t)"k t(t)i (t)"k k i (t)i (t)"k i (t)i (t) m 1 a 1 2f a mf a u (t)"k t(t)u(t)"k k i (t)u(t)"k i (t)u(t) (8) i 1 1 2f mf where u (t) is the input voltage for field motor circuit [V], & i (t), the field motor current [A], R , the resistance of the & & field motor circuit [)], ¸ , inductance of the field motor & circuit [H], u (t), armature terminal voltage source [V], ! i (t), armature current [A], R , armature resistance [)], ! ! ¸ , armature inductance [H], u (t), induced voltage in the ! * armature coil [V], u(t), motor angular speed [rad/sec], q (t), electromagnetic torque of the motor [Nm], q (t), . l load torque [Nm], J , J and J, the load, motor and total l . system moment of inertia [N m sec2], B , B and B, the l . load, motor and total system viscous friction constant [N m sec], t(t)"k i (t), the magnetic flux [N m/A], and 2& k "k k , is known as a motor constant [N m/A2]. . 1 2 In general, a dc motor has three inputs, namely, the voltage applied to the armature coil u (t), the voltage ! applied to the field coil u (t), and the load torque q (t) that & l results from the (assumed) machinery that is being driven by the motor. Accordingly, the system of equations in Eq. (8) represents a model of the dc-motor system as a whole, since it consists of the two electrical subsystems describing the field and armature, the mechanical subsystem describing the inertia and friction of the rotating mass, plus the two coupling equations that connect the mechanical subsystem to the electrical subsystem. Moreover, the dc-motor system can be modelled by the simulation diagram shown in Fig. 3. The problem is to estimate the physical parameters of the plant, knowing its structure and based on the inputoutput data records over a finite time interval. In the experimental study a constant load torque and armature terminal voltage sources are considered. Thus, letting the input field current i (t)"u(t), the armature & current i (t)"x (t) and motor angular speed u(t)" ! 1 x (t), then the system equations become 2 R k u xR (t)"! a x (t)! m u(t)x (t)# a 1 2 ¸ 1 ¸ ¸ a a a B k q xR (t)"! x (t)# m u(t)x (t)! l (9) 2 1 J 2 J J y(t)"x (t) 1
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
or
G
yR (t)"!h y(t)!h u(t)x (t)#h u , 1 2 2 3 a (10) xR (t)"!h x (t)#h u(t) y(t)!h q 2 4 2 5 6 l where h "R /¸ , h "k /¸ , h "1/¸ , h "B/J, 1 a a 2 m a 3 a 4 h "k /J and h "1/J. 5 m 6 Eliminating the motor angular speed x (t) and sim2 plifying, the preceding equation in the input field current u(t) and output armature current y(t) yields
u(t) y¨ (t)!uR (t) yR (t)#(h #h ) u(t) yR (t)!h y(t) uR (t) 1 4 1 #h h u3(t)y(t)#h h u(t) y(t)!q h h u2(t) 2 5 1 4 1 2 6 #u h uR (t)!u h h u(t)"0. (11) a 3 a 3 4 Applying the identity uR (t) yR (t)"1 [d2/dt2Mu(t) y(t)N! 2 y(t) u( (t)!u(t) y( (t)] in Eq. (11) and normalizing its first element yields to the standard form of a nonlinear inputoutput differential operator model (Pearson, 1992) of the form u(t) p2y(t)!1 p2Mu(t) y(t)N#1 y(t) p2u(t) 3 3 #2 (h #h ) u(t) py(t)!2 h y(t) pu(t) 3 1 4 3 1 #2 h h u3(t) y(t)#2 h h u(t) y(t)!2 q h h u2(t) 3 2 5 3 1 4 3 l 2 6 #2 u h pu(t)!2 u h h u(t)"0 (12) 3 a 3 3 a 3 4 or m
1
m
2
+ + g (h) F (u, y)P (p) E (u, y)"0 j jk jk k j/0 k/1
(13)
619
Table 1 The functions g (h), F (u, y), P (p) and E (u, y) for the nonlinear model j jk jk k of Eq. (12) j
k
g (h) j
F (u, y) jk
P (p) jk
E (u, y) k
0 1
1 1 2 1 1 1 1 1 1 1
1 1
u(t) !1 3 1 y(t) 3 2 u(t) 3 !2 y(t) 3 2 u3(t) y(t) 3 2 u(t) y(t) 3 !2 q u2(t) 3 l 2u 3 a !2 u u(t) 3 a
p2 p2 p2 p1 p1 p0 p0 p0 p1 p0
y(t) u(t) y(t) u(t) y(t) u(t) 1 1 1 u(t) 1
2 3 4 5 6 7 8
h #h 1 4 h 1 h h 2 5 h h 1 4 h h 2 6 h 3 h h 3 4
and applying their corresponding HMF spectrum notations discussed above, leads to the HMF-model of the system H ? HM (2)"1 H ? H (2)!1 H ? HM (2) 3 u u y 3 y y u !2 (h #h )H ? HM (1)#2 h H ? HM (1) 3 1 y 3 1 4 u y u !2 h h H ? HM 3!2 h h H ? HM u 3 1 4 u y 3 2 5 y #2 q h h HM 2!2 u h HM (1)#2 u h h HM , 3 l 2 6 u 3 a 3 u 3 a 3 4 u (14) q where H ? H (2) represents : Mu(t)® y(t)N/ (t) dt, and apu m y 0 plying integration-by-parts replacing one of the signal by its infinite Hartley series and using the identity for the product of cas function implies
H ? H (2)" u y n 2n H ((l#m#n!j)u )!H ((!l!m!n#j)u ) = n y 0 y 0 (n#m!j)2u2 cas@ ]1 . + H (lu ) + (!1)j 2 #H ((l!m!n#j)u )#H ((!l#m#n!j)u ) 0 u 0 j 2 y 0 y 0 l/~= j/0 (15)
AB
A B G
where g (h) are given functions of the parameters j h , h , 2 , h with g "1, F (u, y) and E (u, y) are speci1 2 6 0 jk k fied functions of the input-output pair (u, y), and P (p) jk are fixed polynomials of degree n in the differential operator p"d/dt as they are defined in Table 1.
4. Parameter-estimation algorithm The HMF-method uses modulating functions to convert a set of differential equations that represent the real dynamic systems into an exact algebraic form. Now, multiplying both sides of Eq. (12) by / (t), m integrating over the observation time interval [0, ¹ ],
H
Refer to the appendix for the proof. Now, let z(mu )"H ? HM (2), then Eq. (14) can be rewritten as a 0 u y linear regression of the form z(mu )"uT(mu )h#e(mu ), (16) 0 0 0 where uT(mu )"[1 H ? H (2), !1 H ? HM (2), !2 H 3 u 3 y 3 u 0 y u ? HM (1), 2 H ? HM (1), !2 H ? HM u3, !2 H ? HM y, 2 q HM u2 , 3 y 3 u 3 l u y 3 y !2 u HM (1), 2 u HM ] and 3 a u 3 a u h"[1, 1, (h #h ), h , h , h , h h , h h , h , h h ]T. 1 4 1 2 5 1 4 2 6 3 3 4 Furthermore, the parameter vector h is to be estimated from z(mu ) and u(mu ) for m"0, $1, 2 ,$M, which 0 0 can be obtained from the evaluation of input-output data
620
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
through the HMF method. Minimizing a cost function J (h) given J(h)"1 eT(Mu , h)e(Mu , h) 2 0 0 "1 [z(Mu )!W(Mu )h]T[z(Mu )!W(Mu )h] 2 0 0 0 0 (17) with respect to the unknown parameter vector h leads to the least-squares (LS) estimate hK "[WT(Mu ) W(Mu )]~1 WT(Mu ) z(Mu ), 0 0 0 0
(18a)
where WT(Mu )"[u(!Mu )2u(!u )u(0) u(u )2 0 0 0 0 u(Mu )]. 0 zT(Mu )"[z(!Mu )2z(!u ) z(0)z(u )2z(Mu )] 0 0 0 0 0 and eT(Mu )"[e(!Mu )2e(!u )e(0)e(u )2e(Mu )]. 0 0 0 0 0 Moreover, introducing a positive definite symmetric weighting matrix W in the cost function J (h) leads to the weighted LS estimate of h given by hK "[WT(mu )W W(mu )]~1 WT(mu )Wz(mu ). W 0 0 0 0
ac-network and supplies it to the separately excited dc motor as shown in Fig. 2b. Here, two different cases for study have been considered: Case I — Window size [0, ¹ ]"[0, 2.56 sec.], sampling time t "40]10~3 sec. and HT order of st N"64; Case II — Window size [0, ¹ ]"[0, 3.2 sec.], t "50]10~3 sec. with HT order of N"64. A constant st input armature terminal voltage of u "60 [V], a a sinusoidal field current i (t) and a constant load torque f q "1.47 [N m] for case I, as well as q "1.65 [N m] for l l case II, are applied to the plant through the real-time operating-system software. Then, the observed input field current i (t), the output armature current i (t) and the f a motor speed n(t)"60 u(t)/2n [rpm] have been obtained as shown in Figs. 4, 5 and 6, respectively for both cases of study. Because of the noisy nature of the real-time system, a number of measurements with identical operating conditions have been considered for parameter-estimation purposes. Generally, having the input-output responses, the first main procedure is the computation of the HMF spectra of the differential equation describing the plant. Here, Hartley transforms are computed using the standard parabolic rule for better approximation. The second
(18b)
5. Experimental results This section presents experimental results using the approach outlined in the previous sections. In general, the experimental set-up consists of a thyristor-fed 1.1 kW separately excited dc-motor, an eddy-current brake, a magnetic clutch with a mass and a tacho-generator as shown in Fig. 2b. The experimental laboratory set-up is designed for variations in field excitation, moment of inertia and load. A variation of the load torque q can be l obtained by altering the excitation current of the eddycurrent brake i . Voltage, current, speed and torque b measurement instruments are available in the plant. Furthermore, the armature power controlled rectifier contains two thyristor bridges, each with four thristors in fully controlled single-phase configuration, and the field excitation supply unit is a single-phase half-controlled thyristor converter with current controller. The main functional assemblies of the converters are the power section, the trigger section, the control section and the rectifier for supplying power to the motor. The output of the converter can be changed by adjusting the timing of the thyristor firing pulses, and the necessary control pulses are supplied by the trigger set. Further, the output voltage of the current controller influences the timing of the control pulses, and thus also the firing angle. In addition, in the specified control range the converter will operate as a rectifier, i.e. it takes power from the
Fig. 4. Samples of observed input field currents, i (t). f
Fig. 5. Samples of output armature currents, i (t). a
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
621
Table 2 Estimated HMF-model parameter values of the plant Model hK
1
5
Number of runs 10 15
20
Case I: [0, ¹ ]"[0, 2.56sec], t "40]10~3 sec, q "1.47 Nm st l hK #hK 1 4 hK 1 hK hK 2 5 hK hK 1 4 hK hK 2 6 hK 3 hK hK 3 4 Fig. 6. Samples of observed motor speeds of the dc motor plant, n(t).
main procedure is the estimation of the HMF-model parameters using the weighted least squares, with a weighting factor inversely proportional to frequency. Finally, the two main procedures above are performed for a required number of runs, and then the average values of the parameters are calculated. Applying the aforementioned HMF method for case studies I and II gives the estimated HMF-model parameters and physical parameters shown in Tables 2 and 3, respectively for different numbers of runs. In the parameter-estimation study only (2M#1)"17 significant spectral data points were used. The physical parameters of the plant can be estimated using the following relations of the HMF-model parameters, i.e., armature inductance ¸K "1/hK [H]; armature resistance RK "hK /hK a 3 a 1 3 [)]; motor constant kK "(hK hK )/(hK hK ) [N m/A2]; total m 2 5 2 6 system moment of inertia JK "(hK hK )hK /(hK hK )2 2 5 3 2 6 [N m sec2], and viscous friction constant BK "(hK hK ) 1 4 (hK hK )hK /hK (hK hK )2 [N m sec]. 2 5 3 1 2 6 In general, it is not always possible to obtain the physical parameters directly/uniquely for a more complex system. It will depend on the mathematical relation between the physical system and its HMF-model parameters. Note that the main computations of the HMF method, given the input-output data, involve the approximation of the continuous-time signals’ Hartley transforms. These can be computed using a standard parabolic rule such as the extended Simpson’s rule for better approximation, given by
P
T
0
C
¹ x(t)cas(mu t)dt+ x(0) cas(0) 0 3N
A B A A B A D
#4
l¹ mu l¹ N~1 0 + x cas N N l"1, 3, 2
#2
l¹ mu l¹ N~2 0 + x cas N N l"2, 4, 2
#x(¹ ) cas(mu ¹ ) , 0
B B (19)
41.030 39.367 46.492 65.455 183.471 3.189 6.695
42.830 40.884 52.600 79.556 217.152 3.228 7.122
43.190 41.110 67.985 85.475 280.164 3.242 8.550
42.843 40.819 53.905 82.594 227.338 3.248 8.065
42.878 40.844 46.996 83.081 199.342 3.255 8.012
Case II: [0, ¹ ]"[0, 3.2 sec], t "50]10~3 sec, q "1.65 Nm st l hK #hK 1 4 hK 1 Kh hK 2 5 hK hK 1 4 Kh hK 2 6 hK 3 hK hK 3 4
38.849 37.346 41.882 56.124 176.008 3.045 6.018
43.685 41.981 54.028 71.529 200.881 3.309 7.092
43.127 41.385 50.959 72.084 208.535 3.241 7.438
43.463 41.798 47.676 69.618 200.103 3.243 7.238
43.575 41.807 50.834 73.908 208.016 3.223 7.602
JK [Nmsec2]
BK [Nmsec]
Table 3 Estimated physical parameters of the plant Number of runs
¸K a [H]
RK a [)]
kK m [Nm/A2]
Case I: [0, ¹ ]"[0, 2.56sec], t "40]10~3 sec, q "1.47 Nm st l 1 5 10 15 20
0.314 0.310 0.308 0.308 0.307
12.345 12.665 12.680 12.567 12.548
0.25340 0.24223 0.24266 0.23711 0.23576
0.00441 0.00360 0.00281 0.00339 0.00385
0.00732 0.00701 0.00584 0.00686 0.00783
Case II: [0, ¹ ]"[0, 3.2 sec], t "50]10~3 sec, q "1.65 Nm st l 1 5 10 15 20
0.328 0.302 0.308 0.308 0.310
12.265 12.687 12.769 12.889 12.971
0.23796 0.26896 0.24437 0.23826 0.24438
0.00412 0.00443 0.00380 0.00386 0.00379
0.00619 0.00755 0.00662 0.00643 0.00669
where ¹ is the window size and N is the number of samples within the fixed time interval (N"64). The other main computation of the methodology is the approximation of the HMF spectra of the input-output continuous signals, their derivatives, and the different form of products as given in Eqs. (5), (6), (7) and (15), respectively. Finally, comes the standard computation of the weighted least-squares algorithm to estimate the physical parameters of the nonlinear continuous-time dynamics. The main thing to be noted is that, in the computation of the parameter estimation, only M"8 significant spectral coefficients, or the mode number of
622
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626 Table 4 Estimated values of time constants Number of runs
1 5 10 15 20
Electrical time constant, tK []10~3 sec] e
Mechanical time constant, tK []10~3 sec] m
Case I
Case II
Case I
Case II
25.402 24.460 24.325 24.500 24.483
26.777 23.820 24.163 23.925 23.919
601.436 513.902 480.960 494.213 491.617
665.419 586.909 574.122 600.391 565.663
Fig. 7. HT and HMF spectrum of field current, i (t)"u(t). f
Fig. 8. HT and HMF spectrum of armature current, i (t)"y(t). a
Eq. (18b), have been used, even if the window size/record length has been discretized into N"64 subintervals as shown in Figs. 7 and 8 (the HT and HMF spectra versus frequency for the field current and armature current of case study I, respectively). Furthermore, based on the above estimated physical parameters, the electrical time constant, tK "¸K /RK " e a a 1/hK [sec] and the mechanical time constant, 1 tK "JK /BK "hK /(hK hK ) [sec] can be computed as shown in m 1 1 4 Table 4. The significant Hartley transform (HT) and the HMF spectra of the input-output signals for case I are shown in Figs. 7 and 8, respectively. The demonstration plant for variable-speed drives has been designed for general-purpose use, for instance, to understand the operating behavior of the dc machine (generator, motor) when separately excited, in shunt and compound connections, and with/without commutatingfield winding. This multipurpose machine has the following ratings: rated motor power 1.1 kW, voltage 130 V, current 11 A and speed 1500 rpm. According to the original data provided by the producer, the theoretical design values of the mechanical and electrical time constants of the plant, depending on the operating conditions, lie in the range of 3304t 4700 m []10~3 sec] and t +30 []10~3 sec], respectively. The e variation of the mechanical time constant can be achieved by varying the moment of inertia in the experi-
mental set-up, since the mechanical time constant will depend on the total system moment of inertia and the total system viscous friction constant. The load moment of inertia can be varied by changing the weight of the mass as shown in Fig. 2b. The experimental set-up has facilities to alter the mass, and as a result of this situation the mechanical time constant can vary correspondingly. Therefore, according to the operating conditions performed in the experiment, the values given in Table 4 are obtained, depending on different numbers of runs, and they lie within the range according to the mass taken in the experimental set-up. In general, if the minimum mass is coupled, the minimum mechanical time constant can be reached. On the other hand, if all the masses of the set-up are coupled, the maximum mechanical time constant can be reached. In addition, even if it is not possible to relate all values of the estimated physical parameters exactly to the design data of the general-purpose dcmachine set-up, one can clearly observe that the estimated values of the physical parameters are almost kept constant as the number of runs increases from 1 to 5, 10, 15 and 20, as shown in Tables 3 and 4 for both case studies. Moreover, the approximate design value of the armature inductance can be calculated, based on a method described in the literature (Langhoff and Raatz, 1977), and on dc-machine ratings. According to the experimental conditions involved here, an approximate design value of inductance ¸ "3¸ "3(º k)/(I N D )" a .*/ ar ar r a 0.26211 H can be taken into consideration, where º , ar I and N are the rated values of the dc-motor voltage, ar r current and speed, respectively; k is the compensation constant of the machine, k"1000; D "3JP /N and a r r P is the rated power of the motor. Furthermore, it may r be necessary to take into account that the estimated inductance ¸K +¸K #¸K , where ¸K is the motor a amot int amot armature inductance and ¸K is the inductance of the int operating effect or actuating mechanism. On the other hand, the estimated resistance RK +RK #RK #RK , a amot cab int where RK is the motor armature resistance including amot the commutating poles and compensating windings, RK is the resistance of the cables, and RK is the internal cab int resistance of the operating effect or actuating mechanism.
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Essentially, the experimental set-up is designed for multipurpose use, and was devised as a teaching, training aid and research set-up in the field of control engineering. As a result, some improvements and modifications have been implemented over time. Basically, the electrical time constant will depend on the armature resistance and inductance. In general, the deviation of the estimated electrical time constant from its approximate theoretical value may be due to the replacement of old cables, or replacement of actuating mechanisms, motor windings, compensating windings, etc. In general, the results show the performance of the HMF method in estimating the physical parameters of the nonlinear continuous-time dynamics of a thyristordriven dc-motor.
6. Conclusion In this paper, an application of the HMF-method to estimate the parameters of an experimental thyristordriven dc-motor laboratory set-up is carried out. The proposed method has given satisfactory results, and is promising for use in parameter estimation of physicallybased continuous-time nonlinear systems from their input-output data records, observed over a fixed time interval or window size, without the necessity for estimating unknown initial conditions. This recently developed methodology uses a frequency-dependent weighted leastsquares algorithm to identify the physical parameters. The practical implementation of the HMF method is facilitated by the CADACS computer-aided control-system design software. The armature terminal voltage and field motor circuit inputs are set from the PC of the real-time operating-system software of CADACS, and the input-output responses are displayed on the analog instruments, as well as on the PC screen via A/D converters. Then, the measured input-output responses are applied to an HMF algorithm of the dc-motor model to estimate the physical system parameters. The practical application of the HMF method and good results obtained here demonstrate the efficiency of this new approach. The proposed methodology is of considerable practical importance in problems of continuous-time parameter estimation in which input-output measurements are available on a finite interval of time. In general, the main advantage of the Hartley modulating functions method is that it is suitable for use in identifying the physically-based parameters of nonlinear continuoustime systems. Furthermore, the basic benefits of the HMF approach are that it is able to handle noisy data, and it does not need fixed/known initial conditions. In addition, the method employs integration-by-parts and the two-point boundary conditions of the modulating properties to transfer derivatives of data variables into smooth
623
moduating functions, in order to obtain an algebraic equation that is suitable for identification purposes. Generally, these features may make the modulating functions method desirable for use in several real processes. In the light of this encouraging result of the HMF method, further studies will continue on the computational considerations associated with the approach, and on how to implement this technique in the form of batch-scheme recursive estimation for nonlinear continuous-time systems. In principle, the HMF-method can be applied for real-time identification using a batch scheme recursive parameter-estimation technique. This can be implemented by shifting a fixed window size, or recorded batch data, one step forward at each sampling time. The main procedure of the batch scheme identification method is as follows: First, obtain an estimate based on initial batch of data; for instance one can assume the estimated parameters obtained by a single run shown in Tables 2, 3 or 4, then consider the new input-output sample data, and shift the fixed window size one step forward at each sampling time by appending a new column and deleting the first column of the data vector. Finally, the Hartley spectrums can be computed recursively, and the new estimate of the parameters obtained. This study of batchscheme recursive identification for nonlinear continuoustime systems is currently under investigation. Acknowledgements The first author thanks the German academic exchange service for providing the financial support for his stay at the Ruhr-University Bochum. He is also grateful to the Department of Electrical Engineering, Faculty of Technology, Addis Abada University, Addis Ababa, Ethiopia, for sanctioning leave of absence. References Langhoff, J., Raatz, E., 1977. Geregelte Gleichstromantriebe AEGTelefunken-Handbu¨cher. Band 19, Elitera-Verlag, Berlin. Hartley, R.V.L., 1942. A more symmetrical Fourier analysis applied to transmission problems. Proc. Inst. Radio Engrs., 144 —150. Jalili, S.A., Jordan, J.A., Mackie, R.D.L., 1992. Measurement of the parameters of all-pole transfer functions using shifted Hermite modulating functions. Automatica 38, 613— 617. Karanam, V.R., Frick, P.A., Mohler, R.R., 1978. Bilinear system identification by Walsh functions. IEEE Trans. Autom. Cont., AC-23, 709—713. Kung, F.C., Shih, D.H., 1986. Analysis and identification of Hammerstein model nonlinear delay systems using block-pulse function expansion. Int. J. Cont. 43(1), 139—147. Maletinsky, V., 1979. Identification of continuous dynamical systems with spline-type modulating functions method. 5th IFAC Symp. on Ident. and Syst. Par. Est. Germany, 275—281. Patra, A., 1989. General hybrid orthogonal functions and some applications in systems and control. Ph.D. Thesis, Department of Electrical Engineering, Indian Institute of Technology, Kharagpur-721302, India.
624
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Patra, A., Unbehauen, H., 1995. Identification of a class of nonlinear continuous-time systems using Hartley modulating functions. Int. Journal of Control 62(2), 1431—1451. Pearson, A.E., 1992. Explicit parameter identification for a class of nonlinear I/O-differential operator models. Proc. of 31st IEEE CDC, Arizona, 3656—3660. Pearson, A.E., Lee, F.C., 1983. Time limited identification of continuous systems using trigonometric modulating functions. Proc. of 3rd Yale Workshop on Applic. of Adap. Syst., New Haven, CT, pp. 168—173. Pearson, A.E., Lee, F.C., 1985. On the identification of polynomial I-O-differential systems. IEEE Trans. Autom. Cont. AC-30(8), 778—782. Saha, D.C., Roy, B.K., Bapat, V.N., 1991. Some generalizations and further uses of the Poisson moment functionals (PMF ) approach to continuous-time model estimation. Proc. of the 9th IFAC Symp. on Ident. and Par. Estim., Budapest, Hungary, pp. 1324—1333. Schmid, C., 1992. Real-time control with CADACS-PC. In: Jamshidi, M., Herget, C.J. (Eds.) Recent advances in computer-aided control System Engineering, pp. 337—355, Elsevier Science Publishers, B.V., Amsterdam. Sinha, N.K., Rao, G.P., (eds.) 1991. Identification of continuous-time systems methodology and computer implementation. Kluwer Academic Publishers, Dordrecht, The Netherlands. Shinbrot, M., 1957, On the analysis of linear and nonlinear systems. Trans. ASME 79, 547—552. Subrahmanyam, A.V.B., Saha, D.C., Rao, G.P., 1995. Identification of continuous-time MIMO systems via time moments. ControlTheory and Advanced Technology 10(4), 1359 —1378. Tsang, K.M., Billings, S.A., 1994. Identification of continuous-time nonlinear systems using delayed state variable filters. Int. J. Cont. 60(2), 159—180. Unbehauen, H., 1993, CADACS: An integrated design facility for computer-aided design and analysis of control systems. In: Linkens, D.A. (Ed.), CAD for control systems, pp. 469— 487, Marcel Dekken Inc., New York. Unbehauen, H., Rao, G.P., 1987. Identification of continuous systems. North-Holland, Amsterdam. Unbehauen, H., Rao, G.P., 1990. Continuous-time approaches to system identification — a survey. Automatica 26(1), 23—35. Wang, M.L., Chen, K.S., Chou, C.K., 1984. Solution of integral equations via modified Laguerre polynomials. Int. J. Syst. Sci. 15(6), 661—672.
AB
Wang, C.H., Marleau, R.S., 1987. Recursive computational algorithms for the generalized block pulse operational matrices. Int. J. Cont. 45(1), 195—201. Young, P.C., 1981. Parameter estimation of continuous time models — a survey. Automatica 17(1), 23—39.
Appendix – Proofs [A]. Proof of Eq. (15). Let H ? H (i) be the m-th HMF spectral component of y u the i-th derivative of the product of two signals Mu(t)(i) y(t)N given by
P
T H ? H (i)" Mu(t)(i) y(t)N/ (t) dt. (20) u m y 0 Applying integration-by-parts i times while substituting the two end-point boundary conditions of / (t) given by m /(i) (0)"/(i) (¹ )"0, one can transfer the derivative of m m the signals to the derivative of the known modulating function, i.e.
P
AB
]cas((!1)i(n#m!j )u t) dt. 0
(21)
Now, replacing one of the signals by its infinite inverse Hartley series, say for u(t), = u(t)" + H (lu ) cas(lu t) u 0 0 l/~= results in
A BP
= n n in H ? H (i)" + H (lu ) + (!1)j (n#m!j )i ui cas@ y u 0 u 0 j 2 l/~= j/0
AB
P
T H ? H (i)"(!1)i u(t) y(t)M/(i) (t)N dt u y m 0 T n n " u(t) y(t)(!1)i + (!1)j j 0 j/0 ](n#m!j)iui (!1)i cas@(in/2) 0
A B
T 0
= n n in 1 " + H (lu ) + (!1)j (n#m!j )i ui cas@ ] u 0 0 j 2 2 l/~= j/0
(22)
y(t) cas(lu t) cas ((!1)i(n#m!j )u t) dt 0 0
G
P
T
y(t)cas ((l#(!1)i(n#m!j )) u t) dt 0
0
P P P
!
T
#
0 T
#
0 T 0
y(t)cas ((!l!(!1)i(n#m!j )) u t) dt 0 y(t)cas ((l!(!1)i(n#m!j )) u t) dt 0 y(t)cas ((!l#(!1)i(n#m!j )) u t) dt 0
since cas(m) cas(f)"1 [cas(m#f)!cas(!m!f)#cas(m!f)#cas(!m#f)]. 2
H
,
(23)
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
625
Replacing the integral terms of (23) by their corresponding Hartley transforms leads to the preceding equation of the form
AB
A B
= n n in H ? H (i)" + H (lu ) + (!1)j (n#m!j )i ui cas@ u y u 0 0 j 2 l/~= j/0
G
H
1 H ((l#(!1)i(n#m!j ))u )!H ((!l!(!1)i(n#m!j ))u ) y 0 y 0 ] . 2 #H ((l!(!1)i (n#m!j ))u )#H ((!l#(!1)i (n#m!j ))u ) y 0 y 0
(24)
For the case of Eq. (15), H ? H (2), substituting 2 in place of i gives the expression described in Eq. (15). This completes u y the proof for the general case. [B]. Proof of Eq. (7). Let the signals x (t) and x (t) be in the form of x (t)x(i) (t); then the HMF spectrum for such a product is given by 1 2 1 2
P
T
P
T
HM (0, x , xi) (mu )" 1 2 0
"
0
x (t) x(i) (t)/ (t) dt 1 2 m
AB
n n x (t) x(i) (t) + (!1)j 1 2 j 0 j/0
]cas((n#m!j )u t) dt. 0
(25)
Replacing x (t) by its infinite inverse Hartley series, i.e., x (t)"+ = H (lu ) cas (lu t), implies 1 1 0 l/~= x1 0
P
AB
T n n = HM x(0,, xi) (mu )" + Hx (lu ) x(i) (t) cas (lu t) + (!1)j cas((n#m!j ) u t) dt. 1 2 1 2 0 0 0 0 j 0 j/0 l/~=
(26)
Now, utilizing the identity given by cas(m) cas(f) "1 [cas(m#f)!cas(!m!f)#cas(m!f)#cas(!m#f)] into (26) 2 leads to
A BP
n = n HM (0, x , xi) (mu )" + Hx (lu ) + (!1)j 1 2 1 0 0 j l/~= j/0
T 0
x(i) (t) 1 Mcas((l#m#n!j ) u t)!cas ((!l!m!n#j) u t) 2 2 0 0
#cas((l!m!n!j )u t)#cas ((!l#m#n!j ) u t)N dt. 0 0
(27)
Applying integration-by-parts i times while substituting the two end-point boundary conditions of / (t) given by m /(i) (0)"/(i) (¹ )"0, one can transfer the derivative of the signals to the derivative of the known modulating m m function, i.e.
G
P
T
(!1)ix (t) cas(i)((l#m#n!j )u t) dt 2 0
0
P P P
T
(!1)ix (t) cas(i)((!l!m!n#j )u t) dt 2 0 n = n 1 0 HM (0, ] x , xi) (mu )" + Hx (lu ) + (!1)j 1 2 1 0 0 j 2 T l/~= j/0 # (!1)ix (t) cas(i)((l!m!n#j )u t) dt 2 0 0 T # (!1)ix (t) cas(i)((!l#m#n!j )u t) dt 2 0 0
AB
!
H
626
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
AB
A B
n in = n " + Hx (lu ) + (!1)j ui cas@ 1 0 0 j 2 l/~= j/0
1 ] 2
G
(l#m#n!j)i
P
T
0
x (t) cas((!1)i (l#m#n!j )u t) dt 2 0
P
!(!l!m!n#j)i
#(l!m!n#j )i
P
T
0
#(!l#m#n!j)i
T 0
x (t) cas((!1)i (!l!m!n#j )u t) dt 2 0
x (t) cas((!1)i (l!m!n#j )u t) dt 2 0
P
T
0
x (t) cas ((!1)i (!l#m#n!j )u t) dt 2 0
H
(28)
since cas(i)(qu t)"(!1)icas@(in/2)qiui cas((!1)iqu t). Replacement of the integral terms of Eq. (28) by their corre0 0 0 sponding Hartley transform symbols leads to the preceding equation of the form
G
(l#m#n!j )i Hx ((!1)i(l#m#n!j )u ) 2 0 n in = n 1 !(!l!m!n#j )iHx2 ((!1)i(!l!m!n#j )u0 ) HM x(0,, xi) (mu )" + Hx (lu ) + (!1)j ui cas@ ] 1 2 1 0 0 0 j 2 2 #(l!m!n#j )i Hx ((!1)i(l!m!n#j )u ) 2 l/~= j/0 0 #(!l#m#n!j )iHx ((!1)i(!l#m#n!j )u ) 2 0
AB
"Hx (mu ) ? HM x(i) (mu ). 1 2 0 0 This completes the proof of (7).
A B
H
(29)
Experimental physical parameter estimation of a thyristor driven DC-motor using the HMF-method1 S. Daniel-Berhe, H. Unbehauen* Control Engineering Laboratory, Faculty of Electrical Engineering, Ruhr-University Bochum, ESR IC 3/131, D-44780 Bochum, Germany Received 28 April 1997; accepted 17 January 1998
Abstract This paper deals with the physical parameter estimation of the nonlinear continuous-time dynamics of a thyristor-driven dc-motor experimental set-up using the Hartley modulating functions (HMF) method. The approach applies a frequency-weighted leastsquares formulation to estimate the nonlinear continuous-time system parameters based on input-output data records over a finite time interval. The laboratory experimental set-up of the plant is facilitated by a computer-aided control system design software system, CADACS. The method has given satisfactory experimental results, and demonstrates its feasibility for use in parameter estimation of physically-based continuous-time systems. Due to its encouraging results, further studies will continue, in order to implement the methodology in the form of batch scheme recursive estimation for nonlinear continuous-time dynamics. ( 1998 Elsevier Science ¸td. All rights reserved. Keywords: Physical parameter estimation; dc motors; modulating functions; nonlinear; continuous-time systems
1. Introduction Most real, dynamic systems are continuous in time and nonlinear, and can be represented by a set of differential equations in which the parameters are physically defined. For numerical treatment, the dynamics of real systems are usually described by discrete-time models, for instance in the form of difference equations. However, a number of practical and theoretical considerations indicate that continuous-time treatment of these systems is often relevant, leading to continuous-time parameterestimation methods. Furthermore, identification of nonlinear continuous-time systems from available sampled input-output data records is desirable, to enhance the interpretation of the nonlinear dynamic system behavior, and to support the controller design.
* All correspondence should be sent to Professor Dr.-Ing. Heinz Unbehauen. Tel.: ##49 234 700 4071; fax ##49 234 709-4101; e-mail: [email protected] 1 The original version of this paper was presented at the 11th IFAC Symposium on System Identification, which was held in Fukuoka, Japan during 8—11 July 1997.
In recent years the identification of continuous-time linear and nonlinear systems is a topic that has received considerable attention. For instance, the relevance and importance of continuous-time system identification, with detailed reviews, have been reported in the literature by Young (1981), Unbehauen and Rao (1987, 1990). Sinha and Rao (1991), etc. These works discussed the relationship of the estimated parameters to their physically-based continuous-time system descriptions, and reviewed the progress of research on parameter estimation for continuous-time models of dynamic systems over the past few decades. Furthermore, the need for models/methods for nonlinear continuous-time systems identification is due to the fact that a conventional continuous-time linear model holds in only a small neighborhood of the actual equilibrium state, and leads to limitations in applications. Today, quite a number of methods are available for the parameter identification of linear continuous-time systems. Unfortunately these methods cannot be applied for the identification of nonlinear continuous-time systems. The use of standard linear techniques has been mostly restricted to only very special classes of nonlinear systems, such as bilinear systems, Hammerstein systems, etc.
0967-0661/98/$ — See front matter ( 1998 Elsevier Science Ltd. All rights reserved PII S 0 9 6 7 - 0 6 6 1 ( 9 8 ) 0 0 0 3 6 - 7
616
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
In general, some of the familiar approaches of continuous-time systems identification include Poisson moment functionals (Saha et al., 1991), block pulse functions (Kung and Shih, 1986), Walsh functions (Karanam et al., 1978), orthogonal polynomials (Wang and Marteau, 1987), general hybrid orthogonal functions (Patra, 1989), Laguerre polynomials (Wang et al., 1984), time moments (Subrahmanyam et al., 1995), delayed state variable filters (Tsang and Billings, 1994), Shinbrot’s method of moment functional (Shinbrot, 1957), Hermite modulating functions (Jalili et al., 1992), trigonometric modulating functions (Pearson and Lee, 1983), Fourier modulating functions (Pearson and Lee, 1985; Pearson, 1992), splinetype modulating functions (Maletinsky, 1979), Hartley modulating functions methods (Patra and Unbehauen, 1995). Most of the aforementioned methods were first developed for linear continuous-time systems, and further developed for special nonlinear systems, and have shown promising results. Generally, they require further studies for nonlinear continuous-time system identification. Relatively, the literature of moduling functions methods shows very promising results for the identification of a large class of nonlinear continuous-time systems (Unbehauen and Rao, 1987). The concept of modulating functions, introduced by Shinbrot in 1957, is commonly known as Shinbrot’s method of moment functionals. Essentially, the use of modulating functions involves converting a differential expression involving input-output signals on a specified time interval, into a sequence of algebraic equations. Shinbrot explored parameter identification of higher-order nonlinear dynamic systems by applying integral transformation techniques. In particular, Fourier modulating functions methods have shown encouraging performance for the identification of linear, integrable as well as convolvable, nonlinear continuoustime models (Pearson and Lee 1985, Pearson 1992). Basically, the Hartley modulating functions (Patra and Unbehauen, 1995) is closely related to the Fourier modulating functions. Compared to the very efficient Fourier modulating functions method, the HMF method has additional important advantages: the Hartley modulating functions are real-valued, and their Hartley spectra can be computed efficiently with the help of fast algorithms. This new methodology is applicable to a large class of nonlinear continuous-time systems by defining a set of HMF for characterizing the continuous process signals. The approach has adopted a class of nonlinear differential operator models, essentially developed for Fourier modulating functions (Pearson, 1992), which has been applied to parameter identification by formulating an efficient frequency-weighted least-squares algorithm. It seems that modulating functions methods provide the most promising alternative for use in parameter estimation of physically-based nonlinear continuous-time
systems (Unbehauen and Rao, 1987). Therefore, in this paper a new, feasible continuous-time parameter estimation method, based on the Hartley modulating functions approach is applied, which directly relates the estimated parameters to their physically-based continuous-time systems representation. Basically, the demonstration of this experimental set-up is intended to show the feasibility of the modulating functions method for physical parameter estimation of nonlinear continuoustime systems. This paper is organized as follows. In Section 2, the background of the Hartley modulating functions, and the Hartley transform and spectra are discussed. The experimental set-up and plant description are outlined in Section 3. A new parameter-estimation algorithm is discussed in Section 4. Experimental results are documented in Section 5. Finally, in Section 6, conclusions and future research plans are presented.
2. Background The general idea of modulating functions was first introduced by Shinbrot (1957). The application of modulating functions to convert a differential equation into an exact algebraic form is well known in the field of system identification, depending on several choices of modulating functions, for instance Shinbrot’s method of moment functionals (Shinbrot, 1957), spline-type modulating functions (Maletinsky, 1979), Hermite modulating functions (Jalili et al., 1992), trigonometric modulating functions (Pearson and Lee, 1983), Fourier modulating functions (Pearson and Lee, 1985; Pearson, 1992), Hartley modulating functions methods (Patra and Unbehauen, 1995), etc. The Hartley modulating functions have recently been introduced as a new member of the modulating functions family. This modulating function is closely related to the Fourier modulating functions. A member of the family of an nth-order HMF / (t), m"0, $1, $2,2 relative to m a fixed time interval [0, ¹ ] is given by
AB
n n / (t)" + (!1)i cas(n#m!i)u t, m 0 i i/0
(1)
where m is the modulating frequency index, cas ut"cos ut#sin ut and u "2n/¹ is essentially the 0 resolving frequency. Furthermore, / (t) is an HMF of m order n on a fixed interval [0, ¹] if it is sufficiently smoother and satisfies for each m"0, $1, $2,2, $M the two-point boundary conditions dk / (t)"0, dtk m
for t"0 and t"¹, k"0, 1,2, n!1
as depicted in Fig. 1.
(2)
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
617
The approach uses Hartley transforms and spectra. The continuous and discrete Hartley transforms (HT) (Hartley, 1942) of a signal x are given by =
P~= x(t)cas ut dt
H (u)"HMx(t)N" x
(3)
and
G A BH
A B A B
k¹ 2nlk 1 N~1 " + x cas (4) N N N k/0 respectively, for l"0, 1,2, N!1, where N is the number of samples within the fixed interval [0, ¹]. The mth HMF spectral component of a signal x(t) is given by k¹ HK (l)"H x x N
T
P0
HM (mu )" x 0
Fig. 1. The family of 2nd-order HMF, / (t), m"0, 1, 2; n"2 and m t3[0, ¹ ].
x(t)/ (t) dt m
AB AB
n :T x(t)cas(n#m!j)u t dt n 0 0 " + (!1)j j hggggiggggj j/0 Hx ((n`m~j)u0) n n H ((n#m!j)u ). (5) " + (!1)j x 0 j j/0 If x(i)(t), i"1, 2,2, n is the ith derivative of the signal x(t), then its corresponding HMF spectrum is given by
ristor-driven dc-motor system and analog input-output display instruments, as shown in Fig. 2a. The armature terminal and field motor circuit input voltages are set from a PC, using computer-aided design and analysis of control systems (CADACS) software (Unbehauen, 1993; Schmid, 1992) through D/A converters. The software has a capability to link the results of a computer-aided analysis and design procedure with the real-time experimental environment. Consequently, the input-output responses can be displayed on the PC screen via converters and interfacing circuits, as well as on the analog instruments. Fig. 2c shows the standard model of a dc motor. The main aim is to apply the measured input-output data to a suitable HMF algorithm of the dc-motor model, and estimate the physically-based system parameters.
AB
n n HM (i) (mu )" + (!1)j cas@(in/2)(n#m!j )i x 0 j j/0 ]ui H ((!1)i(n#m!j)u ), (6) 0 x 0 where cas@(t)"cos(t)!sin(t). Moreover, if the signals x (t) and x (t) are in the form of x (t)x(i) (t), then the 1 2 1 2 HMF spectrum for such a product is given by
G
(l#m#n!j)iH ((!1)i(l#m#n!j)u ) x2 0 n in = n 1 !(!l!m!n#j)iHx2 ((!1)i(!l!m!n#j)u0) HM (0,i) (mu )" + H (lu ) + (!1)j ui cas@ ] x1,x2 0 x1 0 0 j 2 2 #(l!m!n#j)iH ((!1)i(l!m!n#j)u ) l/~= j/0 x2 0 #(!l#m#n!j)iH ((!1)i(!l#m#n!j)u ) x2 0 "H (mu ) ? HM (i)(mu ), x1 0 x2 0
AB
AB
where the operator ? symbol represents the description in short form for simplicity. The proof of (7) is given in the appendix.
3. Experimental set-up and plant description The configuration of the experimental thyristor-driven dc-motor system is shown in Fig. 2(a—c), along with a detailed circuit of the dc motor and load. Primarily, the experimental set-up of the plant consists of a personal computer (PC) with real-time operating system software and interfacing circuits, digital-to-analog (D/A) and analog-to-digital (A/D) converters, a thy-
H
(7)
In general, dc-motor systems have extensive applications, varying from industrial drivers to medical and other instrumentation controls, because of their versatile speed-load characteristics. A dc-motor system can be represented mathematically by the following differential equations: di (t) f u (t)"R i (t)#¸ f ff f dt di (t) a #k i (t)u(t) u (t)"R i (t)#¸ a aa a dt mf J
du(t) "k i (t)i (t)!q (t)!Bu(t) mf a l dt
618
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Fig. 3. Block diagram of the dc motor and load experimental set-up.
Fig. 2(a—c). Laboratory experimental thyristor-driven dc motor system layout. (a) Scheme of the experimental set-up under study. (b) Part of the thyristor-driven dc motor installation. (c) Equivalent circuit representation of the dc motor and load.
q (t)"k t(t)i (t)"k k i (t)i (t)"k i (t)i (t) m 1 a 1 2f a mf a u (t)"k t(t)u(t)"k k i (t)u(t)"k i (t)u(t) (8) i 1 1 2f mf where u (t) is the input voltage for field motor circuit [V], & i (t), the field motor current [A], R , the resistance of the & & field motor circuit [)], ¸ , inductance of the field motor & circuit [H], u (t), armature terminal voltage source [V], ! i (t), armature current [A], R , armature resistance [)], ! ! ¸ , armature inductance [H], u (t), induced voltage in the ! * armature coil [V], u(t), motor angular speed [rad/sec], q (t), electromagnetic torque of the motor [Nm], q (t), . l load torque [Nm], J , J and J, the load, motor and total l . system moment of inertia [N m sec2], B , B and B, the l . load, motor and total system viscous friction constant [N m sec], t(t)"k i (t), the magnetic flux [N m/A], and 2& k "k k , is known as a motor constant [N m/A2]. . 1 2 In general, a dc motor has three inputs, namely, the voltage applied to the armature coil u (t), the voltage ! applied to the field coil u (t), and the load torque q (t) that & l results from the (assumed) machinery that is being driven by the motor. Accordingly, the system of equations in Eq. (8) represents a model of the dc-motor system as a whole, since it consists of the two electrical subsystems describing the field and armature, the mechanical subsystem describing the inertia and friction of the rotating mass, plus the two coupling equations that connect the mechanical subsystem to the electrical subsystem. Moreover, the dc-motor system can be modelled by the simulation diagram shown in Fig. 3. The problem is to estimate the physical parameters of the plant, knowing its structure and based on the inputoutput data records over a finite time interval. In the experimental study a constant load torque and armature terminal voltage sources are considered. Thus, letting the input field current i (t)"u(t), the armature & current i (t)"x (t) and motor angular speed u(t)" ! 1 x (t), then the system equations become 2 R k u xR (t)"! a x (t)! m u(t)x (t)# a 1 2 ¸ 1 ¸ ¸ a a a B k q xR (t)"! x (t)# m u(t)x (t)! l (9) 2 1 J 2 J J y(t)"x (t) 1
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
or
G
yR (t)"!h y(t)!h u(t)x (t)#h u , 1 2 2 3 a (10) xR (t)"!h x (t)#h u(t) y(t)!h q 2 4 2 5 6 l where h "R /¸ , h "k /¸ , h "1/¸ , h "B/J, 1 a a 2 m a 3 a 4 h "k /J and h "1/J. 5 m 6 Eliminating the motor angular speed x (t) and sim2 plifying, the preceding equation in the input field current u(t) and output armature current y(t) yields
u(t) y¨ (t)!uR (t) yR (t)#(h #h ) u(t) yR (t)!h y(t) uR (t) 1 4 1 #h h u3(t)y(t)#h h u(t) y(t)!q h h u2(t) 2 5 1 4 1 2 6 #u h uR (t)!u h h u(t)"0. (11) a 3 a 3 4 Applying the identity uR (t) yR (t)"1 [d2/dt2Mu(t) y(t)N! 2 y(t) u( (t)!u(t) y( (t)] in Eq. (11) and normalizing its first element yields to the standard form of a nonlinear inputoutput differential operator model (Pearson, 1992) of the form u(t) p2y(t)!1 p2Mu(t) y(t)N#1 y(t) p2u(t) 3 3 #2 (h #h ) u(t) py(t)!2 h y(t) pu(t) 3 1 4 3 1 #2 h h u3(t) y(t)#2 h h u(t) y(t)!2 q h h u2(t) 3 2 5 3 1 4 3 l 2 6 #2 u h pu(t)!2 u h h u(t)"0 (12) 3 a 3 3 a 3 4 or m
1
m
2
+ + g (h) F (u, y)P (p) E (u, y)"0 j jk jk k j/0 k/1
(13)
619
Table 1 The functions g (h), F (u, y), P (p) and E (u, y) for the nonlinear model j jk jk k of Eq. (12) j
k
g (h) j
F (u, y) jk
P (p) jk
E (u, y) k
0 1
1 1 2 1 1 1 1 1 1 1
1 1
u(t) !1 3 1 y(t) 3 2 u(t) 3 !2 y(t) 3 2 u3(t) y(t) 3 2 u(t) y(t) 3 !2 q u2(t) 3 l 2u 3 a !2 u u(t) 3 a
p2 p2 p2 p1 p1 p0 p0 p0 p1 p0
y(t) u(t) y(t) u(t) y(t) u(t) 1 1 1 u(t) 1
2 3 4 5 6 7 8
h #h 1 4 h 1 h h 2 5 h h 1 4 h h 2 6 h 3 h h 3 4
and applying their corresponding HMF spectrum notations discussed above, leads to the HMF-model of the system H ? HM (2)"1 H ? H (2)!1 H ? HM (2) 3 u u y 3 y y u !2 (h #h )H ? HM (1)#2 h H ? HM (1) 3 1 y 3 1 4 u y u !2 h h H ? HM 3!2 h h H ? HM u 3 1 4 u y 3 2 5 y #2 q h h HM 2!2 u h HM (1)#2 u h h HM , 3 l 2 6 u 3 a 3 u 3 a 3 4 u (14) q where H ? H (2) represents : Mu(t)® y(t)N/ (t) dt, and apu m y 0 plying integration-by-parts replacing one of the signal by its infinite Hartley series and using the identity for the product of cas function implies
H ? H (2)" u y n 2n H ((l#m#n!j)u )!H ((!l!m!n#j)u ) = n y 0 y 0 (n#m!j)2u2 cas@ ]1 . + H (lu ) + (!1)j 2 #H ((l!m!n#j)u )#H ((!l#m#n!j)u ) 0 u 0 j 2 y 0 y 0 l/~= j/0 (15)
AB
A B G
where g (h) are given functions of the parameters j h , h , 2 , h with g "1, F (u, y) and E (u, y) are speci1 2 6 0 jk k fied functions of the input-output pair (u, y), and P (p) jk are fixed polynomials of degree n in the differential operator p"d/dt as they are defined in Table 1.
4. Parameter-estimation algorithm The HMF-method uses modulating functions to convert a set of differential equations that represent the real dynamic systems into an exact algebraic form. Now, multiplying both sides of Eq. (12) by / (t), m integrating over the observation time interval [0, ¹ ],
H
Refer to the appendix for the proof. Now, let z(mu )"H ? HM (2), then Eq. (14) can be rewritten as a 0 u y linear regression of the form z(mu )"uT(mu )h#e(mu ), (16) 0 0 0 where uT(mu )"[1 H ? H (2), !1 H ? HM (2), !2 H 3 u 3 y 3 u 0 y u ? HM (1), 2 H ? HM (1), !2 H ? HM u3, !2 H ? HM y, 2 q HM u2 , 3 y 3 u 3 l u y 3 y !2 u HM (1), 2 u HM ] and 3 a u 3 a u h"[1, 1, (h #h ), h , h , h , h h , h h , h , h h ]T. 1 4 1 2 5 1 4 2 6 3 3 4 Furthermore, the parameter vector h is to be estimated from z(mu ) and u(mu ) for m"0, $1, 2 ,$M, which 0 0 can be obtained from the evaluation of input-output data
620
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
through the HMF method. Minimizing a cost function J (h) given J(h)"1 eT(Mu , h)e(Mu , h) 2 0 0 "1 [z(Mu )!W(Mu )h]T[z(Mu )!W(Mu )h] 2 0 0 0 0 (17) with respect to the unknown parameter vector h leads to the least-squares (LS) estimate hK "[WT(Mu ) W(Mu )]~1 WT(Mu ) z(Mu ), 0 0 0 0
(18a)
where WT(Mu )"[u(!Mu )2u(!u )u(0) u(u )2 0 0 0 0 u(Mu )]. 0 zT(Mu )"[z(!Mu )2z(!u ) z(0)z(u )2z(Mu )] 0 0 0 0 0 and eT(Mu )"[e(!Mu )2e(!u )e(0)e(u )2e(Mu )]. 0 0 0 0 0 Moreover, introducing a positive definite symmetric weighting matrix W in the cost function J (h) leads to the weighted LS estimate of h given by hK "[WT(mu )W W(mu )]~1 WT(mu )Wz(mu ). W 0 0 0 0
ac-network and supplies it to the separately excited dc motor as shown in Fig. 2b. Here, two different cases for study have been considered: Case I — Window size [0, ¹ ]"[0, 2.56 sec.], sampling time t "40]10~3 sec. and HT order of st N"64; Case II — Window size [0, ¹ ]"[0, 3.2 sec.], t "50]10~3 sec. with HT order of N"64. A constant st input armature terminal voltage of u "60 [V], a a sinusoidal field current i (t) and a constant load torque f q "1.47 [N m] for case I, as well as q "1.65 [N m] for l l case II, are applied to the plant through the real-time operating-system software. Then, the observed input field current i (t), the output armature current i (t) and the f a motor speed n(t)"60 u(t)/2n [rpm] have been obtained as shown in Figs. 4, 5 and 6, respectively for both cases of study. Because of the noisy nature of the real-time system, a number of measurements with identical operating conditions have been considered for parameter-estimation purposes. Generally, having the input-output responses, the first main procedure is the computation of the HMF spectra of the differential equation describing the plant. Here, Hartley transforms are computed using the standard parabolic rule for better approximation. The second
(18b)
5. Experimental results This section presents experimental results using the approach outlined in the previous sections. In general, the experimental set-up consists of a thyristor-fed 1.1 kW separately excited dc-motor, an eddy-current brake, a magnetic clutch with a mass and a tacho-generator as shown in Fig. 2b. The experimental laboratory set-up is designed for variations in field excitation, moment of inertia and load. A variation of the load torque q can be l obtained by altering the excitation current of the eddycurrent brake i . Voltage, current, speed and torque b measurement instruments are available in the plant. Furthermore, the armature power controlled rectifier contains two thyristor bridges, each with four thristors in fully controlled single-phase configuration, and the field excitation supply unit is a single-phase half-controlled thyristor converter with current controller. The main functional assemblies of the converters are the power section, the trigger section, the control section and the rectifier for supplying power to the motor. The output of the converter can be changed by adjusting the timing of the thyristor firing pulses, and the necessary control pulses are supplied by the trigger set. Further, the output voltage of the current controller influences the timing of the control pulses, and thus also the firing angle. In addition, in the specified control range the converter will operate as a rectifier, i.e. it takes power from the
Fig. 4. Samples of observed input field currents, i (t). f
Fig. 5. Samples of output armature currents, i (t). a
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
621
Table 2 Estimated HMF-model parameter values of the plant Model hK
1
5
Number of runs 10 15
20
Case I: [0, ¹ ]"[0, 2.56sec], t "40]10~3 sec, q "1.47 Nm st l hK #hK 1 4 hK 1 hK hK 2 5 hK hK 1 4 hK hK 2 6 hK 3 hK hK 3 4 Fig. 6. Samples of observed motor speeds of the dc motor plant, n(t).
main procedure is the estimation of the HMF-model parameters using the weighted least squares, with a weighting factor inversely proportional to frequency. Finally, the two main procedures above are performed for a required number of runs, and then the average values of the parameters are calculated. Applying the aforementioned HMF method for case studies I and II gives the estimated HMF-model parameters and physical parameters shown in Tables 2 and 3, respectively for different numbers of runs. In the parameter-estimation study only (2M#1)"17 significant spectral data points were used. The physical parameters of the plant can be estimated using the following relations of the HMF-model parameters, i.e., armature inductance ¸K "1/hK [H]; armature resistance RK "hK /hK a 3 a 1 3 [)]; motor constant kK "(hK hK )/(hK hK ) [N m/A2]; total m 2 5 2 6 system moment of inertia JK "(hK hK )hK /(hK hK )2 2 5 3 2 6 [N m sec2], and viscous friction constant BK "(hK hK ) 1 4 (hK hK )hK /hK (hK hK )2 [N m sec]. 2 5 3 1 2 6 In general, it is not always possible to obtain the physical parameters directly/uniquely for a more complex system. It will depend on the mathematical relation between the physical system and its HMF-model parameters. Note that the main computations of the HMF method, given the input-output data, involve the approximation of the continuous-time signals’ Hartley transforms. These can be computed using a standard parabolic rule such as the extended Simpson’s rule for better approximation, given by
P
T
0
C
¹ x(t)cas(mu t)dt+ x(0) cas(0) 0 3N
A B A A B A D
#4
l¹ mu l¹ N~1 0 + x cas N N l"1, 3, 2
#2
l¹ mu l¹ N~2 0 + x cas N N l"2, 4, 2
#x(¹ ) cas(mu ¹ ) , 0
B B (19)
41.030 39.367 46.492 65.455 183.471 3.189 6.695
42.830 40.884 52.600 79.556 217.152 3.228 7.122
43.190 41.110 67.985 85.475 280.164 3.242 8.550
42.843 40.819 53.905 82.594 227.338 3.248 8.065
42.878 40.844 46.996 83.081 199.342 3.255 8.012
Case II: [0, ¹ ]"[0, 3.2 sec], t "50]10~3 sec, q "1.65 Nm st l hK #hK 1 4 hK 1 Kh hK 2 5 hK hK 1 4 Kh hK 2 6 hK 3 hK hK 3 4
38.849 37.346 41.882 56.124 176.008 3.045 6.018
43.685 41.981 54.028 71.529 200.881 3.309 7.092
43.127 41.385 50.959 72.084 208.535 3.241 7.438
43.463 41.798 47.676 69.618 200.103 3.243 7.238
43.575 41.807 50.834 73.908 208.016 3.223 7.602
JK [Nmsec2]
BK [Nmsec]
Table 3 Estimated physical parameters of the plant Number of runs
¸K a [H]
RK a [)]
kK m [Nm/A2]
Case I: [0, ¹ ]"[0, 2.56sec], t "40]10~3 sec, q "1.47 Nm st l 1 5 10 15 20
0.314 0.310 0.308 0.308 0.307
12.345 12.665 12.680 12.567 12.548
0.25340 0.24223 0.24266 0.23711 0.23576
0.00441 0.00360 0.00281 0.00339 0.00385
0.00732 0.00701 0.00584 0.00686 0.00783
Case II: [0, ¹ ]"[0, 3.2 sec], t "50]10~3 sec, q "1.65 Nm st l 1 5 10 15 20
0.328 0.302 0.308 0.308 0.310
12.265 12.687 12.769 12.889 12.971
0.23796 0.26896 0.24437 0.23826 0.24438
0.00412 0.00443 0.00380 0.00386 0.00379
0.00619 0.00755 0.00662 0.00643 0.00669
where ¹ is the window size and N is the number of samples within the fixed time interval (N"64). The other main computation of the methodology is the approximation of the HMF spectra of the input-output continuous signals, their derivatives, and the different form of products as given in Eqs. (5), (6), (7) and (15), respectively. Finally, comes the standard computation of the weighted least-squares algorithm to estimate the physical parameters of the nonlinear continuous-time dynamics. The main thing to be noted is that, in the computation of the parameter estimation, only M"8 significant spectral coefficients, or the mode number of
622
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626 Table 4 Estimated values of time constants Number of runs
1 5 10 15 20
Electrical time constant, tK []10~3 sec] e
Mechanical time constant, tK []10~3 sec] m
Case I
Case II
Case I
Case II
25.402 24.460 24.325 24.500 24.483
26.777 23.820 24.163 23.925 23.919
601.436 513.902 480.960 494.213 491.617
665.419 586.909 574.122 600.391 565.663
Fig. 7. HT and HMF spectrum of field current, i (t)"u(t). f
Fig. 8. HT and HMF spectrum of armature current, i (t)"y(t). a
Eq. (18b), have been used, even if the window size/record length has been discretized into N"64 subintervals as shown in Figs. 7 and 8 (the HT and HMF spectra versus frequency for the field current and armature current of case study I, respectively). Furthermore, based on the above estimated physical parameters, the electrical time constant, tK "¸K /RK " e a a 1/hK [sec] and the mechanical time constant, 1 tK "JK /BK "hK /(hK hK ) [sec] can be computed as shown in m 1 1 4 Table 4. The significant Hartley transform (HT) and the HMF spectra of the input-output signals for case I are shown in Figs. 7 and 8, respectively. The demonstration plant for variable-speed drives has been designed for general-purpose use, for instance, to understand the operating behavior of the dc machine (generator, motor) when separately excited, in shunt and compound connections, and with/without commutatingfield winding. This multipurpose machine has the following ratings: rated motor power 1.1 kW, voltage 130 V, current 11 A and speed 1500 rpm. According to the original data provided by the producer, the theoretical design values of the mechanical and electrical time constants of the plant, depending on the operating conditions, lie in the range of 3304t 4700 m []10~3 sec] and t +30 []10~3 sec], respectively. The e variation of the mechanical time constant can be achieved by varying the moment of inertia in the experi-
mental set-up, since the mechanical time constant will depend on the total system moment of inertia and the total system viscous friction constant. The load moment of inertia can be varied by changing the weight of the mass as shown in Fig. 2b. The experimental set-up has facilities to alter the mass, and as a result of this situation the mechanical time constant can vary correspondingly. Therefore, according to the operating conditions performed in the experiment, the values given in Table 4 are obtained, depending on different numbers of runs, and they lie within the range according to the mass taken in the experimental set-up. In general, if the minimum mass is coupled, the minimum mechanical time constant can be reached. On the other hand, if all the masses of the set-up are coupled, the maximum mechanical time constant can be reached. In addition, even if it is not possible to relate all values of the estimated physical parameters exactly to the design data of the general-purpose dcmachine set-up, one can clearly observe that the estimated values of the physical parameters are almost kept constant as the number of runs increases from 1 to 5, 10, 15 and 20, as shown in Tables 3 and 4 for both case studies. Moreover, the approximate design value of the armature inductance can be calculated, based on a method described in the literature (Langhoff and Raatz, 1977), and on dc-machine ratings. According to the experimental conditions involved here, an approximate design value of inductance ¸ "3¸ "3(º k)/(I N D )" a .*/ ar ar r a 0.26211 H can be taken into consideration, where º , ar I and N are the rated values of the dc-motor voltage, ar r current and speed, respectively; k is the compensation constant of the machine, k"1000; D "3JP /N and a r r P is the rated power of the motor. Furthermore, it may r be necessary to take into account that the estimated inductance ¸K +¸K #¸K , where ¸K is the motor a amot int amot armature inductance and ¸K is the inductance of the int operating effect or actuating mechanism. On the other hand, the estimated resistance RK +RK #RK #RK , a amot cab int where RK is the motor armature resistance including amot the commutating poles and compensating windings, RK is the resistance of the cables, and RK is the internal cab int resistance of the operating effect or actuating mechanism.
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Essentially, the experimental set-up is designed for multipurpose use, and was devised as a teaching, training aid and research set-up in the field of control engineering. As a result, some improvements and modifications have been implemented over time. Basically, the electrical time constant will depend on the armature resistance and inductance. In general, the deviation of the estimated electrical time constant from its approximate theoretical value may be due to the replacement of old cables, or replacement of actuating mechanisms, motor windings, compensating windings, etc. In general, the results show the performance of the HMF method in estimating the physical parameters of the nonlinear continuous-time dynamics of a thyristordriven dc-motor.
6. Conclusion In this paper, an application of the HMF-method to estimate the parameters of an experimental thyristordriven dc-motor laboratory set-up is carried out. The proposed method has given satisfactory results, and is promising for use in parameter estimation of physicallybased continuous-time nonlinear systems from their input-output data records, observed over a fixed time interval or window size, without the necessity for estimating unknown initial conditions. This recently developed methodology uses a frequency-dependent weighted leastsquares algorithm to identify the physical parameters. The practical implementation of the HMF method is facilitated by the CADACS computer-aided control-system design software. The armature terminal voltage and field motor circuit inputs are set from the PC of the real-time operating-system software of CADACS, and the input-output responses are displayed on the analog instruments, as well as on the PC screen via A/D converters. Then, the measured input-output responses are applied to an HMF algorithm of the dc-motor model to estimate the physical system parameters. The practical application of the HMF method and good results obtained here demonstrate the efficiency of this new approach. The proposed methodology is of considerable practical importance in problems of continuous-time parameter estimation in which input-output measurements are available on a finite interval of time. In general, the main advantage of the Hartley modulating functions method is that it is suitable for use in identifying the physically-based parameters of nonlinear continuoustime systems. Furthermore, the basic benefits of the HMF approach are that it is able to handle noisy data, and it does not need fixed/known initial conditions. In addition, the method employs integration-by-parts and the two-point boundary conditions of the modulating properties to transfer derivatives of data variables into smooth
623
moduating functions, in order to obtain an algebraic equation that is suitable for identification purposes. Generally, these features may make the modulating functions method desirable for use in several real processes. In the light of this encouraging result of the HMF method, further studies will continue on the computational considerations associated with the approach, and on how to implement this technique in the form of batch-scheme recursive estimation for nonlinear continuous-time systems. In principle, the HMF-method can be applied for real-time identification using a batch scheme recursive parameter-estimation technique. This can be implemented by shifting a fixed window size, or recorded batch data, one step forward at each sampling time. The main procedure of the batch scheme identification method is as follows: First, obtain an estimate based on initial batch of data; for instance one can assume the estimated parameters obtained by a single run shown in Tables 2, 3 or 4, then consider the new input-output sample data, and shift the fixed window size one step forward at each sampling time by appending a new column and deleting the first column of the data vector. Finally, the Hartley spectrums can be computed recursively, and the new estimate of the parameters obtained. This study of batchscheme recursive identification for nonlinear continuoustime systems is currently under investigation. Acknowledgements The first author thanks the German academic exchange service for providing the financial support for his stay at the Ruhr-University Bochum. He is also grateful to the Department of Electrical Engineering, Faculty of Technology, Addis Abada University, Addis Ababa, Ethiopia, for sanctioning leave of absence. References Langhoff, J., Raatz, E., 1977. Geregelte Gleichstromantriebe AEGTelefunken-Handbu¨cher. Band 19, Elitera-Verlag, Berlin. Hartley, R.V.L., 1942. A more symmetrical Fourier analysis applied to transmission problems. Proc. Inst. Radio Engrs., 144 —150. Jalili, S.A., Jordan, J.A., Mackie, R.D.L., 1992. Measurement of the parameters of all-pole transfer functions using shifted Hermite modulating functions. Automatica 38, 613— 617. Karanam, V.R., Frick, P.A., Mohler, R.R., 1978. Bilinear system identification by Walsh functions. IEEE Trans. Autom. Cont., AC-23, 709—713. Kung, F.C., Shih, D.H., 1986. Analysis and identification of Hammerstein model nonlinear delay systems using block-pulse function expansion. Int. J. Cont. 43(1), 139—147. Maletinsky, V., 1979. Identification of continuous dynamical systems with spline-type modulating functions method. 5th IFAC Symp. on Ident. and Syst. Par. Est. Germany, 275—281. Patra, A., 1989. General hybrid orthogonal functions and some applications in systems and control. Ph.D. Thesis, Department of Electrical Engineering, Indian Institute of Technology, Kharagpur-721302, India.
624
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
Patra, A., Unbehauen, H., 1995. Identification of a class of nonlinear continuous-time systems using Hartley modulating functions. Int. Journal of Control 62(2), 1431—1451. Pearson, A.E., 1992. Explicit parameter identification for a class of nonlinear I/O-differential operator models. Proc. of 31st IEEE CDC, Arizona, 3656—3660. Pearson, A.E., Lee, F.C., 1983. Time limited identification of continuous systems using trigonometric modulating functions. Proc. of 3rd Yale Workshop on Applic. of Adap. Syst., New Haven, CT, pp. 168—173. Pearson, A.E., Lee, F.C., 1985. On the identification of polynomial I-O-differential systems. IEEE Trans. Autom. Cont. AC-30(8), 778—782. Saha, D.C., Roy, B.K., Bapat, V.N., 1991. Some generalizations and further uses of the Poisson moment functionals (PMF ) approach to continuous-time model estimation. Proc. of the 9th IFAC Symp. on Ident. and Par. Estim., Budapest, Hungary, pp. 1324—1333. Schmid, C., 1992. Real-time control with CADACS-PC. In: Jamshidi, M., Herget, C.J. (Eds.) Recent advances in computer-aided control System Engineering, pp. 337—355, Elsevier Science Publishers, B.V., Amsterdam. Sinha, N.K., Rao, G.P., (eds.) 1991. Identification of continuous-time systems methodology and computer implementation. Kluwer Academic Publishers, Dordrecht, The Netherlands. Shinbrot, M., 1957, On the analysis of linear and nonlinear systems. Trans. ASME 79, 547—552. Subrahmanyam, A.V.B., Saha, D.C., Rao, G.P., 1995. Identification of continuous-time MIMO systems via time moments. ControlTheory and Advanced Technology 10(4), 1359 —1378. Tsang, K.M., Billings, S.A., 1994. Identification of continuous-time nonlinear systems using delayed state variable filters. Int. J. Cont. 60(2), 159—180. Unbehauen, H., 1993, CADACS: An integrated design facility for computer-aided design and analysis of control systems. In: Linkens, D.A. (Ed.), CAD for control systems, pp. 469— 487, Marcel Dekken Inc., New York. Unbehauen, H., Rao, G.P., 1987. Identification of continuous systems. North-Holland, Amsterdam. Unbehauen, H., Rao, G.P., 1990. Continuous-time approaches to system identification — a survey. Automatica 26(1), 23—35. Wang, M.L., Chen, K.S., Chou, C.K., 1984. Solution of integral equations via modified Laguerre polynomials. Int. J. Syst. Sci. 15(6), 661—672.
AB
Wang, C.H., Marleau, R.S., 1987. Recursive computational algorithms for the generalized block pulse operational matrices. Int. J. Cont. 45(1), 195—201. Young, P.C., 1981. Parameter estimation of continuous time models — a survey. Automatica 17(1), 23—39.
Appendix – Proofs [A]. Proof of Eq. (15). Let H ? H (i) be the m-th HMF spectral component of y u the i-th derivative of the product of two signals Mu(t)(i) y(t)N given by
P
T H ? H (i)" Mu(t)(i) y(t)N/ (t) dt. (20) u m y 0 Applying integration-by-parts i times while substituting the two end-point boundary conditions of / (t) given by m /(i) (0)"/(i) (¹ )"0, one can transfer the derivative of m m the signals to the derivative of the known modulating function, i.e.
P
AB
]cas((!1)i(n#m!j )u t) dt. 0
(21)
Now, replacing one of the signals by its infinite inverse Hartley series, say for u(t), = u(t)" + H (lu ) cas(lu t) u 0 0 l/~= results in
A BP
= n n in H ? H (i)" + H (lu ) + (!1)j (n#m!j )i ui cas@ y u 0 u 0 j 2 l/~= j/0
AB
P
T H ? H (i)"(!1)i u(t) y(t)M/(i) (t)N dt u y m 0 T n n " u(t) y(t)(!1)i + (!1)j j 0 j/0 ](n#m!j)iui (!1)i cas@(in/2) 0
A B
T 0
= n n in 1 " + H (lu ) + (!1)j (n#m!j )i ui cas@ ] u 0 0 j 2 2 l/~= j/0
(22)
y(t) cas(lu t) cas ((!1)i(n#m!j )u t) dt 0 0
G
P
T
y(t)cas ((l#(!1)i(n#m!j )) u t) dt 0
0
P P P
!
T
#
0 T
#
0 T 0
y(t)cas ((!l!(!1)i(n#m!j )) u t) dt 0 y(t)cas ((l!(!1)i(n#m!j )) u t) dt 0 y(t)cas ((!l#(!1)i(n#m!j )) u t) dt 0
since cas(m) cas(f)"1 [cas(m#f)!cas(!m!f)#cas(m!f)#cas(!m#f)]. 2
H
,
(23)
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
625
Replacing the integral terms of (23) by their corresponding Hartley transforms leads to the preceding equation of the form
AB
A B
= n n in H ? H (i)" + H (lu ) + (!1)j (n#m!j )i ui cas@ u y u 0 0 j 2 l/~= j/0
G
H
1 H ((l#(!1)i(n#m!j ))u )!H ((!l!(!1)i(n#m!j ))u ) y 0 y 0 ] . 2 #H ((l!(!1)i (n#m!j ))u )#H ((!l#(!1)i (n#m!j ))u ) y 0 y 0
(24)
For the case of Eq. (15), H ? H (2), substituting 2 in place of i gives the expression described in Eq. (15). This completes u y the proof for the general case. [B]. Proof of Eq. (7). Let the signals x (t) and x (t) be in the form of x (t)x(i) (t); then the HMF spectrum for such a product is given by 1 2 1 2
P
T
P
T
HM (0, x , xi) (mu )" 1 2 0
"
0
x (t) x(i) (t)/ (t) dt 1 2 m
AB
n n x (t) x(i) (t) + (!1)j 1 2 j 0 j/0
]cas((n#m!j )u t) dt. 0
(25)
Replacing x (t) by its infinite inverse Hartley series, i.e., x (t)"+ = H (lu ) cas (lu t), implies 1 1 0 l/~= x1 0
P
AB
T n n = HM x(0,, xi) (mu )" + Hx (lu ) x(i) (t) cas (lu t) + (!1)j cas((n#m!j ) u t) dt. 1 2 1 2 0 0 0 0 j 0 j/0 l/~=
(26)
Now, utilizing the identity given by cas(m) cas(f) "1 [cas(m#f)!cas(!m!f)#cas(m!f)#cas(!m#f)] into (26) 2 leads to
A BP
n = n HM (0, x , xi) (mu )" + Hx (lu ) + (!1)j 1 2 1 0 0 j l/~= j/0
T 0
x(i) (t) 1 Mcas((l#m#n!j ) u t)!cas ((!l!m!n#j) u t) 2 2 0 0
#cas((l!m!n!j )u t)#cas ((!l#m#n!j ) u t)N dt. 0 0
(27)
Applying integration-by-parts i times while substituting the two end-point boundary conditions of / (t) given by m /(i) (0)"/(i) (¹ )"0, one can transfer the derivative of the signals to the derivative of the known modulating m m function, i.e.
G
P
T
(!1)ix (t) cas(i)((l#m#n!j )u t) dt 2 0
0
P P P
T
(!1)ix (t) cas(i)((!l!m!n#j )u t) dt 2 0 n = n 1 0 HM (0, ] x , xi) (mu )" + Hx (lu ) + (!1)j 1 2 1 0 0 j 2 T l/~= j/0 # (!1)ix (t) cas(i)((l!m!n#j )u t) dt 2 0 0 T # (!1)ix (t) cas(i)((!l#m#n!j )u t) dt 2 0 0
AB
!
H
626
S. Daniel-Berhe, H. Unbehauen / Control Engineering Practice 6 (1998) 615—626
AB
A B
n in = n " + Hx (lu ) + (!1)j ui cas@ 1 0 0 j 2 l/~= j/0
1 ] 2
G
(l#m#n!j)i
P
T
0
x (t) cas((!1)i (l#m#n!j )u t) dt 2 0
P
!(!l!m!n#j)i
#(l!m!n#j )i
P
T
0
#(!l#m#n!j)i
T 0
x (t) cas((!1)i (!l!m!n#j )u t) dt 2 0
x (t) cas((!1)i (l!m!n#j )u t) dt 2 0
P
T
0
x (t) cas ((!1)i (!l#m#n!j )u t) dt 2 0
H
(28)
since cas(i)(qu t)"(!1)icas@(in/2)qiui cas((!1)iqu t). Replacement of the integral terms of Eq. (28) by their corre0 0 0 sponding Hartley transform symbols leads to the preceding equation of the form
G
(l#m#n!j )i Hx ((!1)i(l#m#n!j )u ) 2 0 n in = n 1 !(!l!m!n#j )iHx2 ((!1)i(!l!m!n#j )u0 ) HM x(0,, xi) (mu )" + Hx (lu ) + (!1)j ui cas@ ] 1 2 1 0 0 0 j 2 2 #(l!m!n#j )i Hx ((!1)i(l!m!n#j )u ) 2 l/~= j/0 0 #(!l#m#n!j )iHx ((!1)i(!l#m#n!j )u ) 2 0
AB
"Hx (mu ) ? HM x(i) (mu ). 1 2 0 0 This completes the proof of (7).
A B
H
(29)