- Email: [email protected]
Identification of Multivariable Canonical State-Space Representations
Proceedings of the 15th IFAC Symposium on System Identification Saint-Malo, France, July 6-8, 2009
Identification of Multivariable Canonical State-Space Representations L. Bako ∗ G. Merc` ere ∗∗ and S. Lecœuche ∗ ∗ Ecole des Mines de Douai, D´epartement Informatique et Automatique, 941, rue Charles Bourseul, 59508 Douai Cedex, France (e-mail: [email protected], [email protected]) ∗∗ Universit´e de Poitiers, Laboratoire d’Automatique et d’Informatique Industrielle, 40 avenue du recteur Pineau, 86022 Poitiers, France (e-mail: [email protected])
Abstract: The paper introduces a method to estimate a linear state-space model from its input-output data. Its starting point is the conversion of the state-space model into an inputoutput ARMAX-like model. In order to conveniently estimate the parameters of this ARMAX model, it is necessary, given the structural characteristics of the initial state-space model, to assume that the dynamic state matrix A is nonderogatory. Under this condition, it is shown that the input-output model and the state-space model from which it is derived have the same order. Based on this important result in linear system realization theory, a canonical statespace representation is finally obtained through a simple construction technique. The developed method can be regarded as an alternative to the standard subspace methods particularly in situations where the state coordinates basis need to be maintained constant with respect to the input-output realization used for the estimation. Keywords: Least-squares algorithm ; State-space representation ; Canonical form ; ARMAX model ; Multi-input/multi-output systems. 1. INTRODUCTION
During the last two decades, interesting developments have been achieved in linear system identification thanks the so-called subspace-based methods (Katayama, 2005). One of the main advantages of this approach is that it gives directly access to minimal state-space representation without requiring the use of canonical parameterizations and the determination of structural indices demanded by classic prediction error methods when multivariable systems are studied. Although the subspace-based methods have proved their efficiency from a practical point of view (Mevel et al., 1999), (Favoreel et al., 2000), (Lovera, 2003), most of the techniques developed until now suffer from two main drawbacks which can be crucial in particular cases, e.g, when hybrid systems or on-line identification framework are involved. More precisely, • the user of the common subspace algorithms has difficulties to fix the basis in which the state-space models are estimated. In fact, this choice is realized during the reduction step (when the system order is estimated) by using generally a Singular Value Decomposition (SVD) which gives rise to particular fully-parametrizations. But the reproducibility of the basis is not ensured when such a tool is employed. This characteristic can be problematic when hybrid systems such as switched or LPV processes are identified since the different sub-models must be expressed in coherent bases. Otherwise, the input-output behav-
978-3-902661-47-0/09/$20.00 © 2009 IFAC
ior of the global model may be distorted (T´oth et al., 2007). • the SVD, interesting for its computational robustness, is numerically cumbersome. Thus, this tool is difficult to use in a recursive estimation framework and alternatives must be developed (Lovera et al., 2000), (Merc`ere et al., 2008) to overcome this problem when it is necessary to update the model in real time. To sort out these two problems, a new SVD-free algorithm for the identification of linear MIMO state-space models is suggested in this paper. It is based on an initial estimation step of an ARMAX model using traditional least squares algorithms followed by an algebraic reconstruction of a canonical state-space representation from the parameters calculated during the first phase. This approach is inspired by the works of Stoica and Jansson (2000) where it is shown that canonical transfer function approach can be seen as an interesting alternative to the subspace methods even when MIMO systems are considered. The introduction of least squares optimization techniques is an important asset of the developed method since it avoids the use of SVD and make its adaptation easier for on-line scenarii. Furthermore, as shown afterwards, the estimated state-space matrices will be given in a fixed coordinate basis. Thus, this algorithm can be easily extended when composite systems are involved. The outline of the paper is as follows. Section 2 is concerned with the presentation of the model, the relative assumptions and a statement of the identification problem. In Section 3, the state-space model is transformed into
1638
10.3182/20090706-3-FR-2004.0192
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 an ARMAX-like model whose parameters are estimated under an appropriate assumption on the A-matrix. Given this latter model, an original but quite simple technique to construct a state-space model is suggested in Section 4. Conditions that guarantee equality of orders between the initial state-space model and its corresponding ARMAXlike model are also derived. 2. PROBLEM FORMULATION AND NOTATIONS Consider the following finite dimensional, causal and linear time-invariant state-space model x(t + 1) = Ax(t) + Bu(t) + w(t) (1a) y(t) = Cx(t) + Du(t) + v(t) (1b) where u(t) ∈ Rnu , y(t) ∈ Rny and x(t) ∈ Rn are respectively the input, the output and the state vectors. The process noise sequence is w(t) ∈ Rn and the output measurement disturbances are collected in v(t) ∈ Rny . (A, B, C, D) are the system matrices relatively to a certain coordinate basis of the state-space. It will be assumed that {x(t)}, {u(t)}, {y(t)}, {w(t)} and {v(t)}, t ∈ Z, are all ergodic and (weakly) stationary stochastic processes (Ljung, 1999), the system is asymptotically stable, controllable and observable and there is no feedback from y to u. The considered identification problem can then be stated N N as follows: given realizations {u(t)}t=1 and {y(t)}t=1 of the input and output processes generated by a system of the form (1) on a finite but sufficiently wide time horizon N , estimate the matrices (A, B, C, D) up to a similarity transformation. To reach this goal, two steps will be considered. Firstly, an input-output ARMAX model will be identified from the input-output measurements by using a classic least squares or instrumental variable method. Then, under the assumption that the matrix A is nonderogatory (Horn and Johnson, 1990), a canonical state-space form will be reconstructed from the estimated ARMAX model parameters. 3. ARMAX MODEL ESTIMATION Consider the output equation (1b). A direct way to get an input-output representation of the system from this statespace form is to get rid of the state. For that, let us start by expressing the state as a linear combination of the past inputs, outputs and a prior value of the state. Particularly, for t > n, it is easy to show that x(t) = An x(t − n) + ∆n un (t − n) + An wn (t − n) (2) with 1 ⊤ un (t − n) = u⊤ (t − n) · · · u⊤ (t − 1) ∈ Rnnu ∆n = An−1 B · · · AB B ∈ Rn×nnu 2 An = An−1 · · · A In ∈ Rn×n . Substituting Eq. (2) into Eq. (1b), we get y(t) = CAn x(t − n) + C∆n un (t − n) + CAn wn (t − n) + Du(t) + v(t). (3) Now, by using the Cayley-Hamilton theorem (Horn and Johnson, 1990), we know that An = −a0 In − a1 A − 1
yn (t − n) and wn (t − n) are built in a similar way.
· · · − an−1 An−1 where the {ai }i∈[0,n−1] elements are the coefficients of the characteristic polynomial of A, i.e., pA (z) = a0 + a1 z + · · · + an−1 z n−1 + z n . Thus, CAn x(t − n) = − a⊤ ⊗ Iny Γn x(t − n) (4) where (5) a⊤ = [a0 a1 · · · an−1 ] i⊤ h (6) Γn = C⊤ (AC)⊤ · · · CAn−1 ⊤
and with ⊗ the Kronecker product (Horn and Johnson, 1990). By using now the recurrent relations (1), it is straightforward to verify that the stacked vectors yn (t−n), un (t − n), vn (t − n) and wn (t − n) are related as follows yn (t − n) = Γn x(t − n) + Hn un (t − n) + Gn wn (t − n) + vn (t − n) (7) with 0 D 0 ··· 0 0 ··· 0 .. .. .. .. C . . . . . , G = Hn = CB n .. .. .. .. .. .. . .0 . . 0 . . CAn−2 B ··· CB D
CAn−2 ···
C 0
Eq. (7) allows us to express Γn x(t − n) as a combination of past inputs, outputs and disturbances, i.e., Γn x(t − n) = yn (t − n) − Hn un (t − n) − Gn wn (t − n) − vn (t − n). (8) Finally, gathering equations (3), (4) and (8), the state vector is eliminated in the output equation and we get y(t) = − a⊤ ⊗ Iny yn (t − n) + Fun+1 (t − n) + η(t) (9) with F = a⊤ ⊗ Iny Hn + C∆n D ∈ Rny ×(n+1)nu η(t) = a⊤ ⊗ Iny Gn + CAn wn (t − n) + a⊤ ⊗ Iny vn (t − n) + v(t). Note that F contain the coefficients of the polynomial matrix 2 that defines the numerator of the transfer function of (1) and a is the coefficient vector of the denominator. The model (9) is in fact quite more complex than a standard ARMAX since it contains at the same time a measurement and a process noise. It is also interesting to note that the noise part of this ARMAX model is composed of particular combinations of the AR and X deterministic parameters. Thus, as soon as a and F are consistently estimated, a model of the stochastic part can be built. In order to simplify the estimation of a and F, it is interesting to notice that Eq. (9) can be rewritten as a classic linear regression. For that, introduce the identity vec (AXB) = B⊤ ⊗ A vec (X) where vec (•) is the vectorization operator (Horn and Johnson, 1990). Then, a + η(t) y(t) = −Φyn (t − n) u⊤ (t − n) ⊗ I n y n+1 vec (F) (10) with Φyn (t − n) = [y(t − n) · · · y(t − 1)] ∈ Rny ×n . Thus, a and F can be estimated using a least squares algorithm such as the instrumental variable method (Ljung, 1999) when the data matrix 2 For an explicit expression of this polynomial matrix, we refer to the proof of Theorem 1 given in Appendix B.
1639
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
Φ(n) =
Φyn (1) .. .
u⊤ n+1 (1) ⊗ Iny .. .
Φyn (N − n) u⊤ n+1 (N − n) ⊗ Iny has full column rank.
requires that the matrix Φ(n) has full column rank. No matter how rich is the excitation input, it can be easily shown that rank(Φ(n)) ≤ nIO + (n + 1)nu ny in the noisefree case (i.e., when η(t) = 0 in Eq. (9)). It follows that it is necessary to satisfy nIO = n in order to have a full column rank matrix Φ(n).
4. FROM ARMAX MODEL TO MINIMAL CANONICAL STATE-SPACE REPRESENTATION 4.1 Condition to get ARMAX and state-space models of the same order The problem of parameterizing state-space models and the equivalence between state-space representation and MIMO canonical transfer functions is not new (see, e.g. (Kailath, 1980; McKelvey, 1995)). One of the main problems is the relation between the order nIO of the I/O model and the one of the corresponding minimal state-space representation. A glance at Eq. (9) and the way it is built from the state-space model (1) may tend to suggest that both representations have the same order. Unfortunately, this is not the case for all the MIMO systems (see (Stoica and Jansson, 2000) for a counterexample). Indeed, we have in general, nIO ≤ nSS ≤ nIO · ny , where nSS = n is the order of the state-space model (1). To see this fact more clearly, let us think of the order as the smallest integer m such that y(t) can be expressed in function of ym (t−m) and um+1 (t−m) as in Eq. (9). Then, intuitively, the order nIO of the input-output model (9) should be equal to the degree ν of the minimal polynomial of A. This remark follows essentially from Eq. (4) which suggests that nIO should be the smallest integer such that AnIO can be expressed in function of Aj , j = 0, . . . , nIO −1, where A is taken from a minimal realization of (1). On the other hand, if we consider the input-output transfer matrix of (1) to be TF(z) = L (z)/pA (z), with L (z) = Cadj(zI − A)B + pA (z)D f11 (z) · · · f1nu (z) .. = ... , . fny 1 (z) · · · fny nu (z) then the order nIO of Eq. (9) corresponds to o n j∈[1,n ] n − deg gcd {fij (z)}i∈[1,nyu] , pA (z) .
In fact, minimality of (A, B, C, D) is not enough to guarantee coprimeness between the polynomials fij (z) and pA (z). This is formally stated by the following theorem which relates the order of model (9) to that of (1). Theorem 1. Assume that the model (1) is minimal, that is, observable and reachable. Then the order nIO of the inputoutput model (9) is equal to the degree ν of the minimal polynomial of A, i.e., nIO = deg(mA (z)) = ν ≤ n. For the proof, see Appendix A. We know from Theorem 1 that, in general, the orders of the models (1) and (9) do not coincide, even when the model (1) is minimal. Therefore, even though it is assumed that the order n is known a priori, it is compulsory to determine nIO before being able to estimate a and F from (10). Indeed, consistent estimation of these parameters
To get equality between both orders, a special property on the matrix A in Eq. (1) is required as formulated by the following corollary of Theorem 1. Corollary 2. Assume that the model (1) is minimal. Then nIO = nSS = n if and only if the matrix A is nonderogatory. 3 In order to explicitly reconstruct a state-space representation from the estimated ARMAX parameters, it is essential to find conditions ensuring an equality relation between nIO and nSS . 4.2 Reconstruction of a minimal canonical state-space model Once the parameters of the model (9) are estimated, the reconstruction of a canonical state-space model can be considered. Before, the following lemma must be stated. Lemma 3. If the system (1) is observable, then there exists a matrix Ψ of the form ⊤ γ 0 ··· 0 .. 0 ... . ∈ Rn×nny , (11) Ψ= . . .. 0 .. 0 · · · 0 γ⊤ with γ ∈ Rny , and verifying rank (ΨΓn ) = n (12) where Γn is the observability matrix of order n, if and only if the matrix A is nonderogatory. A proof of this lemma can be found in (Bako et al., 2009). Note that the condition (12) in Lemma 3 can be interpreted as follows. There exists a γ such that if we let ya (t) = γ ⊤ y(t) be a linear combination of all the outputs, then the state x(t) of system (1) is entirely observable from ya (t). This further means that the system (1) and the MISO system with input u(t) and output ya (t) ∈ R have the same state. Lemma 3 tells us that this is in fact possible if and only if the matrix A is nonderogatory 4 . In view of that lemma, if, e.g, γ ∈ Rny is randomly generated from a uniform distribution, then the required condition (12) is fulfilled with probability one (Bako et al., 2009). By assuming now that the matrix A is nonderogatory (which ensures the equality nIO = nSS ), a minimal statespace realization can be readily obtained as explained in the following proposition: Proposition 4. Assume that the system (1) is minimal and that the matrix A is nonderogatory. Consider a matrix Ψ defined as in (11) and satisfying condition 3
A square matrix is said to be nonderogatory if its minimal polynomial coincides with its characteristic polynomial (Horn and Johnson, 1990). 4 In fact, for any SISO or MISO system to be observable the Amatrix needs to be nonderogatory.
1640
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 • an off-line version of the propagator method suggested initially in (Merc`ere et al., 2003) in a recursive form (see also (Merc`ere et al., 2008)) where a particular linear operator is introduced to rewrite the subspace-based identification problem into a quadratic optimization problem without constraint. For all these methods, the instrumental variable used to reduce the effect of the noise is built from the inputs of the system and the outputs of a simulated biased auxiliary model obtained with the ordinary MOESP algorithm (Verhaegen, 1994) applied on noisy data. Method in (Bako et al., 2008)
colnu ([ Y1
···
with Xi ∈ R
ny ו
and Yi ∈ R
, i ∈ [1, n].
As opposed to classic techniques leading to canonical forms for multivariable systems (Luenberger, 1967), (Mayne, 1972), (Kailath, 1980), the developed approach does not need any observability or controllability indices or an initial step leading to a non minimal realization estimate followed by a reduction procedure. Here is proposed a relatively straightforward method which gives access to minimal MIMO state-space representation from a standard vector difference equation. 5. SIMULATION EXAMPLE In order to illustrate the performances of the algorithm, the following system is considered # # " " C=
0.8736 0 0 −0.3881 , 0 0
−0.3814 0.6134 −0.3495 0.4044 −0.0624 −0.7160 , B = −0.5787 −0.5476 −0.1790
h
i
0.9397 0 1.1787 , 0 0 −1.3274
D=
h
−0.1
−0.5
0 0.5 Real part
1
i
0.5463 −0.5293 . 0 −2.4003
The input is a zero-mean white Gaussian noise with unit variance. The disturbances v and w are both zero-mean white Gaussian noise set such that the noise signal ratio equals 20dB in both cases. A Monte Carlo simulation with 100 realizations of noise and input is carried out. This new algorithm is compared with three other methods: • another SVD-free identification method developed in (Bako et al., 2008) where the subspace identification problem is recasted into a simple least squares problem by introducing a particular user-defined matrix which does not change the rank of the extended observability matrix when multiplying this latter matrix on the left hand side, • the traditional PO MOESP algorithm (Verhaegen, 1994),
0.1 0 −0.1
−1
−0.5
0 0.5 Real part
1
Method in (Merc` ere et al., 2008)
PO-MOESP
0.3 0.2
0.1 0 −0.1 −0.2 −1
Our algorithm
−0.2
0.2
For the proof, see appendix B.
A=
0
0.3
Yn
•×nu
0.2
0.1
Imaginary part
Xn ] ,
Xn
0.2
−1
"Y # 1 .. , Yn ]) = . Imaginary part
rowny
0.3
−0.2
n−1
1 0 ··· 0 " X #! 1 .. = [ X1 ··· .
0.3 Imaginary part
Imaginary part
(12). Then a minimal canonical state-space representation can be obtained from the estimated parameters a⊤ = [a0 · · · an−1 ] and F of the ARMAX model (9) as: 0 1 ··· 0 .. . . .. . 0 , . . (13) A= 0 0 ··· 1 −a0 −a1 · · · −an−1 D = F (:, nnu + 1 : (n + 1)nu ) , (14) B = ΨQ, (15) −1 ⊤ C = rowny (Q) Ω⊤ , (16) n Ωn Ω n where −1 colnu F(:, 1 : nnu ) − a⊤ ⊗ D , Q = K ⊗ Iny Ωn = B AB · · · An−1 B ∈ Rn×nnu , a1 · · · an−1 1 .. . 1 0 ∈ Rn×n , K= a
0.1 0 −0.1 −0.2
−0.5
0 0.5 Real part
1
−1
−0.5
0 0.5 Real part
1
Fig. 1. Estimated (o) and real poles (+) of the system. Figure 1 compares the position of the estimated poles of the models given by the four aforementioned algorithms with respect to the real system. A first conclusion is that all these methods are quite accurate since the estimated poles are centered on the true ones and the variances of the estimates are relatively small. We can however notice that, on this example, the new algorithm leads, without requiring the use of an SVD, to similar results as the PO MOESP algorithm. Therefore, the approach developed in this article seems to be an interesting alternative to the classic subspace-based techniques. 6. CONCLUSION In this paper, a simple SVD-free method has been proposed to get a state-space model, given a collection of its input-output measurements. In comparison with existing subspace identification techniques, the main motivation of this new method is to prevent the state representation basis from shifting according to the estimation data. This feature is highly desirable in, e.g., recursive identification of state-space models or in the estimation of composite models with linear submodels. From a realization theory point of view, this work highlights the behavior of certain types of linear MIMO systems that can indeed be handled as MISO systems in terms of state dynamics. Future research will focus on generalizing this procedure to more arbitrary models than the particular cases where the Amatrix is nonderogatory.
1641
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 j∈[1,n ] gcd {lij (z)}i∈[1,nyu] , pA (z)
Appendix A. PROOF OF THEOREM 1 We first state the following lemma. Lemma 5. Let A ∈ Rn×n and B = [b1 · · · bm ] ∈ Rn×m be some arbitrary matrices. Let mA (z) be the mini mal polynomial of A and define An−1 B · · · AB B = ∆n (A, B) ∈ Rn×nm , A b1 .. .. nm×nm ¯ ¯ A= , b = . ∈ Rnm . . ∈R A
bm
Then the following holds. If rank ∆n (A, B) = n, then ¯ = ν, where ν = deg(mA (z)). ¯ b) rank ∆n (A, Proof. Assume that rank ∆n (A, B) = n. We then need ¯ b) ¯ is equal to to show that the rank of the matrix ∆n (A, the degree ν of the minimal polynomial of A. We proceed by contradiction. Assume that there exists j < ν such ¯ , ¯ is a linear combination of the vectors A ¯ ib ¯ jb that A i = 0, . . . , j − 1. Then, one can find a set of scalars ¯+ ¯ = αj−1 A ¯ j−1 b ¯ jb α0 , · · · , αj−1 not all zero satisfying A ¯ · · · + α0 b. By Denoting with p the polynomial of degree j defined by p(z) = z j − αj−1 z j−1 − · · · − α0 , we get ¯ = 0. This is equivalent to p(A)bi = 0, i = 1, . . . , m. ¯ b p(A) As a consequence, p(A)B = 0 and p(A)∆n (A, B) = 0 because, as a polynomial of A, p(A) commutes with any power of A. Since ∆n (A, B) has full row rank, it follows that p(A) = 0. Therefore, if we call mA (z) the minimal polynomial of A, then, by definition of mA (z) we get that mA (z) divides p(z). This contradicts the fact that deg {p(z)} = j < deg {mA (z)} = ν. In conclusion, ¯ b) ¯ = ν. rank ∆n (A,
From this lemma, it can be concluded that under control¯ = n if ¯ b) lability assumption of system (1), rank ∆n (A, and only if the A-matrix is nonderogatory. For the proof of Theorem 1, we will need also the following Theorem given in (Barnett, 1973): Theorem 6. Let bi (z) = bi1 z n−1 + · · · bin , i = 1, · · · , p, pA (z) = a0 + a1 z + · · · + an−1 z n−1 + z n . Then the degree of the gcd of pA (z), b1 (z), · · · , bp (z) is equal to n − rank (Γn (Ac , M)) with 0 1 ··· 0 b1n · · · b11 .. . . . .. , A = . . .. . M = ... c . 0 0 ··· 1 bpn · · · bp1 −a0 −a1 · · · −an−1 With this material, Theorem 1 can be proved as follows:
Proof. The input-output transfert matrix is given by Cadj(zI − A)B + pA (z)D C(zI − A)−1 B + D = , pA (z) where adj(·) is the cofactor matrix of (·) and pA (z) = det(zI − A) is the characteristic polynomial of A. Let us introduce the following notation j∈[1,n ]
Cadj(zI − A)B = [lij (z)]i∈[1,nyu] , where lij (z), the numerator polynomial that corresponds to the transfert between the i-th output and the j-th input. To prove Theorem 1, we need to show that the degree of
is equal to n − ν, with ν = deg(mA (z)). To proceed, we start by pointing out that, as far as we are concerned with seeking the common factors between the polynomials lij (z) and pA (z), we can, without loss of generality, assume that D = 0. Now we compute the coefficients of the polynomials lij (z). The expressions of these coefficients can be directly obtained from the term F of Eq. (9) (i.e., the matrix that multiplies un (t−n)). At this step, it is important to notice that F = (a⊤ ⊗ Iny )Hn + C∆n = C∆n (L ⊗ Inu ), with 1 0 · · · 0 .. . an−1 1 L= . ∈ Rn×n . . ... ... 0 . a1 · · · an−1 1 Let bj be j-th column of B and c⊤ i be the i-th row of C. The n coefficients of the n − 1 degree polynomial lij (z) are given by the following row vector n−1 lij = c⊤ bj · · · Abj bj L. i A | {z } ∆n (A,bj )
Let
and then
l1j Nj = ... = C∆n (A, bj )L lny j
C∆n (A, b1 ) N1 .. nn ×n N = ... = L ∈ R y . C∆n (A, bnu ) Nnu as a concatenation of all the coefficient vectors of the polynomial matrix Cadj(zI− A)B. According to Theorem 6, we just need to show that rank (P) = ν, where N NAc P= .. .
NAcn−1 and Ac a companion matrix (see Eq. (13)). Notice that the matrix P is equal to C∆n (A,b1 )L C∆n (A,b1 ) C∆n (A,b1 )AL C∆n (A,b1 )LAc .. .. . . n−1 C∆n (A,b1 )LAn−1 C∆n (A,b1 )AL c .. .. = P1 = L, . . C∆n (A,bnu )L C∆n (A,bnu ) C∆n (A,bn )LAc C∆n (A,bn )AL u u .. .. . . n−1 n−1 C∆n (A,bnu )LAc
C∆n (A,bnu )AL
up to some rows permutations, where −an−1 1 .. .. . . AL = LAc L−1 = −a 0
1642
1
−a0
··· .. . ··· 0 ···
0 .. . . 1
0
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 For any k ≥ 0 and any bj , we have ∆n (A, bj )AkL = Ak ∆n (A, bj ). By making use of this fact, it follows that Γn .. ¯ ¯ b)L, P1 = ∆n (A, . Γn
with Γn = Γn (A, C). Since the system is observable, the matrix that contains Γn in the right hand side of the previous equation, has full column rank. Therefore, Lemma 5 allows us to conclude that rank P = rank P1 = ¯ = ν since L is nonsingular. ¯ b) rank ∆n (A, Appendix B. PROOF OF PROPOSITION 4
The reconstruction of a canonical form from the parameters a and F is in fact only based on shrewd calculations. For D, it is easy to see that D = F (:, nnu + 1 : (n + 1)nu ) . Knowing D, it is interesting to notice that F (:, 1 : nnu ) − a⊤ ⊗ D = a⊤ ⊗ Iny Hn + C∆n − a⊤ ⊗ D ⊤ = k⊤ 1 ⊗ Iny Γn B · · · kn ⊗ Iny Γn B
where k⊤ i is the ith row of K. Thus, by using the operator colnu defined in Proposition 4 and noticing that K ⊗ Iny is invertible, we get −1 Γn B = K ⊗ Iny colnu F (:, 1 : nnu ) − a⊤ ⊗ D . Introducing matrix Ψ, it holds that ¯ ΨΓn B = TB = B −1 = Ψ K ⊗ Iny colnu F (:, 1 : nnu ) − a⊤ ⊗ D since ΨΓn = T ∈ Rn×n can be viewed as a similarity transformation if rank (ΨΓn ) = n. The next step consists in reconstructing the state matrix A. This can be done by ¯ = TAT−1 knowing T, i.e. calculating A ¯ = ΨΓn A (ΨΓn )−1 A ⊤ γ CA CA −1 −1 .. = Ψ ... (ΨΓn ) = (ΨΓn ) . . CAn γ ⊤ CAn Now, using again the Cayley-Hamilton theorem, we have 0 ⊤ 1 ··· 0 ⊤ γ C γ CA .. .. . . .. .. . . . = . . 0 0 · · · 1 ⊤ n−1 ⊤ n γ CA γ CA −a0 −a1 · · · −an−1 = Ac ΨΓn . where Ac is a companion matrix (see Eq. (13)). Hence, ¯ = Ac ΨΓn (ΨΓn )−1 = Ac . A ¯ B ¯ and Γn B, it is easy to see that Finally, knowing A, ¯ =C ¯ B ¯ nB ¯ n−1 B ¯ ¯ ··· A rowny (Γn B) = rowny Γ ¯ which leads to the extraction of C.
REFERENCES
L. Bako, G. Merc`ere, and S. Lecœuche. On-line structured subspace identification with application to switched linear systems. International Journal of Control, 2009. Accepted for publication. S. Barnett. Matrices, polynomials and linear time invariant systems. IEEE Transactions of Automatic Control, 18:1–10, 1973. W. Favoreel, B. De Moor, and P. Van Overschee. Subspace state space system identification for industrial processes. Journal of Process Control, 10:149–155, 2000. R. Horn and C. Johnson. Matrix analysis. Cambridge University Press, 1990. T. Kailath. Linear Systems. Prentice Hall, Engelwood Cliffs, 1980. T. Katayama. Subspace methods for system identification. Springer, 2005. L. Ljung. System identification. Theory for the user. Prentice Hall Information and System Sciences Series, Upper Saddle River, 2nd edition, 1999. M. Lovera. Identification of MIMO state space models for helicopters dynamics. In Proceedings of the 13th IFAC Symposium on System Identification, Rotterdam, The Netherlands, August 2003. M. Lovera, T. Gustafsson, and M. Verhaegen. Recursive subspace identification of linear and non linear Wiener state space models. Automatica, 36:1639–1650, 2000. D. Luenberger. Canonical forms for linear multivariable systems. IEEE Transactions on Automatic Control, 12: 290–293, 1967. D. Mayne. A canonical model for identification of multivariable linear systems. IEEE Transactions on Automatic Control, 17:728–729, 1972. T. McKelvey. Identification of state space models from time and frequency data. PhD thesis, Link¨oping University, Link¨oping, Sweden, 1995. G. Merc`ere, S. Lecœuche, and C. Vasseur. A new recursive method for subspace identification of noisy systems: EIVPM. In Proceedings of the 13th IFAC Symposium on System Identification, Rotterdam, The Netherlands, August 2003. G. Merc`ere, L. Bako, and S. Lecœuche. Propagator-based methods for recursive subspace model identification. Signal Processing, 88:468–491, 2008. L. Mevel, L. Hermans, and H. Van der Auweraer. On the application of subspace-based fault detection methods to industrial structures. Mechanical Systems and Signal Processing, 13:823–838, 1999. P. Stoica and M. Jansson. MIMO system identification: state-space and subspace approximations versus transfer function and instrumental variables. IEEE Transactions on Signal Processing, 48:3087–3099, 2000. R. T´oth, F. Felici, P. Heuberger, and P. Van den Hof. Discrete-time LPV I/O and state-space representations. differences of behavior and pitfalls of interpolation. In Proceeding of the European Control Conference, Kos, Greece, July 2007. M. Verhaegen. Identification of the deterministic part of MIMO state space models given in innovations form from input output data. Automatica, 30:61–74, 1994.
L. Bako, G. Merc`ere, and S. Lecœuche. A least squares approach to the subspace identification problem. In In the 47th IEEE Conference on Decision and Control, Cancun, Mexico, December 2008.
1643
Identification of Multivariable Canonical State-Space Representations L. Bako ∗ G. Merc` ere ∗∗ and S. Lecœuche ∗ ∗ Ecole des Mines de Douai, D´epartement Informatique et Automatique, 941, rue Charles Bourseul, 59508 Douai Cedex, France (e-mail: [email protected], [email protected]) ∗∗ Universit´e de Poitiers, Laboratoire d’Automatique et d’Informatique Industrielle, 40 avenue du recteur Pineau, 86022 Poitiers, France (e-mail: [email protected])
Abstract: The paper introduces a method to estimate a linear state-space model from its input-output data. Its starting point is the conversion of the state-space model into an inputoutput ARMAX-like model. In order to conveniently estimate the parameters of this ARMAX model, it is necessary, given the structural characteristics of the initial state-space model, to assume that the dynamic state matrix A is nonderogatory. Under this condition, it is shown that the input-output model and the state-space model from which it is derived have the same order. Based on this important result in linear system realization theory, a canonical statespace representation is finally obtained through a simple construction technique. The developed method can be regarded as an alternative to the standard subspace methods particularly in situations where the state coordinates basis need to be maintained constant with respect to the input-output realization used for the estimation. Keywords: Least-squares algorithm ; State-space representation ; Canonical form ; ARMAX model ; Multi-input/multi-output systems. 1. INTRODUCTION
During the last two decades, interesting developments have been achieved in linear system identification thanks the so-called subspace-based methods (Katayama, 2005). One of the main advantages of this approach is that it gives directly access to minimal state-space representation without requiring the use of canonical parameterizations and the determination of structural indices demanded by classic prediction error methods when multivariable systems are studied. Although the subspace-based methods have proved their efficiency from a practical point of view (Mevel et al., 1999), (Favoreel et al., 2000), (Lovera, 2003), most of the techniques developed until now suffer from two main drawbacks which can be crucial in particular cases, e.g, when hybrid systems or on-line identification framework are involved. More precisely, • the user of the common subspace algorithms has difficulties to fix the basis in which the state-space models are estimated. In fact, this choice is realized during the reduction step (when the system order is estimated) by using generally a Singular Value Decomposition (SVD) which gives rise to particular fully-parametrizations. But the reproducibility of the basis is not ensured when such a tool is employed. This characteristic can be problematic when hybrid systems such as switched or LPV processes are identified since the different sub-models must be expressed in coherent bases. Otherwise, the input-output behav-
978-3-902661-47-0/09/$20.00 © 2009 IFAC
ior of the global model may be distorted (T´oth et al., 2007). • the SVD, interesting for its computational robustness, is numerically cumbersome. Thus, this tool is difficult to use in a recursive estimation framework and alternatives must be developed (Lovera et al., 2000), (Merc`ere et al., 2008) to overcome this problem when it is necessary to update the model in real time. To sort out these two problems, a new SVD-free algorithm for the identification of linear MIMO state-space models is suggested in this paper. It is based on an initial estimation step of an ARMAX model using traditional least squares algorithms followed by an algebraic reconstruction of a canonical state-space representation from the parameters calculated during the first phase. This approach is inspired by the works of Stoica and Jansson (2000) where it is shown that canonical transfer function approach can be seen as an interesting alternative to the subspace methods even when MIMO systems are considered. The introduction of least squares optimization techniques is an important asset of the developed method since it avoids the use of SVD and make its adaptation easier for on-line scenarii. Furthermore, as shown afterwards, the estimated state-space matrices will be given in a fixed coordinate basis. Thus, this algorithm can be easily extended when composite systems are involved. The outline of the paper is as follows. Section 2 is concerned with the presentation of the model, the relative assumptions and a statement of the identification problem. In Section 3, the state-space model is transformed into
1638
10.3182/20090706-3-FR-2004.0192
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 an ARMAX-like model whose parameters are estimated under an appropriate assumption on the A-matrix. Given this latter model, an original but quite simple technique to construct a state-space model is suggested in Section 4. Conditions that guarantee equality of orders between the initial state-space model and its corresponding ARMAXlike model are also derived. 2. PROBLEM FORMULATION AND NOTATIONS Consider the following finite dimensional, causal and linear time-invariant state-space model x(t + 1) = Ax(t) + Bu(t) + w(t) (1a) y(t) = Cx(t) + Du(t) + v(t) (1b) where u(t) ∈ Rnu , y(t) ∈ Rny and x(t) ∈ Rn are respectively the input, the output and the state vectors. The process noise sequence is w(t) ∈ Rn and the output measurement disturbances are collected in v(t) ∈ Rny . (A, B, C, D) are the system matrices relatively to a certain coordinate basis of the state-space. It will be assumed that {x(t)}, {u(t)}, {y(t)}, {w(t)} and {v(t)}, t ∈ Z, are all ergodic and (weakly) stationary stochastic processes (Ljung, 1999), the system is asymptotically stable, controllable and observable and there is no feedback from y to u. The considered identification problem can then be stated N N as follows: given realizations {u(t)}t=1 and {y(t)}t=1 of the input and output processes generated by a system of the form (1) on a finite but sufficiently wide time horizon N , estimate the matrices (A, B, C, D) up to a similarity transformation. To reach this goal, two steps will be considered. Firstly, an input-output ARMAX model will be identified from the input-output measurements by using a classic least squares or instrumental variable method. Then, under the assumption that the matrix A is nonderogatory (Horn and Johnson, 1990), a canonical state-space form will be reconstructed from the estimated ARMAX model parameters. 3. ARMAX MODEL ESTIMATION Consider the output equation (1b). A direct way to get an input-output representation of the system from this statespace form is to get rid of the state. For that, let us start by expressing the state as a linear combination of the past inputs, outputs and a prior value of the state. Particularly, for t > n, it is easy to show that x(t) = An x(t − n) + ∆n un (t − n) + An wn (t − n) (2) with 1 ⊤ un (t − n) = u⊤ (t − n) · · · u⊤ (t − 1) ∈ Rnnu ∆n = An−1 B · · · AB B ∈ Rn×nnu 2 An = An−1 · · · A In ∈ Rn×n . Substituting Eq. (2) into Eq. (1b), we get y(t) = CAn x(t − n) + C∆n un (t − n) + CAn wn (t − n) + Du(t) + v(t). (3) Now, by using the Cayley-Hamilton theorem (Horn and Johnson, 1990), we know that An = −a0 In − a1 A − 1
yn (t − n) and wn (t − n) are built in a similar way.
· · · − an−1 An−1 where the {ai }i∈[0,n−1] elements are the coefficients of the characteristic polynomial of A, i.e., pA (z) = a0 + a1 z + · · · + an−1 z n−1 + z n . Thus, CAn x(t − n) = − a⊤ ⊗ Iny Γn x(t − n) (4) where (5) a⊤ = [a0 a1 · · · an−1 ] i⊤ h (6) Γn = C⊤ (AC)⊤ · · · CAn−1 ⊤
and with ⊗ the Kronecker product (Horn and Johnson, 1990). By using now the recurrent relations (1), it is straightforward to verify that the stacked vectors yn (t−n), un (t − n), vn (t − n) and wn (t − n) are related as follows yn (t − n) = Γn x(t − n) + Hn un (t − n) + Gn wn (t − n) + vn (t − n) (7) with 0 D 0 ··· 0 0 ··· 0 .. .. .. .. C . . . . . , G = Hn = CB n .. .. .. .. .. .. . .0 . . 0 . . CAn−2 B ··· CB D
CAn−2 ···
C 0
Eq. (7) allows us to express Γn x(t − n) as a combination of past inputs, outputs and disturbances, i.e., Γn x(t − n) = yn (t − n) − Hn un (t − n) − Gn wn (t − n) − vn (t − n). (8) Finally, gathering equations (3), (4) and (8), the state vector is eliminated in the output equation and we get y(t) = − a⊤ ⊗ Iny yn (t − n) + Fun+1 (t − n) + η(t) (9) with F = a⊤ ⊗ Iny Hn + C∆n D ∈ Rny ×(n+1)nu η(t) = a⊤ ⊗ Iny Gn + CAn wn (t − n) + a⊤ ⊗ Iny vn (t − n) + v(t). Note that F contain the coefficients of the polynomial matrix 2 that defines the numerator of the transfer function of (1) and a is the coefficient vector of the denominator. The model (9) is in fact quite more complex than a standard ARMAX since it contains at the same time a measurement and a process noise. It is also interesting to note that the noise part of this ARMAX model is composed of particular combinations of the AR and X deterministic parameters. Thus, as soon as a and F are consistently estimated, a model of the stochastic part can be built. In order to simplify the estimation of a and F, it is interesting to notice that Eq. (9) can be rewritten as a classic linear regression. For that, introduce the identity vec (AXB) = B⊤ ⊗ A vec (X) where vec (•) is the vectorization operator (Horn and Johnson, 1990). Then, a + η(t) y(t) = −Φyn (t − n) u⊤ (t − n) ⊗ I n y n+1 vec (F) (10) with Φyn (t − n) = [y(t − n) · · · y(t − 1)] ∈ Rny ×n . Thus, a and F can be estimated using a least squares algorithm such as the instrumental variable method (Ljung, 1999) when the data matrix 2 For an explicit expression of this polynomial matrix, we refer to the proof of Theorem 1 given in Appendix B.
1639
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
Φ(n) =
Φyn (1) .. .
u⊤ n+1 (1) ⊗ Iny .. .
Φyn (N − n) u⊤ n+1 (N − n) ⊗ Iny has full column rank.
requires that the matrix Φ(n) has full column rank. No matter how rich is the excitation input, it can be easily shown that rank(Φ(n)) ≤ nIO + (n + 1)nu ny in the noisefree case (i.e., when η(t) = 0 in Eq. (9)). It follows that it is necessary to satisfy nIO = n in order to have a full column rank matrix Φ(n).
4. FROM ARMAX MODEL TO MINIMAL CANONICAL STATE-SPACE REPRESENTATION 4.1 Condition to get ARMAX and state-space models of the same order The problem of parameterizing state-space models and the equivalence between state-space representation and MIMO canonical transfer functions is not new (see, e.g. (Kailath, 1980; McKelvey, 1995)). One of the main problems is the relation between the order nIO of the I/O model and the one of the corresponding minimal state-space representation. A glance at Eq. (9) and the way it is built from the state-space model (1) may tend to suggest that both representations have the same order. Unfortunately, this is not the case for all the MIMO systems (see (Stoica and Jansson, 2000) for a counterexample). Indeed, we have in general, nIO ≤ nSS ≤ nIO · ny , where nSS = n is the order of the state-space model (1). To see this fact more clearly, let us think of the order as the smallest integer m such that y(t) can be expressed in function of ym (t−m) and um+1 (t−m) as in Eq. (9). Then, intuitively, the order nIO of the input-output model (9) should be equal to the degree ν of the minimal polynomial of A. This remark follows essentially from Eq. (4) which suggests that nIO should be the smallest integer such that AnIO can be expressed in function of Aj , j = 0, . . . , nIO −1, where A is taken from a minimal realization of (1). On the other hand, if we consider the input-output transfer matrix of (1) to be TF(z) = L (z)/pA (z), with L (z) = Cadj(zI − A)B + pA (z)D f11 (z) · · · f1nu (z) .. = ... , . fny 1 (z) · · · fny nu (z) then the order nIO of Eq. (9) corresponds to o n j∈[1,n ] n − deg gcd {fij (z)}i∈[1,nyu] , pA (z) .
In fact, minimality of (A, B, C, D) is not enough to guarantee coprimeness between the polynomials fij (z) and pA (z). This is formally stated by the following theorem which relates the order of model (9) to that of (1). Theorem 1. Assume that the model (1) is minimal, that is, observable and reachable. Then the order nIO of the inputoutput model (9) is equal to the degree ν of the minimal polynomial of A, i.e., nIO = deg(mA (z)) = ν ≤ n. For the proof, see Appendix A. We know from Theorem 1 that, in general, the orders of the models (1) and (9) do not coincide, even when the model (1) is minimal. Therefore, even though it is assumed that the order n is known a priori, it is compulsory to determine nIO before being able to estimate a and F from (10). Indeed, consistent estimation of these parameters
To get equality between both orders, a special property on the matrix A in Eq. (1) is required as formulated by the following corollary of Theorem 1. Corollary 2. Assume that the model (1) is minimal. Then nIO = nSS = n if and only if the matrix A is nonderogatory. 3 In order to explicitly reconstruct a state-space representation from the estimated ARMAX parameters, it is essential to find conditions ensuring an equality relation between nIO and nSS . 4.2 Reconstruction of a minimal canonical state-space model Once the parameters of the model (9) are estimated, the reconstruction of a canonical state-space model can be considered. Before, the following lemma must be stated. Lemma 3. If the system (1) is observable, then there exists a matrix Ψ of the form ⊤ γ 0 ··· 0 .. 0 ... . ∈ Rn×nny , (11) Ψ= . . .. 0 .. 0 · · · 0 γ⊤ with γ ∈ Rny , and verifying rank (ΨΓn ) = n (12) where Γn is the observability matrix of order n, if and only if the matrix A is nonderogatory. A proof of this lemma can be found in (Bako et al., 2009). Note that the condition (12) in Lemma 3 can be interpreted as follows. There exists a γ such that if we let ya (t) = γ ⊤ y(t) be a linear combination of all the outputs, then the state x(t) of system (1) is entirely observable from ya (t). This further means that the system (1) and the MISO system with input u(t) and output ya (t) ∈ R have the same state. Lemma 3 tells us that this is in fact possible if and only if the matrix A is nonderogatory 4 . In view of that lemma, if, e.g, γ ∈ Rny is randomly generated from a uniform distribution, then the required condition (12) is fulfilled with probability one (Bako et al., 2009). By assuming now that the matrix A is nonderogatory (which ensures the equality nIO = nSS ), a minimal statespace realization can be readily obtained as explained in the following proposition: Proposition 4. Assume that the system (1) is minimal and that the matrix A is nonderogatory. Consider a matrix Ψ defined as in (11) and satisfying condition 3
A square matrix is said to be nonderogatory if its minimal polynomial coincides with its characteristic polynomial (Horn and Johnson, 1990). 4 In fact, for any SISO or MISO system to be observable the Amatrix needs to be nonderogatory.
1640
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 • an off-line version of the propagator method suggested initially in (Merc`ere et al., 2003) in a recursive form (see also (Merc`ere et al., 2008)) where a particular linear operator is introduced to rewrite the subspace-based identification problem into a quadratic optimization problem without constraint. For all these methods, the instrumental variable used to reduce the effect of the noise is built from the inputs of the system and the outputs of a simulated biased auxiliary model obtained with the ordinary MOESP algorithm (Verhaegen, 1994) applied on noisy data. Method in (Bako et al., 2008)
colnu ([ Y1
···
with Xi ∈ R
ny ו
and Yi ∈ R
, i ∈ [1, n].
As opposed to classic techniques leading to canonical forms for multivariable systems (Luenberger, 1967), (Mayne, 1972), (Kailath, 1980), the developed approach does not need any observability or controllability indices or an initial step leading to a non minimal realization estimate followed by a reduction procedure. Here is proposed a relatively straightforward method which gives access to minimal MIMO state-space representation from a standard vector difference equation. 5. SIMULATION EXAMPLE In order to illustrate the performances of the algorithm, the following system is considered # # " " C=
0.8736 0 0 −0.3881 , 0 0
−0.3814 0.6134 −0.3495 0.4044 −0.0624 −0.7160 , B = −0.5787 −0.5476 −0.1790
h
i
0.9397 0 1.1787 , 0 0 −1.3274
D=
h
−0.1
−0.5
0 0.5 Real part
1
i
0.5463 −0.5293 . 0 −2.4003
The input is a zero-mean white Gaussian noise with unit variance. The disturbances v and w are both zero-mean white Gaussian noise set such that the noise signal ratio equals 20dB in both cases. A Monte Carlo simulation with 100 realizations of noise and input is carried out. This new algorithm is compared with three other methods: • another SVD-free identification method developed in (Bako et al., 2008) where the subspace identification problem is recasted into a simple least squares problem by introducing a particular user-defined matrix which does not change the rank of the extended observability matrix when multiplying this latter matrix on the left hand side, • the traditional PO MOESP algorithm (Verhaegen, 1994),
0.1 0 −0.1
−1
−0.5
0 0.5 Real part
1
Method in (Merc` ere et al., 2008)
PO-MOESP
0.3 0.2
0.1 0 −0.1 −0.2 −1
Our algorithm
−0.2
0.2
For the proof, see appendix B.
A=
0
0.3
Yn
•×nu
0.2
0.1
Imaginary part
Xn ] ,
Xn
0.2
−1
"Y # 1 .. , Yn ]) = . Imaginary part
rowny
0.3
−0.2
n−1
1 0 ··· 0 " X #! 1 .. = [ X1 ··· .
0.3 Imaginary part
Imaginary part
(12). Then a minimal canonical state-space representation can be obtained from the estimated parameters a⊤ = [a0 · · · an−1 ] and F of the ARMAX model (9) as: 0 1 ··· 0 .. . . .. . 0 , . . (13) A= 0 0 ··· 1 −a0 −a1 · · · −an−1 D = F (:, nnu + 1 : (n + 1)nu ) , (14) B = ΨQ, (15) −1 ⊤ C = rowny (Q) Ω⊤ , (16) n Ωn Ω n where −1 colnu F(:, 1 : nnu ) − a⊤ ⊗ D , Q = K ⊗ Iny Ωn = B AB · · · An−1 B ∈ Rn×nnu , a1 · · · an−1 1 .. . 1 0 ∈ Rn×n , K= a
0.1 0 −0.1 −0.2
−0.5
0 0.5 Real part
1
−1
−0.5
0 0.5 Real part
1
Fig. 1. Estimated (o) and real poles (+) of the system. Figure 1 compares the position of the estimated poles of the models given by the four aforementioned algorithms with respect to the real system. A first conclusion is that all these methods are quite accurate since the estimated poles are centered on the true ones and the variances of the estimates are relatively small. We can however notice that, on this example, the new algorithm leads, without requiring the use of an SVD, to similar results as the PO MOESP algorithm. Therefore, the approach developed in this article seems to be an interesting alternative to the classic subspace-based techniques. 6. CONCLUSION In this paper, a simple SVD-free method has been proposed to get a state-space model, given a collection of its input-output measurements. In comparison with existing subspace identification techniques, the main motivation of this new method is to prevent the state representation basis from shifting according to the estimation data. This feature is highly desirable in, e.g., recursive identification of state-space models or in the estimation of composite models with linear submodels. From a realization theory point of view, this work highlights the behavior of certain types of linear MIMO systems that can indeed be handled as MISO systems in terms of state dynamics. Future research will focus on generalizing this procedure to more arbitrary models than the particular cases where the Amatrix is nonderogatory.
1641
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 j∈[1,n ] gcd {lij (z)}i∈[1,nyu] , pA (z)
Appendix A. PROOF OF THEOREM 1 We first state the following lemma. Lemma 5. Let A ∈ Rn×n and B = [b1 · · · bm ] ∈ Rn×m be some arbitrary matrices. Let mA (z) be the mini mal polynomial of A and define An−1 B · · · AB B = ∆n (A, B) ∈ Rn×nm , A b1 .. .. nm×nm ¯ ¯ A= , b = . ∈ Rnm . . ∈R A
bm
Then the following holds. If rank ∆n (A, B) = n, then ¯ = ν, where ν = deg(mA (z)). ¯ b) rank ∆n (A, Proof. Assume that rank ∆n (A, B) = n. We then need ¯ b) ¯ is equal to to show that the rank of the matrix ∆n (A, the degree ν of the minimal polynomial of A. We proceed by contradiction. Assume that there exists j < ν such ¯ , ¯ is a linear combination of the vectors A ¯ ib ¯ jb that A i = 0, . . . , j − 1. Then, one can find a set of scalars ¯+ ¯ = αj−1 A ¯ j−1 b ¯ jb α0 , · · · , αj−1 not all zero satisfying A ¯ · · · + α0 b. By Denoting with p the polynomial of degree j defined by p(z) = z j − αj−1 z j−1 − · · · − α0 , we get ¯ = 0. This is equivalent to p(A)bi = 0, i = 1, . . . , m. ¯ b p(A) As a consequence, p(A)B = 0 and p(A)∆n (A, B) = 0 because, as a polynomial of A, p(A) commutes with any power of A. Since ∆n (A, B) has full row rank, it follows that p(A) = 0. Therefore, if we call mA (z) the minimal polynomial of A, then, by definition of mA (z) we get that mA (z) divides p(z). This contradicts the fact that deg {p(z)} = j < deg {mA (z)} = ν. In conclusion, ¯ b) ¯ = ν. rank ∆n (A,
From this lemma, it can be concluded that under control¯ = n if ¯ b) lability assumption of system (1), rank ∆n (A, and only if the A-matrix is nonderogatory. For the proof of Theorem 1, we will need also the following Theorem given in (Barnett, 1973): Theorem 6. Let bi (z) = bi1 z n−1 + · · · bin , i = 1, · · · , p, pA (z) = a0 + a1 z + · · · + an−1 z n−1 + z n . Then the degree of the gcd of pA (z), b1 (z), · · · , bp (z) is equal to n − rank (Γn (Ac , M)) with 0 1 ··· 0 b1n · · · b11 .. . . . .. , A = . . .. . M = ... c . 0 0 ··· 1 bpn · · · bp1 −a0 −a1 · · · −an−1 With this material, Theorem 1 can be proved as follows:
Proof. The input-output transfert matrix is given by Cadj(zI − A)B + pA (z)D C(zI − A)−1 B + D = , pA (z) where adj(·) is the cofactor matrix of (·) and pA (z) = det(zI − A) is the characteristic polynomial of A. Let us introduce the following notation j∈[1,n ]
Cadj(zI − A)B = [lij (z)]i∈[1,nyu] , where lij (z), the numerator polynomial that corresponds to the transfert between the i-th output and the j-th input. To prove Theorem 1, we need to show that the degree of
is equal to n − ν, with ν = deg(mA (z)). To proceed, we start by pointing out that, as far as we are concerned with seeking the common factors between the polynomials lij (z) and pA (z), we can, without loss of generality, assume that D = 0. Now we compute the coefficients of the polynomials lij (z). The expressions of these coefficients can be directly obtained from the term F of Eq. (9) (i.e., the matrix that multiplies un (t−n)). At this step, it is important to notice that F = (a⊤ ⊗ Iny )Hn + C∆n = C∆n (L ⊗ Inu ), with 1 0 · · · 0 .. . an−1 1 L= . ∈ Rn×n . . ... ... 0 . a1 · · · an−1 1 Let bj be j-th column of B and c⊤ i be the i-th row of C. The n coefficients of the n − 1 degree polynomial lij (z) are given by the following row vector n−1 lij = c⊤ bj · · · Abj bj L. i A | {z } ∆n (A,bj )
Let
and then
l1j Nj = ... = C∆n (A, bj )L lny j
C∆n (A, b1 ) N1 .. nn ×n N = ... = L ∈ R y . C∆n (A, bnu ) Nnu as a concatenation of all the coefficient vectors of the polynomial matrix Cadj(zI− A)B. According to Theorem 6, we just need to show that rank (P) = ν, where N NAc P= .. .
NAcn−1 and Ac a companion matrix (see Eq. (13)). Notice that the matrix P is equal to C∆n (A,b1 )L C∆n (A,b1 ) C∆n (A,b1 )AL C∆n (A,b1 )LAc .. .. . . n−1 C∆n (A,b1 )LAn−1 C∆n (A,b1 )AL c .. .. = P1 = L, . . C∆n (A,bnu )L C∆n (A,bnu ) C∆n (A,bn )LAc C∆n (A,bn )AL u u .. .. . . n−1 n−1 C∆n (A,bnu )LAc
C∆n (A,bnu )AL
up to some rows permutations, where −an−1 1 .. .. . . AL = LAc L−1 = −a 0
1642
1
−a0
··· .. . ··· 0 ···
0 .. . . 1
0
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 For any k ≥ 0 and any bj , we have ∆n (A, bj )AkL = Ak ∆n (A, bj ). By making use of this fact, it follows that Γn .. ¯ ¯ b)L, P1 = ∆n (A, . Γn
with Γn = Γn (A, C). Since the system is observable, the matrix that contains Γn in the right hand side of the previous equation, has full column rank. Therefore, Lemma 5 allows us to conclude that rank P = rank P1 = ¯ = ν since L is nonsingular. ¯ b) rank ∆n (A, Appendix B. PROOF OF PROPOSITION 4
The reconstruction of a canonical form from the parameters a and F is in fact only based on shrewd calculations. For D, it is easy to see that D = F (:, nnu + 1 : (n + 1)nu ) . Knowing D, it is interesting to notice that F (:, 1 : nnu ) − a⊤ ⊗ D = a⊤ ⊗ Iny Hn + C∆n − a⊤ ⊗ D ⊤ = k⊤ 1 ⊗ Iny Γn B · · · kn ⊗ Iny Γn B
where k⊤ i is the ith row of K. Thus, by using the operator colnu defined in Proposition 4 and noticing that K ⊗ Iny is invertible, we get −1 Γn B = K ⊗ Iny colnu F (:, 1 : nnu ) − a⊤ ⊗ D . Introducing matrix Ψ, it holds that ¯ ΨΓn B = TB = B −1 = Ψ K ⊗ Iny colnu F (:, 1 : nnu ) − a⊤ ⊗ D since ΨΓn = T ∈ Rn×n can be viewed as a similarity transformation if rank (ΨΓn ) = n. The next step consists in reconstructing the state matrix A. This can be done by ¯ = TAT−1 knowing T, i.e. calculating A ¯ = ΨΓn A (ΨΓn )−1 A ⊤ γ CA CA −1 −1 .. = Ψ ... (ΨΓn ) = (ΨΓn ) . . CAn γ ⊤ CAn Now, using again the Cayley-Hamilton theorem, we have 0 ⊤ 1 ··· 0 ⊤ γ C γ CA .. .. . . .. .. . . . = . . 0 0 · · · 1 ⊤ n−1 ⊤ n γ CA γ CA −a0 −a1 · · · −an−1 = Ac ΨΓn . where Ac is a companion matrix (see Eq. (13)). Hence, ¯ = Ac ΨΓn (ΨΓn )−1 = Ac . A ¯ B ¯ and Γn B, it is easy to see that Finally, knowing A, ¯ =C ¯ B ¯ nB ¯ n−1 B ¯ ¯ ··· A rowny (Γn B) = rowny Γ ¯ which leads to the extraction of C.
REFERENCES
L. Bako, G. Merc`ere, and S. Lecœuche. On-line structured subspace identification with application to switched linear systems. International Journal of Control, 2009. Accepted for publication. S. Barnett. Matrices, polynomials and linear time invariant systems. IEEE Transactions of Automatic Control, 18:1–10, 1973. W. Favoreel, B. De Moor, and P. Van Overschee. Subspace state space system identification for industrial processes. Journal of Process Control, 10:149–155, 2000. R. Horn and C. Johnson. Matrix analysis. Cambridge University Press, 1990. T. Kailath. Linear Systems. Prentice Hall, Engelwood Cliffs, 1980. T. Katayama. Subspace methods for system identification. Springer, 2005. L. Ljung. System identification. Theory for the user. Prentice Hall Information and System Sciences Series, Upper Saddle River, 2nd edition, 1999. M. Lovera. Identification of MIMO state space models for helicopters dynamics. In Proceedings of the 13th IFAC Symposium on System Identification, Rotterdam, The Netherlands, August 2003. M. Lovera, T. Gustafsson, and M. Verhaegen. Recursive subspace identification of linear and non linear Wiener state space models. Automatica, 36:1639–1650, 2000. D. Luenberger. Canonical forms for linear multivariable systems. IEEE Transactions on Automatic Control, 12: 290–293, 1967. D. Mayne. A canonical model for identification of multivariable linear systems. IEEE Transactions on Automatic Control, 17:728–729, 1972. T. McKelvey. Identification of state space models from time and frequency data. PhD thesis, Link¨oping University, Link¨oping, Sweden, 1995. G. Merc`ere, S. Lecœuche, and C. Vasseur. A new recursive method for subspace identification of noisy systems: EIVPM. In Proceedings of the 13th IFAC Symposium on System Identification, Rotterdam, The Netherlands, August 2003. G. Merc`ere, L. Bako, and S. Lecœuche. Propagator-based methods for recursive subspace model identification. Signal Processing, 88:468–491, 2008. L. Mevel, L. Hermans, and H. Van der Auweraer. On the application of subspace-based fault detection methods to industrial structures. Mechanical Systems and Signal Processing, 13:823–838, 1999. P. Stoica and M. Jansson. MIMO system identification: state-space and subspace approximations versus transfer function and instrumental variables. IEEE Transactions on Signal Processing, 48:3087–3099, 2000. R. T´oth, F. Felici, P. Heuberger, and P. Van den Hof. Discrete-time LPV I/O and state-space representations. differences of behavior and pitfalls of interpolation. In Proceeding of the European Control Conference, Kos, Greece, July 2007. M. Verhaegen. Identification of the deterministic part of MIMO state space models given in innovations form from input output data. Automatica, 30:61–74, 1994.
L. Bako, G. Merc`ere, and S. Lecœuche. A least squares approach to the subspace identification problem. In In the 47th IEEE Conference on Decision and Control, Cancun, Mexico, December 2008.
1643