- Email: [email protected]
Box–Jenkins identification revisited—Part III
Automatica 43 (2007) 868 – 875 www.elsevier.com/locate/automatica
Brief paper
Box–Jenkins identification revisited—Part III Multivariable systems夡 R. Pintelon a,∗ , J. Schoukens a , P. Guillaume b a Electrical Measurement Department, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium b Department of Mechanical Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
Received 31 October 2005; received in revised form 29 May 2006; accepted 6 November 2006 Available online 12 March 2007
Abstract Part I of this series of three papers handles the identification of single input single output Box–Jenkins models on arbitrary frequency grids in an open and closed loop setting. Part II discusses the computational aspects and illustrates the theory on simulations and a real life problem. This paper extends the results of Parts I and II to multiple input multiple output systems. Contrary to the classical time domain approach, the presented technique does not require symbolic calculus for multiple output polynomial Box–Jenkins models. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Multivariable; Parametric noise model; System identification; Feedback; Discrete-time; Continuous-time; Frequency domain
1. Introduction Recently (Ljung, 1993; McKelvey, 2002; Pintelon & Schoukens, 2005; Pintelon, Rolain, & Schoukens, 2005) a framework for identifying continuous-time and discrete-time Box–Jenkins (BJ) models in a given frequency band (part(s) imaginary axis or unit circle) has been developed for single input, single output (SISO) systems. The use of continuoustime noise models has been motivated in Pintelon, Schoukens, and Guillaume (2006), and an in depth analysis of the theoretical and computational aspects of maximum likelihood (ML) estimation in an open and closed loop setting has been given in Pintelon and Schoukens (2005) and Pintelon et al. (2005). The ML solution for multivariable systems in open loop can be found in McKelvey and Ljung (1997). This paper generalizes all these results to multivariable systems in closed loop, and is intended for those applications where one has to identify transfer function models in a given frequency band using the operational perturbations. Examples are modal analysis of bridges, buildings and airplanes subject to wind excitation, and 夡 This paper was not present at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor Brett Ninness under the direction of Editor Torsten Söderström. ∗ Corresponding author. Tel.: +32 2 629 29 44; fax: +32 2 629 28 50. E-mail address: [email protected] (R. Pintelon).
0005-1098/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.11.007
modelling in the process industry. The following important issues are handled: (i) maximum likelihood (ML) Box–Jenkins (BJ) identification in feedback; (ii) discussion of the asymptotic properties (consistency, asymptotic efficiency, asymptotic normality, robustness, etc.) of the ML estimator; (iii) numerical stable Newton–Gauss minimization of the ML cost function; (iv) generation of starting values; (v) calculation of the Cramér–Rao lower bound. Finally, we also report the results of a real life application: flight flutter analysis. 2. BJ identification in feedback 2.1. Closed loop framework Consider the closed loop setup of Fig. 1 Y (k) = G(k )U (k) + V (k), U (k) = R(k) − M(k )Y (k), V (k) = H (k )E(k)
(1)
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
M (Ω)
R(k)
U(k)
Y(k) G (Ω) V(k) = H (Ω)E(k)
Fig. 1. Closed loop setup where = z−1 for discrete-time systems and = s for continuous-time systems.
with R(k), U (k), and Y (k) the nu × 1 reference, nu × 1 input, and ny × 1 output spectra, respectively; G(), H (), and M() the ny × nu plant, ny × ny noise, and nu × ny controller transfer functions, respectively; E(k) the ny × 1 driving white noise source; k = zk−1 = exp(−j2fk /fs ), with fs the sampling frequency, for discrete-time systems; k = sk = j2fk for continuous-time systems; and k the frequency index. Assmption 1. 1. The input/output data are generated by (1). 2. The controller transfer function M() is known. 3. E(k) is independent (over the frequency index k), circular complex (E{E(k)E T (k)} = 0, with E{ } the expected value, and superscript T the transpose operator) normally distributed noise with zero mean (E{E(k)} = 0) and covariance matrix E{E(k)E H (k)} = (superscript H is the complex conjugate transpose operator). 4. The ny × ny covariance matrix is real valued. 5. E(k) is independent of the reference signal R(k). 6. The input U (k) and output Y (k) of the plant are observed (measured) without errors. 7. The frequency domain data U (k), Y (k) are available at frequencies fk , k ∈ K = {1, 2, . . . , F }. Discussion: • If the raw data are time domain signals, then Assumption 1.1 is asymptotically (for the number of time domain samples N going to infinity) valid after a discrete Fourier transform (DFT) 1 X(k) = √ N
N−1
x(t)zk−t
(2)
t=0
with x = r, u, . . . and X = R, U, . . . (Pintelon & Schoukens, 2001), and the frequency domain data in Assumption 1.7 are available at DFT frequencies fk =kf s /N , k ∈ K, where K ⊆ {0, 1, 2, . . . , N/2} and ord(K) = F . For finite N Eq. (1) is extended with transients terms that take into account the initial and final conditions of the experiment (see the discussion of Section 2.2). These transient terms vanish as O(N −1/2 ) w.r.t. the main contributions (Pintelon & Schoukens, 2001).
869
• Assumption 1.3 is valid after a DFT for normally distributed white noise e(t). It is asymptotically (N → ∞) valid after a DFT for a general class of non-Gaussian time domain errors (see Brillinger, 1981; Pintelon & Schoukens, 2001). • The delays between the input noise channels are included in the noise transfer function H (); otherwise the covariance matrix of E(k) would be complex valued and frequency dependent, and Assumption 1.4 would be violated. • In the open loop case (M = 0 in Fig. 1) the output noise V (k) may also contain output measurement errors. 2.2. Gaussian ML cost function Under Assumption 1, Y (k) is independent (over k), circular complex normally distributed. To construct the likelihood function (U (k) is known exactly) it is sufficient to calculate the mean and covariance of Y (k) given the model parameters and the covariance matrix of the driving white noise source E(k). Since the process noise V (k) is correlated with the input of the plant U (k), and since it is independent of the reference signal R(k), the expected values in the mean and covariance calculation of Y (k) should be conditioned on R(k). The latter is known since the controller is known and since Y (k) and U (k) are observed without errors. Using (1) we find that Y (k) = (Iny + G(k )M(k ))−1 G(k )R(k) + (Iny + G(k )M(k ))−1 H (k )E(k)
(3)
and E{Y (k)|R(k), , } = (Iny + G(k , )M0 (k ))−1 G(k , )R(k) ≡ Y (k, ), Cov(Y (k)|R(k), , ) = S(k , )S H (k , ) ≡ CY (k) (), (4) S(k , ) = (Iny + G(k , )M0 (k ))−1 H (k , ). Hence, within a -independent additive constant, the negative Gaussian log-likelihood function is given by VNLL (, , Z) = log det(CY (k) ()) + k∈K
k∈K
(Y (k) − Y (k, ))H CY−1 (k) ()(Y (k) − Y (k, )),
(5)
where Z is a vector containing the input/output data (see Pintelon & Schoukens, 2001, p. 437, Eq. (14-14)). At DC (f = 0) and Nyquist (f = fs /2) the sums in (5) are multiplied by 21 . Using S −1 (k , )(Y (k) − Y (k, )) = H −1 (k , )(Y (k) − G(k , )U (k)) ≡ (k , ), log det(CY (k) ()) = log | det(S(k , ))|2 + log det(). (6)
870
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
Eq. (5) can be written as
Discussion:
VNLL (, , Z) = log | det(S(k , ))|2 k∈K
+ F log det() +
• If the raw data are time domain signals, then the model Y (k) = G(k , )U (k) + H (k , )E(k)
H (k , )−1 (k , ).
(7)
k∈K
For the open loop case (M0 =0 ⇒ S =H ) Eq. (7) divided by F reduces to Eq. (9) in McKelvey and Ljung (1997). Minimizing the real part of (7) w.r.t. (VNLL is a real number) gives 1 H (k , ) (k , ) , (8) () = Re F k∈K
k∈K
1 + F log det Re (k , )H (k , ) F k∈K
+ F ny .
(9)
For the open loop case (S = H ) Eq. (9) divided by F reduces to Eq. (12) in McKelvey and Ljung (1997), except that in the latter () is not constrained to be real and, hence, is sub-optimal. Dividing (9) by F and taking the exponential function gives, within a -independent multiplicative constant, 1 2 H (k , ) (k , ) , VF (, Z) = |hF ()| det Re F hF () = exp
is asymptotically (number of time domain samples N → ∞) valid. To improve the finite sample behaviour of the estimates the plant and noise transient terms are added in (11) (k , ) = H −1 (k , )(Y (k) − G(k , )U (k) − TG (k , ) − TH (k , )),
(13)
where, for example,
where Re(x) stands for the real part of x. Combining (7) and (8) gives log | det(S(k , ))|2 VNLL (, Z) =
(12)
k∈K
1 log det(S(k , )) . F
G = A−1 B,
TG = A−1 IG ,
H = D −1 C,
TH = D −1 IH
(14)
with A, B, C, D, IG and IH matrix polynomials of the appropriate matrix dimensions (see Pintelon & Schoukens, 2001). The coefficients of IG and IH decrease as an O(N −1/2 ) and, hence, for N sufficiently large the transients terms can be neglected in (13). • If the following conditions are satisfied: (i) k = zk−1 with zk = exp(j2k/N ) and k ∈ K = {0, 1, . . . , N − 1}, (ii) det(S(z−1 , )) has all its roots within the unit circle, and (iii) lim det(S(z−1 , )) = 1,
(15)
z→∞
then (11) simplifies to VN (, Z) = det
(10)
N−1 1 −1 (zk , )H (zk−1 , ) N k=0
k∈K
+ O(|max |N /N )
Taking the ny th root of |hF the factor can be brought into the determinant giving 1 H (k , ) (k , ) , VF (, Z) = det Re F ()|2 ,
k∈K
(k , ) = gF ()(k , ), 1 log det(S(k , )) , gF () = exp F ny
(11)
k∈K
where S(k , ) and (k , ) are defined in (4) and (6). This results in the following theorem:
(16)
with max the dominant pole of log det(S(z−1 , )) (proof: first note that the sum over the full unit circle is real, and next follow the lines of Appendix A.3 in Pintelon & Schoukens, 2005). Using (13), the first term in the right-hand side of (16) equals exactly the classical time domain result (see Ljung, 1999, p. 219, Eq. (7.96)) det
N−1 1 −1 (zk , )H (zk−1 , ) N k=0
= det
N−1 1 (t, )T (t, ) N
(17)
t=0
Theorem 1 (ML cost function). Under Assumption 1, the Gaussian ML cost function is given by (11), where at DC (f = 0) and Nyquist (f = fS /2) the summands are multiplied by 21 .
with (t, )=H −1 (q, )(y(t)−G(q, )u(t)) and q the backward shift operator (qx(t) = x(t − 1)).
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
2.3. ML estimator Consider, for example, a common denominator representation of the plant and noise models nb Br r B(, ) G(, ) = = r=0 na r , A(, ) r=0 ar nig Igr r IG (, ) TG (, ) = = r=0 na r , A(, ) r=0 ar nc Cr r C(, ) H (, ) = = r=0 nd r , D(, ) r=0 dr nih Ihr r IH (, ) = r=0 TH (, ) = nd r D(, ) r=0 dr
(18)
with A and D polynomials, IG and IH vector polynomials of size ny × 1, and B and C, matrix polynomials of size ny × nu and ny × ny , respectively. The parametric models (18), where contains all the numerator and denominator polynomials coefficients, are overparametrized. Indeed, multiplying the numerators and denominators in (18) by a non-zero real number leaves the parametric models unchanged. This imposes a constraint on the coefficients of the A and D polynomials. Further, the fact that E(k) in (3) is not observed, and that H (k , )E(k) = H (k , )T T −1 E(k) for any regular ny ×ny matrix T , imposes n2y additional constraints on the coefficients of the matrix polynomial C. The following parameter constraints are often used: z-domain
s-domain
ARMAX, BJ, OE :
a0 = 1
or ana = 1
ARMA, ARMAX, BJ :
C0 = Iny
or Cnc = Iny
BJ :
d0 = 1
or dnd = 1,
(19)
where for CT models the constraints are imposed on the normalized model parameters; other choices are, however, possible (see Section 3.1). Theorem 2 (ML estimator). Under Assumption 1 the ML ˆ estimator (Z) of model structure (12), where G(, ) and H (, ) are, for example, parametrized as in (18), and with a vector containing the plant and noise transfer function coefficients, minimizes (11) subject to the constraints (19). Discussion • It can easily be verified that the cost function VF (, Z) (11) contains exactly the same parameter ambiguities as the plant and noise models in (18). Using the results of Pintelon, Schoukens, Vandersteen, and Rolain (1999), it follows that ˆ and H (, ) ˆ 1/2 (), ˆ with the estimated models G(, ) ˆ the minimizer of (11), are independent of the particular parameter constraint(s) chosen. • Using H −1 (, ) = D(, ) adj (C(, ))/ det(C(, )), with adj(C) the adjoint matrix of C, it can be seen that the cost function VF (, Z) (11) depends on D(k , ) and
871
det(C(, )) via |D(k , )|2 and det(C(k , ))|2 only. It shows that the cost function remains the same when the roots of D(, ) that are not in common with A(, ) and the roots of det(C(, )) are mirrored w.r.t. the stability border ( → −, ¯ where x¯ denotes the complex conjugate of x, for continuous-time systems, and → −1 for discrete-time systems). This is not true for the zeroes of the polynomial matrix adj(C(, )), unless the MIMO problem reduces to ny independent SISO problems. Note that neither H nor H −1 should be stable for the calculations in the frequency domain. For time domain calculations H (D) and H −1 (det(C)) can be stabilized without modifying the statistical properties of the estimate. The stabilization of H −1 , however, requires symbolic calculations for multiple output systems. This is the bottle neck of the time domain approach. It can only be avoided via a specially parametrized state space model that excludes common dynamics in the plant and noise models. • Under some additional standard assumptions on the true model, the model set and the cost function (see Pintelon & Schoukens, 2005, Assumptions 9–11) it can be shown that the ML estimator in Theorem 2 is strongly consistent, asymptotically normally distributed, and asymptotically efficient (see Caines, 1988). • Following exactly the same lines of Pintelon and Schoukens (2005) the robustness of the ML estimator ˆ in Theorem 2 can be shown. The estimate (Z) remains strongly consistent and asymptotically normally distributed for non-Gaussian, mixing frequency domain noise V (k) (the mixing condition allows for a correlation over the frequency). Moreover, using the definition 1 VF () = det Re E (k , ) H (k , ) F k∈K
(20) the general expressions for the convergence rate and the ˆ asymptotic covariance matrix of (Z) remain valid (see Pintelon & Schoukens, 2005, Theorem 3, Properties 1–5). 3. Minimization of the ML cost function 3.1. Newton–Gauss algorithm The minimizer of (11) ()) with VF (, Z) = det( F 1 H () = Re (k , ) (k , ) F
(21)
k=1
is found as the solution of jVF (, Z)/ji = 0, i = 1, 2, . . . , dim(). Since j log det(X)/jX = X −T it follows that
jX j det(X) , (22) = det(X) trace X −1 ji ji
872
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
where X = (). Hence, jVF (, Z) =0 ji
⇒
j () −1 =0 trace () ji
(23)
(JL ())re = −(L ())re ,
for i = 1, 2, . . . , dim(). Using
j () 1 j (k , ) H = Re 2 herm , (k , ) ji F ji k
−T /2 () −1/2 (), −1 () = Y herm (X)Y H = herm (Y XY H ), trace(herm()) = Re(trace()),
(24)
with herm(X) = (X + X H )/2, Eq. (23) can be written as j () −1 trace () ji 2 j (k , ) = Re trace −1/2 () F ji k
na
ar2
+
r=0
ARMA, ARMAX, BJ :
nb
trace(BrT Br ) = 1
r=0 nc
CrT Cr = Iny
r=0
BJ :
, ) j ( 2 k ( −1/2 () (k , ))H = Re −1/2 () F ji
nd
dr2 = 1
(31)
r=0
because they are numerically more robust than (19) (constraints (19) result in an ill conditioned set of Eqs. (30) if, for example, in the z-domain the true value of a0 is zero).
k
2 j () H = Re (L() , ()) L() F ji
(30)
where ( )re puts the real and imaginary parts on top of each other, via a QR-factorization or an SVD decomposition (Golub & Van Loan, 1996). The iterative Newton–Gauss algorithm is implemented as follows. First we calculate (JL ())re w.r.t. all model parameters (constraints (19) are not imposed), next we solve (30) using the pseudo-inverse of (JL ())re (the dimension of the null space of (JL ())re equals one for OE, two for ARMA(X), and three for BJ), and finally we impose the following constraints on the updated model parameters + : ARMAX, BJ, OE :
× ( −1/2 () (k , ))H
(), and J()=j ()/j, with JL ()=L()J(), L ()=L() and where (), L() and () are defined in (21) and (26), respectively. A numerical stable implementation consists in solving
(25) 3.2. Starting values
where −1/2
L() = IF ⊗
()
and
() = [ T (1 , ), . . . , T (k , ), . . . , T (F , )]T
(26)
with ⊗ the Kronecker product. From (25) we conclude that the minimizer of (21) is found as the solution of the following set of equations:
j () Re (L() ())H L() = 0. (27) j The derivative of the left-hand side of (27) w.r.t. equals
j j () H Re (L() ()) L() j j
j () H j () = Re L() L() + O( ). (28) j j Neglecting the O( ) terms in (28), the Newton–Raphson algorithm for finding the solution of (27) reduces to the Newton–Gauss method (Fletcher, 1991) Re(JLH ()JL ()) = −Re(JLH ()L ())
(29)
First starting values for the plant model are calculated by minimizing the following weighted non-linear least squares cost function: 1 H (k , )Cˆ Y −1 (k , ) subject to constraints(31) (32) F k∈K
w.r.t. , where Cˆ Y is a non-parametric estimate of the noise model (see Appendix A), and where (32) is started from a weighted linear least squares version of (32). Next starting values for the noise model are obtained in three steps: (i) calcuˆ (k), lation of the output residuals Vˆ (k) = Y (k) − G(k , )U with ˆ the minimizer of (32); (ii) estimate of the denominator polynomial D via a SISO-ARMA fit of ⎛
⎞1/2 ny |Vˆ[p] (k)|2 ⎝ ⎠ var(Vˆ[p] )
(33)
p=1
with var(Vˆ[p] ) = (1/F ) k∈K |Vˆ[p] (k)|2 ; (iii) estimate of the numerator polynomial matrix C via a MIMO-MA fit of the ˆ Vˆ (k), with D(k , ) ˆ the solution of filtered residuals D(k , ) step (ii).
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
4. Calculation of the Cramér–Rao lower bound The ultimate goal is to calculate uncertainty bounds on invariants of model structure (12) starting from the Cramér–Rao lower bound on the estimated model parameters. Examples of invariants are the plant transfer function G(, ) and its poles and zeroes, and the noise power spectrum H (, )H H (, ). Note that the uncertainty calculation of H H H only requires the uncertainty of H 1/2 and not that of H and separately. Hence, the calculation of the Cramér–Rao lower bound can be simplified significantly by imposing the constraint = Iny . This is done by right multiplication of the estimated matrix coˆ giving C ˆ Since H 1/2 r = Cˆ r 1/2 (). efficients Cˆ r by 1/2 () can still be right multiplied by an orthogonal matrix T without changing the noise power spectrum H H H , the symmetry of r should also be imposed. This one of the matrix coefficients C is done as follows: first one matrix coefficient is decomposed 0 = U0 0 V T , and next all in singular values, for example, C 0 r are right multiplied by T0 = V0 U T . Summatrix coefficients C 0 marized the Cramér–Rao lower bound of the model parameters subject to the constraints: z-domain
s-domain
ARMA, ARMAX, BJ, OE :
= Iny
or = Iny
ARMAX, BJ, OE :
a0 = 1
or ana = 1
ARMA, ARMAX, BJ :
C0 = C0T
or
BJ :
d0 = 1
or dnd = 1
(34)
Cnc = CnTc
Using constraints (34), the negative log-likelihood function is, within a -independent constant, given by (7) with = Iny . Applying the definition of the Cramér–Rao lower bound CR (0 ). CR−1 (0 ) = j2 E{VNLL (, Iny , Z)}/j20
CR−1 (0 ) = 2(J1 )Tre (J1 )re + (J2 )Tre (J2 )re ,
vec(WF )
1/2
Vk = S0−1 (k )S(k , )S H (k , )S0−H (k ), where vec( ) stacks the columns on top of each other, 1/2 X0 (k ) = X(k , 0 ), and SR (k) is a square root of H E{R(k)R (k)}. If the expected value in (35) is taken w.r.t. to 1/2 the disturbing noise E(k) only, then (36) with SR (k) = R(k) is the Cramér–Rao lower bound for the particular realization of the reference signal R(k). Finally, the Cramér–Rao lower bound on any model parameter related quantity f () is obtained via
20
100
-40 4
5
6
7
0 -100 -200
8
4
40
200
20
100 Phase (°)
Amplitude (dB)
Frequency (Hz)
0 -20 -40 4
6 Frequency (Hz)
8
vec(VF ),
Wk = H0 (k )G(k , )(Inu + M0 (k )G0 (k ))−1 SR (k),
200
-20
(36)
where ( )re stacks the real and imaginary parts on top of each other. J1 and J2 are given by ⎡ ⎡ ⎤ ⎤ vec(W1 ) vec(V1 ) ⎢ ⎢ ⎥ ⎥ vec(W2 ) ⎥ vec(V2 ) ⎥ j ⎢ j ⎢ ⎢ ⎢ ⎥ ⎥ and, J1 = ⎥ , J2 = j ⎢ ⎥ j0 ⎢ 0 ⎣ ... ⎦ ⎣ ... ⎦
40
0
(35)
we find after some calculations
Phase (°)
Amplitude (dB)
will be calculated.
873
5
6 7 Frequency (Hz)
8
6 Frequency (Hz)
8
0 -100 -200
4
Fig. 2. Flight flutter data: measured frequency response function (gray bullets), plant model (bold solid line), and standard deviation plant model (bold dashed line). Left column: amplitude, right column: phase.
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
20
20
0
0
Amplitude (dB)
Amplitude (dB)
874
-20 -40
-20 -40 -60
-60 4
5
6
7
8
4
5
4
5
20
20
0
0
Amplitude (dB)
Amplitude (dB)
Frequency (Hz)
-20 -40 -60
6 7 Frequency (Hz)
8
-20 -40 -60
4
5
6
7
8
Frequency (Hz)
6
7
8
Frequency (Hz)
Fig. 3. Flight flutter data: measured noise power spectrum (gray bullets), modelled noise power spectrum (bold solid line), and standard deviation modelled noise power spectrum (bold dashed line).
CR(f (0 )) = f (0 )CR(0 )f H (0 ), where denotes the derivative w.r.t. (see Caines, 1988; Pintelon & Schoukens, 2001). 5. Measurement example: flight flutter analysis Although the asymptotic properties (consistency, asymptotic covariance matrix, robustness, etc.) of the ML estimator has been verified on numerous simulation examples, we prefer to report the results of a real life application: flight flutter analysis. To excite the air plane during flight, a perturbation signal is injected in the control loop of the flap mechanism at the tip side of the right wing. The angle perturbation of the flap is used as a measure of the applied force, and the acceleration is measured at two different positions on the left wing. Beside the applied perturbation the air plane is also excited during flight by the natural turbulence. The resulting turbulent forces acting on the air plane cannot be measured and are assumed to be white in the frequency band of interest. Since the input/output signals are lowpass filtered before sampling, a continuous-time plant and noise model is the natural choice. The signals are measured during about 128 s at the sampling rate fs = 256 Hz, giving N = 32 768 data points per channel. Fig. 2 shows the 2 × 1 force to acceleration frequency response function (FRF) in the band [4.00, 7.99 Hz] (DTF lines k = 513, 514, . . . , 1023 ⇒ F = 511). From Figs. 2 and 3 it can be seen that a continuoustime ARMAX model ((18) with = s and D = A) of order na = 8, nb = 8, nig = 9, and nc = 8 explains the measurements
very well. The non-parametric FRF and noise power spectra shown in these figures (gray bullets) are calculated as ˆ ˆ ˆ k ) = Y (k) − TG (sk , ) − TH (sk , ) , G(s U (k) ˆ (k) − TG (sk , ) ˆ Vˆ (k)Vˆ H (k) with Vˆ (k) = Y (k) − G(sk , )U ˆ − TH (sk , ),
(37)
where ˆ is the ML estimate. 6. Conclusion It has been shown that all results for SISO BJ identification on parts of the imaginary axis (s-domain) or unit circle (zdomain) gracefully extend to the multivariable case. Although the computational aspects of the ML estimator (Newton–Gauss minimization, parameter constraints, starting values, Cramér–Rao lower bound) have been discussed for one particular parametrization of the transfer function models (common denominator), exactly the same lines can be followed for other parametrization (for example, left matrix fraction description). Acknowledgement This work is sponsored by the Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (GOA-ILiNoS)
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
and the Belgian Program on Interuniversity Poles of Attraction initiated by the Belgian State, Prime Minister’s Office, Science Policy programming (IUAP V/22). Appendix A. Non-parametric noise model The data records of N input/output samples are split in M blocks of equal length N/M , and the input/output DFT spectra U [r] (k), Y [r] (k) of each block are calculated (r =1, 2, . . . , N/M ). To reduce the leakage effect the spectra are differenced (Schoukens, Rolain, & Pintelon, 2006) YD[r] (k + 0.5) = Y [r] (k + 1) − Y [r] (k), UD[r] (k + 0.5) = U [r] (k + 1) − U [r] (k)
(38)
and M is kept as small as possible. Next, the noise covariance matrix is estimated as
ˆH Cˆ y (k) = (SˆYD YD (k) − SˆYD UD (k)Sˆ −1 (39) UD UD (k)S YD UD ), 2 where =M/(M −nu ) is a bias correction factor, 21 accounts for [r] [r] H the differencing (38), and where SˆXY = M −1 M r=1 X (Y ) . From (39) it follows that the minimal value of M equals nu +1. Finally, in order to guarantee the existence of the second order ˆ moments of Cˆ −1 Y (k) for a given input, CY (k) is averaged over Q neighbouring frequencies such that Q = (ny + 3)/(nu + 1) (Mahata, Pintelon, & Schoukens, 2004). To reduce the variability of the estimate we mostly choose Q > (ny + 3)/(nu + 1). References Brillinger, D. R. (1981). Time series: data analysis and theory. New York: McGraw-Hill. Caines, P. E. (1988). Linear stochastic systems. New York: Wiley. Fletcher, R. (1991). Practical methods of optimization. 2nd ed., New York: Wiley. Golub, G. H., & Van Loan, C. F. (1996). Matrix computations. 3rd ed., Baltimore: John Hopkins University Press. Ljung, L., (1993). Some results on identifying linear systems using frequency domain data. In Proceedings 32nd IEEE conference on decision and control, (Vol. 4, pp. 3534–3538) San Antonio, Texas, USA, December 15–17. Ljung, L. (1999). System identification: theory for the users. Upper Saddle River (USA): Prentice-Hall. Mahata, K., Pintelon, R., & Schoukens, J., (2004). Frequency domain identification using non-parametric noise models. In Proceedings of the 43rd IEEE conference on decision and control, (Vol. 1, pp. 821–826) Atlantis, Bahamas, December 14–17. McKelvey, T. (2002). Frequency domain identification methods. Circuits and Systems Signal Processing, 21(1), 39–55.
875
McKelvey, T., & Ljung, L., (1997). Frequency domain maximum likelihood identification. In 11th IFAC symposium on system identification, (Vol. 4, pp. 1741–1746) Kitakyushu, Japan, 8–11, July, preprints. Pintelon, R., Rolain, Y., & Schoukens, J. (2005). Box–Jenkins identification revisited—Part II: Applications. Automatica, 42(1), 77–84. Pintelon, R., & Schoukens, J. (2001). System identification: A frequency domain approach. New York (USA): IEEE Press. Pintelon, R., & Schoukens, J. (2005). Box–Jenkins identification revisited—Part I: Theory. Automatica, 42(1), 63–75. Pintelon, R., Schoukens, J., & Guillaume, P. (2006). Continuous-time noise modelling from sampled data. IEEE Transactions on Instrumentation and Measurement, 55(6), 2253–2258. Pintelon, R., Schoukens, J., Vandersteen, G., & Rolain, Y. (1999). Identification of invariants of (over)parametrized models: Finite sample results. IEEE Transactions on Automatic Control, AC-44(5), 1073–1077. Schoukens, J., Rolain, Y., & Pintelon, R. (2006). Leakage reduction in frequency response function measurements. IEEE Transactions on Instrumentation and Measurement, 55(6), 2286–2291. Rik Pintelon was born in Gent, Belgium, on December 4, 1959. He received the degree of electrical engineer (burgerlijk ingenieur) in July 1982, the degree of doctor in applied sciences in January 1988, and the qualification to teach at university level (geaggregeerde voor het hoger onderwijs) in April 1994, all from the Vrije Universiteit Brussel (VUB), Brussels, Belgium. From October 1982 till September 2000 he was a researcher of the Fund for Scientific Research—Flanders at the VUB. Since October 2000 he is professor at the VUB in the Electrical Measurement Department (ELEC). His main research interests are in the field of parameter estimation/system identification, and signal processing. Johan Schoukens is born in Belgium in 1957. He received the degree of engineer in 1980 and the degree of doctor in applied sciences in 1985, both from the Vrije Universiteit Brussel, Brussels, Belgium. He is presently professor at the Vrije Universiteit Brussel. The prime factors of his interest are in the field of system identification for linear and nonlinear systems, and growing tomatoes in his green house.
Patrick Guillaume was born in Anderlecht, Belgium, in 1963. He received the degree in mechanical-electrotechnical engineering in 1987 from the Vrije Universiteit Brussel, Belgium. In 1987 he joined the Department of Fundamental Electricity and Intrumentation (ELEC) of the Vrije Universiteit Brussel and received the Ph.D. degree in 1992. In 1996 he joined the Department of Mechanical Engineering (MECH) where he is currently lecturer. His main research interests are situated in the field of system identification, signal processing and experimental modal analysis.
Brief paper
Box–Jenkins identification revisited—Part III Multivariable systems夡 R. Pintelon a,∗ , J. Schoukens a , P. Guillaume b a Electrical Measurement Department, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium b Department of Mechanical Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
Received 31 October 2005; received in revised form 29 May 2006; accepted 6 November 2006 Available online 12 March 2007
Abstract Part I of this series of three papers handles the identification of single input single output Box–Jenkins models on arbitrary frequency grids in an open and closed loop setting. Part II discusses the computational aspects and illustrates the theory on simulations and a real life problem. This paper extends the results of Parts I and II to multiple input multiple output systems. Contrary to the classical time domain approach, the presented technique does not require symbolic calculus for multiple output polynomial Box–Jenkins models. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Multivariable; Parametric noise model; System identification; Feedback; Discrete-time; Continuous-time; Frequency domain
1. Introduction Recently (Ljung, 1993; McKelvey, 2002; Pintelon & Schoukens, 2005; Pintelon, Rolain, & Schoukens, 2005) a framework for identifying continuous-time and discrete-time Box–Jenkins (BJ) models in a given frequency band (part(s) imaginary axis or unit circle) has been developed for single input, single output (SISO) systems. The use of continuoustime noise models has been motivated in Pintelon, Schoukens, and Guillaume (2006), and an in depth analysis of the theoretical and computational aspects of maximum likelihood (ML) estimation in an open and closed loop setting has been given in Pintelon and Schoukens (2005) and Pintelon et al. (2005). The ML solution for multivariable systems in open loop can be found in McKelvey and Ljung (1997). This paper generalizes all these results to multivariable systems in closed loop, and is intended for those applications where one has to identify transfer function models in a given frequency band using the operational perturbations. Examples are modal analysis of bridges, buildings and airplanes subject to wind excitation, and 夡 This paper was not present at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor Brett Ninness under the direction of Editor Torsten Söderström. ∗ Corresponding author. Tel.: +32 2 629 29 44; fax: +32 2 629 28 50. E-mail address: [email protected] (R. Pintelon).
0005-1098/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.11.007
modelling in the process industry. The following important issues are handled: (i) maximum likelihood (ML) Box–Jenkins (BJ) identification in feedback; (ii) discussion of the asymptotic properties (consistency, asymptotic efficiency, asymptotic normality, robustness, etc.) of the ML estimator; (iii) numerical stable Newton–Gauss minimization of the ML cost function; (iv) generation of starting values; (v) calculation of the Cramér–Rao lower bound. Finally, we also report the results of a real life application: flight flutter analysis. 2. BJ identification in feedback 2.1. Closed loop framework Consider the closed loop setup of Fig. 1 Y (k) = G(k )U (k) + V (k), U (k) = R(k) − M(k )Y (k), V (k) = H (k )E(k)
(1)
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
M (Ω)
R(k)
U(k)
Y(k) G (Ω) V(k) = H (Ω)E(k)
Fig. 1. Closed loop setup where = z−1 for discrete-time systems and = s for continuous-time systems.
with R(k), U (k), and Y (k) the nu × 1 reference, nu × 1 input, and ny × 1 output spectra, respectively; G(), H (), and M() the ny × nu plant, ny × ny noise, and nu × ny controller transfer functions, respectively; E(k) the ny × 1 driving white noise source; k = zk−1 = exp(−j2fk /fs ), with fs the sampling frequency, for discrete-time systems; k = sk = j2fk for continuous-time systems; and k the frequency index. Assmption 1. 1. The input/output data are generated by (1). 2. The controller transfer function M() is known. 3. E(k) is independent (over the frequency index k), circular complex (E{E(k)E T (k)} = 0, with E{ } the expected value, and superscript T the transpose operator) normally distributed noise with zero mean (E{E(k)} = 0) and covariance matrix E{E(k)E H (k)} = (superscript H is the complex conjugate transpose operator). 4. The ny × ny covariance matrix is real valued. 5. E(k) is independent of the reference signal R(k). 6. The input U (k) and output Y (k) of the plant are observed (measured) without errors. 7. The frequency domain data U (k), Y (k) are available at frequencies fk , k ∈ K = {1, 2, . . . , F }. Discussion: • If the raw data are time domain signals, then Assumption 1.1 is asymptotically (for the number of time domain samples N going to infinity) valid after a discrete Fourier transform (DFT) 1 X(k) = √ N
N−1
x(t)zk−t
(2)
t=0
with x = r, u, . . . and X = R, U, . . . (Pintelon & Schoukens, 2001), and the frequency domain data in Assumption 1.7 are available at DFT frequencies fk =kf s /N , k ∈ K, where K ⊆ {0, 1, 2, . . . , N/2} and ord(K) = F . For finite N Eq. (1) is extended with transients terms that take into account the initial and final conditions of the experiment (see the discussion of Section 2.2). These transient terms vanish as O(N −1/2 ) w.r.t. the main contributions (Pintelon & Schoukens, 2001).
869
• Assumption 1.3 is valid after a DFT for normally distributed white noise e(t). It is asymptotically (N → ∞) valid after a DFT for a general class of non-Gaussian time domain errors (see Brillinger, 1981; Pintelon & Schoukens, 2001). • The delays between the input noise channels are included in the noise transfer function H (); otherwise the covariance matrix of E(k) would be complex valued and frequency dependent, and Assumption 1.4 would be violated. • In the open loop case (M = 0 in Fig. 1) the output noise V (k) may also contain output measurement errors. 2.2. Gaussian ML cost function Under Assumption 1, Y (k) is independent (over k), circular complex normally distributed. To construct the likelihood function (U (k) is known exactly) it is sufficient to calculate the mean and covariance of Y (k) given the model parameters and the covariance matrix of the driving white noise source E(k). Since the process noise V (k) is correlated with the input of the plant U (k), and since it is independent of the reference signal R(k), the expected values in the mean and covariance calculation of Y (k) should be conditioned on R(k). The latter is known since the controller is known and since Y (k) and U (k) are observed without errors. Using (1) we find that Y (k) = (Iny + G(k )M(k ))−1 G(k )R(k) + (Iny + G(k )M(k ))−1 H (k )E(k)
(3)
and E{Y (k)|R(k), , } = (Iny + G(k , )M0 (k ))−1 G(k , )R(k) ≡ Y (k, ), Cov(Y (k)|R(k), , ) = S(k , )S H (k , ) ≡ CY (k) (), (4) S(k , ) = (Iny + G(k , )M0 (k ))−1 H (k , ). Hence, within a -independent additive constant, the negative Gaussian log-likelihood function is given by VNLL (, , Z) = log det(CY (k) ()) + k∈K
k∈K
(Y (k) − Y (k, ))H CY−1 (k) ()(Y (k) − Y (k, )),
(5)
where Z is a vector containing the input/output data (see Pintelon & Schoukens, 2001, p. 437, Eq. (14-14)). At DC (f = 0) and Nyquist (f = fs /2) the sums in (5) are multiplied by 21 . Using S −1 (k , )(Y (k) − Y (k, )) = H −1 (k , )(Y (k) − G(k , )U (k)) ≡ (k , ), log det(CY (k) ()) = log | det(S(k , ))|2 + log det(). (6)
870
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
Eq. (5) can be written as
Discussion:
VNLL (, , Z) = log | det(S(k , ))|2 k∈K
+ F log det() +
• If the raw data are time domain signals, then the model Y (k) = G(k , )U (k) + H (k , )E(k)
H (k , )−1 (k , ).
(7)
k∈K
For the open loop case (M0 =0 ⇒ S =H ) Eq. (7) divided by F reduces to Eq. (9) in McKelvey and Ljung (1997). Minimizing the real part of (7) w.r.t. (VNLL is a real number) gives 1 H (k , ) (k , ) , (8) () = Re F k∈K
k∈K
1 + F log det Re (k , )H (k , ) F k∈K
+ F ny .
(9)
For the open loop case (S = H ) Eq. (9) divided by F reduces to Eq. (12) in McKelvey and Ljung (1997), except that in the latter () is not constrained to be real and, hence, is sub-optimal. Dividing (9) by F and taking the exponential function gives, within a -independent multiplicative constant, 1 2 H (k , ) (k , ) , VF (, Z) = |hF ()| det Re F hF () = exp
is asymptotically (number of time domain samples N → ∞) valid. To improve the finite sample behaviour of the estimates the plant and noise transient terms are added in (11) (k , ) = H −1 (k , )(Y (k) − G(k , )U (k) − TG (k , ) − TH (k , )),
(13)
where, for example,
where Re(x) stands for the real part of x. Combining (7) and (8) gives log | det(S(k , ))|2 VNLL (, Z) =
(12)
k∈K
1 log det(S(k , )) . F
G = A−1 B,
TG = A−1 IG ,
H = D −1 C,
TH = D −1 IH
(14)
with A, B, C, D, IG and IH matrix polynomials of the appropriate matrix dimensions (see Pintelon & Schoukens, 2001). The coefficients of IG and IH decrease as an O(N −1/2 ) and, hence, for N sufficiently large the transients terms can be neglected in (13). • If the following conditions are satisfied: (i) k = zk−1 with zk = exp(j2k/N ) and k ∈ K = {0, 1, . . . , N − 1}, (ii) det(S(z−1 , )) has all its roots within the unit circle, and (iii) lim det(S(z−1 , )) = 1,
(15)
z→∞
then (11) simplifies to VN (, Z) = det
(10)
N−1 1 −1 (zk , )H (zk−1 , ) N k=0
k∈K
+ O(|max |N /N )
Taking the ny th root of |hF the factor can be brought into the determinant giving 1 H (k , ) (k , ) , VF (, Z) = det Re F ()|2 ,
k∈K
(k , ) = gF ()(k , ), 1 log det(S(k , )) , gF () = exp F ny
(11)
k∈K
where S(k , ) and (k , ) are defined in (4) and (6). This results in the following theorem:
(16)
with max the dominant pole of log det(S(z−1 , )) (proof: first note that the sum over the full unit circle is real, and next follow the lines of Appendix A.3 in Pintelon & Schoukens, 2005). Using (13), the first term in the right-hand side of (16) equals exactly the classical time domain result (see Ljung, 1999, p. 219, Eq. (7.96)) det
N−1 1 −1 (zk , )H (zk−1 , ) N k=0
= det
N−1 1 (t, )T (t, ) N
(17)
t=0
Theorem 1 (ML cost function). Under Assumption 1, the Gaussian ML cost function is given by (11), where at DC (f = 0) and Nyquist (f = fS /2) the summands are multiplied by 21 .
with (t, )=H −1 (q, )(y(t)−G(q, )u(t)) and q the backward shift operator (qx(t) = x(t − 1)).
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
2.3. ML estimator Consider, for example, a common denominator representation of the plant and noise models nb Br r B(, ) G(, ) = = r=0 na r , A(, ) r=0 ar nig Igr r IG (, ) TG (, ) = = r=0 na r , A(, ) r=0 ar nc Cr r C(, ) H (, ) = = r=0 nd r , D(, ) r=0 dr nih Ihr r IH (, ) = r=0 TH (, ) = nd r D(, ) r=0 dr
(18)
with A and D polynomials, IG and IH vector polynomials of size ny × 1, and B and C, matrix polynomials of size ny × nu and ny × ny , respectively. The parametric models (18), where contains all the numerator and denominator polynomials coefficients, are overparametrized. Indeed, multiplying the numerators and denominators in (18) by a non-zero real number leaves the parametric models unchanged. This imposes a constraint on the coefficients of the A and D polynomials. Further, the fact that E(k) in (3) is not observed, and that H (k , )E(k) = H (k , )T T −1 E(k) for any regular ny ×ny matrix T , imposes n2y additional constraints on the coefficients of the matrix polynomial C. The following parameter constraints are often used: z-domain
s-domain
ARMAX, BJ, OE :
a0 = 1
or ana = 1
ARMA, ARMAX, BJ :
C0 = Iny
or Cnc = Iny
BJ :
d0 = 1
or dnd = 1,
(19)
where for CT models the constraints are imposed on the normalized model parameters; other choices are, however, possible (see Section 3.1). Theorem 2 (ML estimator). Under Assumption 1 the ML ˆ estimator (Z) of model structure (12), where G(, ) and H (, ) are, for example, parametrized as in (18), and with a vector containing the plant and noise transfer function coefficients, minimizes (11) subject to the constraints (19). Discussion • It can easily be verified that the cost function VF (, Z) (11) contains exactly the same parameter ambiguities as the plant and noise models in (18). Using the results of Pintelon, Schoukens, Vandersteen, and Rolain (1999), it follows that ˆ and H (, ) ˆ 1/2 (), ˆ with the estimated models G(, ) ˆ the minimizer of (11), are independent of the particular parameter constraint(s) chosen. • Using H −1 (, ) = D(, ) adj (C(, ))/ det(C(, )), with adj(C) the adjoint matrix of C, it can be seen that the cost function VF (, Z) (11) depends on D(k , ) and
871
det(C(, )) via |D(k , )|2 and det(C(k , ))|2 only. It shows that the cost function remains the same when the roots of D(, ) that are not in common with A(, ) and the roots of det(C(, )) are mirrored w.r.t. the stability border ( → −, ¯ where x¯ denotes the complex conjugate of x, for continuous-time systems, and → −1 for discrete-time systems). This is not true for the zeroes of the polynomial matrix adj(C(, )), unless the MIMO problem reduces to ny independent SISO problems. Note that neither H nor H −1 should be stable for the calculations in the frequency domain. For time domain calculations H (D) and H −1 (det(C)) can be stabilized without modifying the statistical properties of the estimate. The stabilization of H −1 , however, requires symbolic calculations for multiple output systems. This is the bottle neck of the time domain approach. It can only be avoided via a specially parametrized state space model that excludes common dynamics in the plant and noise models. • Under some additional standard assumptions on the true model, the model set and the cost function (see Pintelon & Schoukens, 2005, Assumptions 9–11) it can be shown that the ML estimator in Theorem 2 is strongly consistent, asymptotically normally distributed, and asymptotically efficient (see Caines, 1988). • Following exactly the same lines of Pintelon and Schoukens (2005) the robustness of the ML estimator ˆ in Theorem 2 can be shown. The estimate (Z) remains strongly consistent and asymptotically normally distributed for non-Gaussian, mixing frequency domain noise V (k) (the mixing condition allows for a correlation over the frequency). Moreover, using the definition 1 VF () = det Re E (k , ) H (k , ) F k∈K
(20) the general expressions for the convergence rate and the ˆ asymptotic covariance matrix of (Z) remain valid (see Pintelon & Schoukens, 2005, Theorem 3, Properties 1–5). 3. Minimization of the ML cost function 3.1. Newton–Gauss algorithm The minimizer of (11) ()) with VF (, Z) = det( F 1 H () = Re (k , ) (k , ) F
(21)
k=1
is found as the solution of jVF (, Z)/ji = 0, i = 1, 2, . . . , dim(). Since j log det(X)/jX = X −T it follows that
jX j det(X) , (22) = det(X) trace X −1 ji ji
872
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
where X = (). Hence, jVF (, Z) =0 ji
⇒
j () −1 =0 trace () ji
(23)
(JL ())re = −(L ())re ,
for i = 1, 2, . . . , dim(). Using
j () 1 j (k , ) H = Re 2 herm , (k , ) ji F ji k
−T /2 () −1/2 (), −1 () = Y herm (X)Y H = herm (Y XY H ), trace(herm()) = Re(trace()),
(24)
with herm(X) = (X + X H )/2, Eq. (23) can be written as j () −1 trace () ji 2 j (k , ) = Re trace −1/2 () F ji k
na
ar2
+
r=0
ARMA, ARMAX, BJ :
nb
trace(BrT Br ) = 1
r=0 nc
CrT Cr = Iny
r=0
BJ :
, ) j ( 2 k ( −1/2 () (k , ))H = Re −1/2 () F ji
nd
dr2 = 1
(31)
r=0
because they are numerically more robust than (19) (constraints (19) result in an ill conditioned set of Eqs. (30) if, for example, in the z-domain the true value of a0 is zero).
k
2 j () H = Re (L() , ()) L() F ji
(30)
where ( )re puts the real and imaginary parts on top of each other, via a QR-factorization or an SVD decomposition (Golub & Van Loan, 1996). The iterative Newton–Gauss algorithm is implemented as follows. First we calculate (JL ())re w.r.t. all model parameters (constraints (19) are not imposed), next we solve (30) using the pseudo-inverse of (JL ())re (the dimension of the null space of (JL ())re equals one for OE, two for ARMA(X), and three for BJ), and finally we impose the following constraints on the updated model parameters + : ARMAX, BJ, OE :
× ( −1/2 () (k , ))H
(), and J()=j ()/j, with JL ()=L()J(), L ()=L() and where (), L() and () are defined in (21) and (26), respectively. A numerical stable implementation consists in solving
(25) 3.2. Starting values
where −1/2
L() = IF ⊗
()
and
() = [ T (1 , ), . . . , T (k , ), . . . , T (F , )]T
(26)
with ⊗ the Kronecker product. From (25) we conclude that the minimizer of (21) is found as the solution of the following set of equations:
j () Re (L() ())H L() = 0. (27) j The derivative of the left-hand side of (27) w.r.t. equals
j j () H Re (L() ()) L() j j
j () H j () = Re L() L() + O( ). (28) j j Neglecting the O( ) terms in (28), the Newton–Raphson algorithm for finding the solution of (27) reduces to the Newton–Gauss method (Fletcher, 1991) Re(JLH ()JL ()) = −Re(JLH ()L ())
(29)
First starting values for the plant model are calculated by minimizing the following weighted non-linear least squares cost function: 1 H (k , )Cˆ Y −1 (k , ) subject to constraints(31) (32) F k∈K
w.r.t. , where Cˆ Y is a non-parametric estimate of the noise model (see Appendix A), and where (32) is started from a weighted linear least squares version of (32). Next starting values for the noise model are obtained in three steps: (i) calcuˆ (k), lation of the output residuals Vˆ (k) = Y (k) − G(k , )U with ˆ the minimizer of (32); (ii) estimate of the denominator polynomial D via a SISO-ARMA fit of ⎛
⎞1/2 ny |Vˆ[p] (k)|2 ⎝ ⎠ var(Vˆ[p] )
(33)
p=1
with var(Vˆ[p] ) = (1/F ) k∈K |Vˆ[p] (k)|2 ; (iii) estimate of the numerator polynomial matrix C via a MIMO-MA fit of the ˆ Vˆ (k), with D(k , ) ˆ the solution of filtered residuals D(k , ) step (ii).
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
4. Calculation of the Cramér–Rao lower bound The ultimate goal is to calculate uncertainty bounds on invariants of model structure (12) starting from the Cramér–Rao lower bound on the estimated model parameters. Examples of invariants are the plant transfer function G(, ) and its poles and zeroes, and the noise power spectrum H (, )H H (, ). Note that the uncertainty calculation of H H H only requires the uncertainty of H 1/2 and not that of H and separately. Hence, the calculation of the Cramér–Rao lower bound can be simplified significantly by imposing the constraint = Iny . This is done by right multiplication of the estimated matrix coˆ giving C ˆ Since H 1/2 r = Cˆ r 1/2 (). efficients Cˆ r by 1/2 () can still be right multiplied by an orthogonal matrix T without changing the noise power spectrum H H H , the symmetry of r should also be imposed. This one of the matrix coefficients C is done as follows: first one matrix coefficient is decomposed 0 = U0 0 V T , and next all in singular values, for example, C 0 r are right multiplied by T0 = V0 U T . Summatrix coefficients C 0 marized the Cramér–Rao lower bound of the model parameters subject to the constraints: z-domain
s-domain
ARMA, ARMAX, BJ, OE :
= Iny
or = Iny
ARMAX, BJ, OE :
a0 = 1
or ana = 1
ARMA, ARMAX, BJ :
C0 = C0T
or
BJ :
d0 = 1
or dnd = 1
(34)
Cnc = CnTc
Using constraints (34), the negative log-likelihood function is, within a -independent constant, given by (7) with = Iny . Applying the definition of the Cramér–Rao lower bound CR (0 ). CR−1 (0 ) = j2 E{VNLL (, Iny , Z)}/j20
CR−1 (0 ) = 2(J1 )Tre (J1 )re + (J2 )Tre (J2 )re ,
vec(WF )
1/2
Vk = S0−1 (k )S(k , )S H (k , )S0−H (k ), where vec( ) stacks the columns on top of each other, 1/2 X0 (k ) = X(k , 0 ), and SR (k) is a square root of H E{R(k)R (k)}. If the expected value in (35) is taken w.r.t. to 1/2 the disturbing noise E(k) only, then (36) with SR (k) = R(k) is the Cramér–Rao lower bound for the particular realization of the reference signal R(k). Finally, the Cramér–Rao lower bound on any model parameter related quantity f () is obtained via
20
100
-40 4
5
6
7
0 -100 -200
8
4
40
200
20
100 Phase (°)
Amplitude (dB)
Frequency (Hz)
0 -20 -40 4
6 Frequency (Hz)
8
vec(VF ),
Wk = H0 (k )G(k , )(Inu + M0 (k )G0 (k ))−1 SR (k),
200
-20
(36)
where ( )re stacks the real and imaginary parts on top of each other. J1 and J2 are given by ⎡ ⎡ ⎤ ⎤ vec(W1 ) vec(V1 ) ⎢ ⎢ ⎥ ⎥ vec(W2 ) ⎥ vec(V2 ) ⎥ j ⎢ j ⎢ ⎢ ⎢ ⎥ ⎥ and, J1 = ⎥ , J2 = j ⎢ ⎥ j0 ⎢ 0 ⎣ ... ⎦ ⎣ ... ⎦
40
0
(35)
we find after some calculations
Phase (°)
Amplitude (dB)
will be calculated.
873
5
6 7 Frequency (Hz)
8
6 Frequency (Hz)
8
0 -100 -200
4
Fig. 2. Flight flutter data: measured frequency response function (gray bullets), plant model (bold solid line), and standard deviation plant model (bold dashed line). Left column: amplitude, right column: phase.
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
20
20
0
0
Amplitude (dB)
Amplitude (dB)
874
-20 -40
-20 -40 -60
-60 4
5
6
7
8
4
5
4
5
20
20
0
0
Amplitude (dB)
Amplitude (dB)
Frequency (Hz)
-20 -40 -60
6 7 Frequency (Hz)
8
-20 -40 -60
4
5
6
7
8
Frequency (Hz)
6
7
8
Frequency (Hz)
Fig. 3. Flight flutter data: measured noise power spectrum (gray bullets), modelled noise power spectrum (bold solid line), and standard deviation modelled noise power spectrum (bold dashed line).
CR(f (0 )) = f (0 )CR(0 )f H (0 ), where denotes the derivative w.r.t. (see Caines, 1988; Pintelon & Schoukens, 2001). 5. Measurement example: flight flutter analysis Although the asymptotic properties (consistency, asymptotic covariance matrix, robustness, etc.) of the ML estimator has been verified on numerous simulation examples, we prefer to report the results of a real life application: flight flutter analysis. To excite the air plane during flight, a perturbation signal is injected in the control loop of the flap mechanism at the tip side of the right wing. The angle perturbation of the flap is used as a measure of the applied force, and the acceleration is measured at two different positions on the left wing. Beside the applied perturbation the air plane is also excited during flight by the natural turbulence. The resulting turbulent forces acting on the air plane cannot be measured and are assumed to be white in the frequency band of interest. Since the input/output signals are lowpass filtered before sampling, a continuous-time plant and noise model is the natural choice. The signals are measured during about 128 s at the sampling rate fs = 256 Hz, giving N = 32 768 data points per channel. Fig. 2 shows the 2 × 1 force to acceleration frequency response function (FRF) in the band [4.00, 7.99 Hz] (DTF lines k = 513, 514, . . . , 1023 ⇒ F = 511). From Figs. 2 and 3 it can be seen that a continuoustime ARMAX model ((18) with = s and D = A) of order na = 8, nb = 8, nig = 9, and nc = 8 explains the measurements
very well. The non-parametric FRF and noise power spectra shown in these figures (gray bullets) are calculated as ˆ ˆ ˆ k ) = Y (k) − TG (sk , ) − TH (sk , ) , G(s U (k) ˆ (k) − TG (sk , ) ˆ Vˆ (k)Vˆ H (k) with Vˆ (k) = Y (k) − G(sk , )U ˆ − TH (sk , ),
(37)
where ˆ is the ML estimate. 6. Conclusion It has been shown that all results for SISO BJ identification on parts of the imaginary axis (s-domain) or unit circle (zdomain) gracefully extend to the multivariable case. Although the computational aspects of the ML estimator (Newton–Gauss minimization, parameter constraints, starting values, Cramér–Rao lower bound) have been discussed for one particular parametrization of the transfer function models (common denominator), exactly the same lines can be followed for other parametrization (for example, left matrix fraction description). Acknowledgement This work is sponsored by the Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (GOA-ILiNoS)
R. Pintelon et al. / Automatica 43 (2007) 868 – 875
and the Belgian Program on Interuniversity Poles of Attraction initiated by the Belgian State, Prime Minister’s Office, Science Policy programming (IUAP V/22). Appendix A. Non-parametric noise model The data records of N input/output samples are split in M blocks of equal length N/M , and the input/output DFT spectra U [r] (k), Y [r] (k) of each block are calculated (r =1, 2, . . . , N/M ). To reduce the leakage effect the spectra are differenced (Schoukens, Rolain, & Pintelon, 2006) YD[r] (k + 0.5) = Y [r] (k + 1) − Y [r] (k), UD[r] (k + 0.5) = U [r] (k + 1) − U [r] (k)
(38)
and M is kept as small as possible. Next, the noise covariance matrix is estimated as
ˆH Cˆ y (k) = (SˆYD YD (k) − SˆYD UD (k)Sˆ −1 (39) UD UD (k)S YD UD ), 2 where =M/(M −nu ) is a bias correction factor, 21 accounts for [r] [r] H the differencing (38), and where SˆXY = M −1 M r=1 X (Y ) . From (39) it follows that the minimal value of M equals nu +1. Finally, in order to guarantee the existence of the second order ˆ moments of Cˆ −1 Y (k) for a given input, CY (k) is averaged over Q neighbouring frequencies such that Q = (ny + 3)/(nu + 1) (Mahata, Pintelon, & Schoukens, 2004). To reduce the variability of the estimate we mostly choose Q > (ny + 3)/(nu + 1). References Brillinger, D. R. (1981). Time series: data analysis and theory. New York: McGraw-Hill. Caines, P. E. (1988). Linear stochastic systems. New York: Wiley. Fletcher, R. (1991). Practical methods of optimization. 2nd ed., New York: Wiley. Golub, G. H., & Van Loan, C. F. (1996). Matrix computations. 3rd ed., Baltimore: John Hopkins University Press. Ljung, L., (1993). Some results on identifying linear systems using frequency domain data. In Proceedings 32nd IEEE conference on decision and control, (Vol. 4, pp. 3534–3538) San Antonio, Texas, USA, December 15–17. Ljung, L. (1999). System identification: theory for the users. Upper Saddle River (USA): Prentice-Hall. Mahata, K., Pintelon, R., & Schoukens, J., (2004). Frequency domain identification using non-parametric noise models. In Proceedings of the 43rd IEEE conference on decision and control, (Vol. 1, pp. 821–826) Atlantis, Bahamas, December 14–17. McKelvey, T. (2002). Frequency domain identification methods. Circuits and Systems Signal Processing, 21(1), 39–55.
875
McKelvey, T., & Ljung, L., (1997). Frequency domain maximum likelihood identification. In 11th IFAC symposium on system identification, (Vol. 4, pp. 1741–1746) Kitakyushu, Japan, 8–11, July, preprints. Pintelon, R., Rolain, Y., & Schoukens, J. (2005). Box–Jenkins identification revisited—Part II: Applications. Automatica, 42(1), 77–84. Pintelon, R., & Schoukens, J. (2001). System identification: A frequency domain approach. New York (USA): IEEE Press. Pintelon, R., & Schoukens, J. (2005). Box–Jenkins identification revisited—Part I: Theory. Automatica, 42(1), 63–75. Pintelon, R., Schoukens, J., & Guillaume, P. (2006). Continuous-time noise modelling from sampled data. IEEE Transactions on Instrumentation and Measurement, 55(6), 2253–2258. Pintelon, R., Schoukens, J., Vandersteen, G., & Rolain, Y. (1999). Identification of invariants of (over)parametrized models: Finite sample results. IEEE Transactions on Automatic Control, AC-44(5), 1073–1077. Schoukens, J., Rolain, Y., & Pintelon, R. (2006). Leakage reduction in frequency response function measurements. IEEE Transactions on Instrumentation and Measurement, 55(6), 2286–2291. Rik Pintelon was born in Gent, Belgium, on December 4, 1959. He received the degree of electrical engineer (burgerlijk ingenieur) in July 1982, the degree of doctor in applied sciences in January 1988, and the qualification to teach at university level (geaggregeerde voor het hoger onderwijs) in April 1994, all from the Vrije Universiteit Brussel (VUB), Brussels, Belgium. From October 1982 till September 2000 he was a researcher of the Fund for Scientific Research—Flanders at the VUB. Since October 2000 he is professor at the VUB in the Electrical Measurement Department (ELEC). His main research interests are in the field of parameter estimation/system identification, and signal processing. Johan Schoukens is born in Belgium in 1957. He received the degree of engineer in 1980 and the degree of doctor in applied sciences in 1985, both from the Vrije Universiteit Brussel, Brussels, Belgium. He is presently professor at the Vrije Universiteit Brussel. The prime factors of his interest are in the field of system identification for linear and nonlinear systems, and growing tomatoes in his green house.
Patrick Guillaume was born in Anderlecht, Belgium, in 1963. He received the degree in mechanical-electrotechnical engineering in 1987 from the Vrije Universiteit Brussel, Belgium. In 1987 he joined the Department of Fundamental Electricity and Intrumentation (ELEC) of the Vrije Universiteit Brussel and received the Ph.D. degree in 1992. In 1996 he joined the Department of Mechanical Engineering (MECH) where he is currently lecturer. His main research interests are situated in the field of system identification, signal processing and experimental modal analysis.