A fast algorithm for subspae state-spae system identiation via exploitation of the displaement struture Niola Mastronardi a d 2, Daniel Kressner b 3, Vasile Sima a e 4, Paul Van Dooren 3 5, and Sabine Van Huel a 1 2 3 ; ; ; ; ; ; ; ; ; ; Department of Eletrial Engineering, ESAT-SISTA/COSIC, Katholieke Universiteit Leuven, Kardinaal Merierlaan 94, 3001 Heverlee, Belgium b Department of Mathematis, University of Chemnitz, Chemnitz, Germany Department of Mathematial Engineering, Universite Catholique de Louvain, Avenue Georges Lemaitre 4 B-1348 Louvain-la-Neuve, Belgium d Dipartimento di Matematia, Universita della Basiliata, via N. Sauro 85, 85100 Potenza, Italy e National Institute for Researh & Development in Informatis, Bd. Maresal Averesu, Nr. 8{10, 71316 Buharest 1, Romania a Abstrat Two reent approahes [15,16℄ in subspae identiation problems require the omputation of the R fator of the QR fatorization of a blok{Hankel matrix H , whih, in general has a huge number of rows. Sine the data are perturbed by noise, the involved matrix H is, in general, full rank. It is well known that, from a theoretial point of view, the R fator of the QR fatorization of H is equivalent to the Cholesky fator of the orrelation matrix H T H , apart from a multipliation by a sign matrix. In [12℄ a fast Cholesky fatorization of the orrelation matrix, exploiting the blok{Hankel struture of H , is desribed. In this paper we onsider a fast algorithm to ompute the R fator based on the generalized Shur algorithm. The proposed algorithm allows to handle the rank{deient ase. Key words: generalized Shur algorithm, Hankel and blok{Hankel matries, subspae identiation, QR deomposition, singular value deomposition. Preprint submitted to Elsevier Preprint 1 Introdution Subspae based system identiation has beome very popular in the last deade [4℄. The suess of this state{spae identiation approah is mainly due to the fat that it relies on a simple matrix deomposition for whih reliable numerial algorithms are available. Its major drawbak, on the other hand, is that large \data" matries are involved and that this may lead to high omputing and storage osts. We now briey reall the basi formulation of the problem. Let uk and yk be the m{dimensional input vetor and the l{dimensional output vetor, respetively, of the linear time{invariant state{ spae model xk+1 = Axk + Buk + wk ; yk = Cxk + Duk + vk ; where xk is the n{dimensional state vetor at time k, fwk g and fvk g are state and output disturbanes or noise sequenes, and A, B , C and D are unknown real matries of appropriate dimensions. For non{sequential data proessing, one hooses N 2(m + l)s and onstruts the N 2(m + l)s matrix H = [U Ts;N Y Ts;N ℄, where U s;N and Y s;N are blok{ 2 2 2 2 This author is a Senior Researh Assoiate with the F.W.O. (Fund for Sienti Researh-Flanders). 2 This work is supported by UE Programme \Training and Mobility of Researhers" projet (ontrat ERBFMRXCT970160) entitled \Advaned Signal Proessing of Medial Magneti Resonane Imaging and Spetrosopy". 3 This work is supported by the Belgian Programme on Interuniversity Poles of Attration (IUAP-4/2 & 24), initiated by the Belgian State, Prime Minister's OÆe for Siene, and by a Conerted Researh Ation (GOA) projet of the Flemish Community, entitled \Mathematial Engineering for Information and Communiations Systems Tehnology". 4 This work is supported by the European Community BRITE-EURAM III Themati Networks Programme NICONET (projet BRRT{CT97-5040). 5 This work is partly supported by the National Siene Foundation ontrat CCR97-96315. 1 2 Hankel matries dened in terms of the input and output data, respetively : 2 U2s;N = 6 6 6 6 6 6 6 6 6 6 6 6 4 2 Y2s;N = 6 6 6 6 6 6 6 6 6 6 6 6 4 u1 u2 u2 u3 u3 u4 ... ... u3 : : : uN u4 : : : uN +1 u5 : : : uN +2 ... ... ... 7 7 7 7 7 7 7; 7 7 7 7 7 5 ... u2s u2s+1 u2s+2 : : : uN +2s y1 y2 y2 y3 y3 y4 3 1 3 y3 : : : yN y4 : : : yN +1 y5 : : : yN +2 ... 7 7 7 7 7 7 7: 7 7 7 7 7 5 ... y2s y2s+1 y2s+2 : : : yN +2s 1 Then the R fator of a QR fatorization H = QR is used for data ompression. In [12℄ a fast Cholesky fatorization of the orrelation matrix, exploiting the blok{Hankel struture of H , is desribed. In this paper we onsider a fast algorithm to ompute the R fator of the QR fatorization of H based on the generalized Shur algorithm, exploiting its displaement struture. The paper is organized as follows. In x2 the generalized Shur algorithm to ompute the Cholesky fator of a symmetri positive denite matrix is desribed and in x3 this algorithm is applied to the matrix H . The rank{deient ase is desribed in x4 and some numerial experiments are reported in x5. 2 The Shur algorithm for positive semi-denite matries 2.1 The Cholesky fator and generator of A We summarize here the key properties of the generalized Shur algorithm to ompute the Cholesky fator of a (symmetri) positive semi-denite (psd) matrix, whih will be used in the next setion. More details an be found in [7,8℄. Let A be a psd matrix of order n, then we dene its displaement rA = A Z T AZ; using a generalized shift matrix Z of the same dimension. Here we only require the shift matrix Z to be stritly upper triangular (and hene nilpotent) and 3 we speialize later on to a partiular hoie of Z . We all the rank of rA the displaement rank of A and we assume that it is signiantly smaller than n. Let the symmetri matrix rA have p positive eigenvalues and q =: p negative eigenvalues then it has a fatorization 2 0 I rA =: GT G; =: 64 p 0 Iq 3 7 5; 2 G =: 64 3 Gp 7 5: Gq (1) The matrix G is alled the generator of A and sine Z is nilpotent one an reonstrut A via the formula A= nX1 i=0 (Z i)T GT GZ i : The generator and intermediate results derived from transformations of the generator, allow to reonstrut the Cholesky fator of the psd matrix A : 3 2 6 r1;1 r1;2 : : : rn;n 7 7 6 7 6 r 6 2;2 : : : r2;n 7 7 6 A = RT R; R = 6 . . . ... 6 6 4 rn;n 7: 7 7 5 Notie that if A has rank r < n then so will the fator R whih will have its last n r rows equal to zero : 2 A= RT R; R = 6 r1;1 6 6 6 6 4 3 : : : : : : rn;n 7 ... ... rr;r : :rr;n 7 7 7: 7 5 If the leading r r prinipal submatrix of A is nonsingular then ri;i; i = 1; : : : ; r are all nonzero. Otherwise the \prole" of the trapezoidal matrix R indents to the right eah time the nullity of the i i prinipal submatrix of A inreases. In our appliation, A is a produt of the type H T H whih learly is a psd matrix. Its Cholesky fator R is, up to a sign matrix D = diag(1; : : : ; 1), also the RH fator of the QH RH deomposition of H : RH = DR. Hene both the problem of omputing the RH fator of the QR fatorization of H and that of omputing the Cholesky fator of H T H are equivalent. We disuss now the omputation of the Cholesky fator R of A starting from the generator of rA. 4 One easily shows that the generator G is not unique. We say that the generator G~ is proper if its rst olumn is zero exept possibly its leading element. The following theorem holds for proper generators [7℄. Theorem 1 Let 3 2 a11 a12 7 A = 64 5 a21 A22 be a positive semi-denite matrix with proper generator 2 3 2 6 g~1;1 g~1;2 g~1;n 7 6 6 7 6 6 7 6 6 0 g~2;2 g~2;n 7 : 6 ~G = 66 . . 7=6 .. 77 66 .. 6 .. . 7 6 6 4 5 4 0 g~;2 g~;n G~ 1 G~ 2 3 7 7 7 7 7; 7 7 7 5 then G~ 1 is the rst row of the Cholesky fator R of A. Furthermore the gen^ where erator matrix for the Shur omplement A^ = A G~ T1 G~ 1 is given by G; 2 60 A^ = 6 6 6 6 6 60 6 4 0 A22 a21 a a12 1 11 3 2 7 7 7 7 7; 7 7 7 5 6 6 6 6 6 6 6 6 4 G^ = G~ 1 Z G~ 2 3 7 7 7 7 7: 7 7 7 5 We observe that the rst olumn of G^ is zero, whih needs to be the ase sine the rst olumn and row of A^ are zero. The generalized Shur algorithm just onsists of a reursive use of this Theorem : via a transformation (dened below) the generator G of the urrent ~ This yields the urrent row of the Cholesky matrix A is put in proper form G: fator and the generator of the Shur omplement is trivially obtained from a shift Z applied to the rst row of the generator. We refer to [7℄ for more details. The omplexity of this algorithm is that of the transformation sine the shift Z does not imply any operations. In the next setion we desribe briey the onstrution of . 5 2.2 Redution of the generator to proper form The rst row of the Cholesky fator R of A is thus obtained from a proper generator G of A. Reduing a non-proper generator of A to a proper generator ~ is obtained by applying a transformation to the generator G. In order G; not to hange the produt GT G it suÆes to hoose to be -unitary, i.e. : T = ; sine then G~ T G~ = (G)T (G) = GT G: Typially the matrix is onstruted as follows : 2 =: 6 6 6 6 Ip 6 6 6 6 4 3 2 7 6 7 6 7 6 7 6 7:6 7 6 7 6 7 6 5 6 4 1 Iq 3 Hp 1 7 7 7 7 7 7: 7 7 7 Hq 5 (2) The bloks Hp and Hq of the seond fator are p p and q q Householder transformations reduing the rst olumn of G as follows 2 2 6 Hp 4 32 Hq 76 54 3 Gp Gq 7 5 = 6 x11 6 6 6 6 6 6 6 6 6 6 6 6 6 6 y11 6 6 6 6 6 6 6 6 4 x12 0 ... X 0 22 y12 0 ... Y 0 6 22 3 7 7 7 7 7 7 7 7 7 7 7 7 7: 7 7 7 7 7 7 7 7 7 7 5 (3) The rst fator :only transforms the rows ontaining x and y and eliminates y provided = y =x is smaller than one in modulus : 11 11 11 11 11 3 3 2 2 3 2 ~11 x~ 12 7 6x 6 7 6 x11 x12 7 5: 5=4 5:4 4 0 y~12 y11 y12 This 2 2 transformation is onstruted from =: 1=p1 , =: =p1 and is alled aphyperboli rotation sine it satises = 1 [1℄. Also note that x~ = x 1 . When it is implemented in fatored form : 2 2 11 2 2 2 11 2 6 4 1 0 p 1 2 32 76 56 4 p 32 3 0 77 6 1 7 5; 54 0 1 0 1 1 1 2 one shows that the generalized Shur algorithm is bakward stable and that it has the same omplexity as the unfatored implementation [13℄. It follows already from (3) that the (1; 1) element of the psd matrix A equals a =x y 0. Therefore if a 6= 0 the above transformation an be performed. On the other hand, if a = 0 then the whole row a must be zero sine otherwise A would not be psd. Sine a = x x y y this also implies that [x ; x ℄ = [y ; y ℄ and that both these rows an just be deleted from the generator [5℄. In other words, if a = 0 a simpliation an be performed to the urrent generator G~ . The omplexity of the redution of G to proper form G~ is essentially that of the Householder transformations Hp and Hq whih osts 4(p + q)n ops. If r = rank A steps are performed, this algorithm thus requires a total of 4r(p + q)n ops. We point out that this is an overestimate sine the number of nonzero olumns n of the generator dereases at eah step and that potentially the number of rows p + q may derease as well. 2 11 11 2 11 11 11 12 12 11 12 11 11 12 11 12 12 11 3 Fast omputation of the R fator of the QR fatorization of H We show here how to ompute the gererator G of A = H T H where H 2 R N; m l s is the blok Hankel matrix desribed in the rst setion with bloks 2( + ) 7 of sizes 1 m and 1 l : 2 uT 6 1 6 6 T 6 u2 6 6 6 6 4 uTN uT 2 uT2s ::: yT ... ... ... ... H= . . . .. . . . . 1 yT 2 ::: 2 1 7 7 7 7 7: 7 7 7 5 ... ... ... ... ... . . . . . . yT : : : : : : uTN +2s 3 y2Ts yNT : : : : : : yNT +2s 1 The shift matries used in this ontext are the matries 6 0m Im 6 6 0m 6 : Zm = 66 6 6 4 ... . . . Im 0m 3 2 3 2 6 0l Il 7 6 7 6 7 0l 6 7 : 7 ; Zl = 6 6 7 6 7 6 4 7 5 7 . . . 777 : . . . Il 777 ; Z = Zm Zl : (4) 7 5 0l The following theorem then gives a ontrution of a generator for A. Theorem 2 Given the QR fatorization of the rst blok olumns : 2 T 6 u1 6 6 T 6 u2 6 6 .. 6 . 6 4 uTN 2 3 3 y1T 7 6 q1T 7 7 6 7 y2T 777 666 q2T 777 .. 77 = 66 .. 77 R1 . 7 6 . 7 5 4 5 yNT qNT (5) where R1 an be assumed upper trapezoidal of row rank k m + l and qi 2 R k ; dene the produt Cu;1 : : : Cu;2s Cy;1 : : : Cy;2s = q : : : qN H: 1 (6) Then a generator G for H T H is given by = Ik Ik G = Gu Gy ; +1 8 +1 (7) where 2 Gu = 6 Cu;1 6 6 6 6 6 6 6 4 0 0 0 2 3 3 Cu;2 : : : Cu;2s 7 6 Cy;1 Cy;2 : : : Cy;2s 7 6 7 7 6 0 yT 7 T uTN +1 : : : uTN +2s 1 777 : : : y 6 N +1 N +2s 1 7 6 7: ; G = y 6 7 7 6 7 7 Cu;2 : : : Cu;2s 7 6 0 Cy;2 : : : Cy;2s 7 4 5 5 uT1 : : : uT2s 1 0 y1T : : : y2Ts 1 Proof : In order to prove the result we onsider the displaement matrix rH T H : 2 T 6 U2s;N U2s;N 4 Y2s;N U2Ts;N ZmT U2s;N U2Ts;N Zm U2s;N Y2Ts;N ZlT Y2s;N U2Ts;N Zm Y2s;N Y2Ts;N 3 ZmT U2s;N Y2Ts;N Zl 7 5 ZlT Y2s;N Y2Ts;N Zl (8) whih ought to be equal to 3 2 T T 6 Gu Gu Gu Gy 7 5: 4 GTy Gu GTy Gy (9) It follows from (5,6) that R = [Cu; Cy; ℄ and hene 1 1 1 [Cu; Cy; ℄ Cu; : : : Cu; s Cy; : : : Cy; s = 1 1 T 1 2 1 2 RT 1 q1 : : : qN H; whih are the rst blok rows of the sub-bloks of (8). This thus veries the rst blok rows and blok olumns of the equality between (8) and (9). The rest easily follows from the blok Hankel struture of H . Note that if the rst blok olumns of H in (5) have full rank then R is square invertible and k = m + l. If moreover the whole matrix H has full olumn rank, then the generalized Shur algorithm will not enounter any singularities. But sine the low rank ase is of partiular interest here, singularities in the generalized Shur algorithm will be enountered and hene lead to a lower omplexity of the algorithm. The above theorem also shows that the displaement rank of H T H is at most 2(k + 1) 2(m + l + 1), with the same number of positive and negative generators. Hene the generalized Shur algorithm to ompute the R fator requires about (8Nrk) ops. To ompute the generator G of H T H , the QR 1 9 Table 1 Numerial results for Example 1 RM RS bakward error RM bakward error RS numerial rank # ops # ops 31660 1:51 10 9489 5:20 10 16 5 15 fatorization (5) requires (6N (m + l) ) ops and the produt (6) requires less than (4Nk(m + l)s) ops. We reall that k (m + l) and r 2(m + l)s but that equality is obtained when no rank deieny is deteted. The most time onsuming steps are then learly the generalized Shur algorithm and the produt (6). 2 4 The generalized Shur algorithm for rank{deient matries Our desription of the generalized Shur algorithm allows to handle rank deient matries H T H: In this ase we an drop some rows: of the generator during the algorithm. For this we need a tolerane, say Æ = kH T H k where is the requested relative auray. Referring to the desription of setion 2, we test if x y Æ: We then hek as well if the leading row a of the urrent Shur omplement is small. If so, the urrently omputed row of the Cholesky fator is negletable and we delete the two orresponding rows of the generator. It is possible that a is muh larger than Æ although a Æ: In this ase the deletion of a row of the Cholesky fator will yield residual errors kH T H RT Rk of the same size. This is analyzed in this setion. From the rst example we an onlude that the desribed proedure works aurately when it is applied to a matrix H with a suÆiently large gap between signiant singular values and negligible ones. On the other hand, a loss of auray in the omputed fator R is observed when the distribution of the small singular values of H shows a uniform and slow derease. The relative auray is hosen equal to 10 in both examples. Example 1 Consider the matrix H = [U T jY T ℄; with Y = U , where the rst 2 11 2 11 12 12 11 13 row and the last olumn of U are [ 40 39 38 : : : 3 2 1 2 2 3 ℄; [ 3 2 2 1 2 3 4 5 6 7 ℄T ; respetively. The rank of the matrix H is 5 and kH T H k = 6:31 10 : 1 5 In Table 1 the results of the omputation of the R fator of the matrix H by means of the standard QR and the generalized Shur algorithm are shown. We denote by RM ; RS ; bakward error R ; and numerial rank, the R fator of 10 5 10 0 10 −5 10 −10 10 −15 10 −20 10 −25 10 −30 10 −35 10 0 2 4 6 8 10 12 14 16 18 20 Fig. 1. Distribution of the singular values, in logarithmi sale, of the matrix on- sidered in Example 1 Table 2 Numerial results for Example 2 R M # ops 12:49 106 R S bakward error RM bakward error RS numerial rank # ops 40:61 104 1:27 10 2:92 10 14 2 18 the QR fatorization of H omputed by the matlab funtion triu(qr(H)) and by the generalized Shur algorithm, the bakward error of H T H dened as kH T H RT Rk ; kH T H k 1 1 and the rank of H deteted by the generalized Shur algorithm, respetively. In this ase, the R fator is aurately omputed by the generalized Shur algorithm, beause of the big dierene between the signiant singular values and the negligible ones of H (see Figure 1). Example 2 This is the fourth appliation onsidered in the next setion. In Figure 2 we an see that the distribution of the small singular values of the involved matrix H slightly dereases. We point out that the orrelation matrix H T H omputed by matlab is not numerially s.p.d. beause of the nearly rank deieny of H: Furthermore kH T H k1 = 3:99 104 : So, in this ase the fast Cholesky fatorization, exploiting the blok{Hankel struture of H and desribed in [12℄, an not be used. In Table 2 we an see that, although the generalized Shur algorithm is very fast w.r.t. the standard QR algorithm, the ahieved auray is not satisfatory. 11 3 10 2 10 1 10 0 10 −1 10 −2 10 −3 10 −4 10 −5 10 −6 10 −7 10 0 10 20 30 40 50 60 70 80 Fig. 2. Distribution of the singular values, in logarithmi sale, of the matrix on- sidered in Example 2 5 Numerial results In this setion results omputing the R matrix by means of the generalized Shur algorithm are summarized. The data sets onsidered are publily available on the DAISY web site http://www.esat.kuleuven.a.be/sista/daisy At eah iteration of the generalized Shur algorithm, two Householder matries and one modied hyperboli rotation are omputed in order to redue the generator in proper form. All the numerial results have been obtained on a Sun workstation Ultra 5 using Matlab 5.3. Table 3 gives a summary desription of the appliations onsidered in our omparison, indiating the number of inputs m; the number of outputs l; the number of blok rows s; the total number of data samples used t and the number of rows of H . In Table 4 some results for the omputation of the R fator of the QR fatorization of H are presented. Rel. residual denotes kjRM j jRS jk : kjRM jk 1 1 The results in Table 4 are omparable with those desribed in [12℄, where the R fator is obtained onsidering the Cholesky fatorization of the orrelation matrix H T H , and exploiting the blok-Hankel struture of H . The analysis of the fourth appliation is desribed in Example 2 of the previous setion. 12 Table 3 Summary desription of appliations. Appl. # Appliation m l s t N 1 Glass tubes 2 2 20 1401 1361 2 Labo dryer 1 1 15 1000 3 Glass oven 3 6 10 1247 1227 4 Mehanial utter 1 1 20 1024 960 5 Flexible robot arm 1 1 20 1024 984 6 Evaporator 3 3 10 6305 6285 7 CD player arm 2 2 15 2048 2018 8 Ball and beam 1 1 20 1000 9 Wall temperature 2 1 20 6800 1640 Table 4 Comparative results for the omputation of the Appl. # Appliation 1 Glass tubes 2 Labo dryer 3 Glass oven 4 Mehanial utter R M R R S # ops 6:76 107 2:61 106 7:63 107 1:25 107 5 Flexible robot arm 1:25 107 6 Evaporator 7 CD player arm 8 Ball and beam 9 Wall temperature 1:82 108 5:77 107 1:22 107 4:67 107 960 fator. # ops 7:01 106 970 3:36 105 6:38 106 4:89 103 4:76 105 1:11 107 2:59 106 4:67 105 1:64 106 bakward error RS 2:20 10 8:30 10 3:73 10 2:92 10 4:13 10 6:26 10 5:01 10 7:59 10 2:45 10 15 15 15 2 15 15 15 15 14 Rel.residual 8:73 10 9:48 10 7:91 10 13 12 0:39 100 3:38 10 5:13 10 2:08 10 6:79 10 3:76 10 6 Conlusions In this paper the generalized Shur algorithm to ompute the R fator of the QR fatorization of blok{Hankel matries, arising in some subspae identiation problems, is desribed. It is shown that the generalized Shur algorithm is signiantly faster than the lasssial QR fatorization. A rank{revealing implementation of the generalized Shur algorithm in ase of rank{deient matries is also disussed. 13 14 5 14 8 13 12 Algorithmi details and numerial results have been presented. Referenes [1℄ A.W. Bojanzyk, R.P. Brent, P. Van Dooren and F.R. De Hoog, A note on downdating the Cholesky fatorization, SIAM J. Si. Stat. Comput., 1 (1980), pp. 210{220. [2℄ Y.M. Cho, G.Xu and T. Kailath, Fast reursive identiation of state spae models via exploitation of displaement struture, Automatia, Vol.30. No.1 (1994), pp. 45{59. [3℄ [4℄ [5℄ , Fast parallel algorithms for QR and triangular fatorization, SIAM J. Si. and Stat. Comp., Vol. 8, No. 6 (1987), pp. 899{913. J. Chun, T. Kailath and H. Lev{ari , Numerial Algorithms for Subspae State-Spae System Identiation { An Overview, in Applied and Computational Control, Signals and Ciruits, Vol. 1, Chapter 6, Birkhauser, Boston, 1999. B. De Moor, P. Van Overshee and W. Favoreel K.A. Gallivan, S. Thirumalai, P. Van Dooren and V. Vermaut, High performane algorithms for Toeplitz and blok Toeplitz matries, Linear Algebra Appl. 241-243 (1996), pp. 343-388. [6℄ [7℄ [8℄ [9℄ T. Kailath, S. Kung and M. Morf, Displaement ranks of matries and linear equations, J. Math. Anal. Appl., 68 (1979), pp. 395{407. , Displaement struture: theory and appliations, SIAM Review, 37 (1995), pp. 297{386. T. Kailath and A.H. Sayed T. Kailath, Displaement struture and array algorithms, in Fast Reliable Algorithms for Matries with Struture, T. Kailath and A. H. Sayed, Ed., SIAM, Philadelphia, 1999. , Onand o-line identiation of linear state spae models, Int. J. Control, 49 (1989), pp. 219-232. M. Moonen, B. De Moor, L. Vandenberghe and J. Vandewalle [10℄ C.M. Rader and A.O. Steihardt, Hyperboli Householder transforms, SIAM J. Matrix Anal. Appl., 9 (1988), pp. 269{290. [11℄ V. Sima, Subspae{based algorithms for multivariable system identiation, Studies in Informatis and Control, Vol.5, No.4 (1996), pp. 335{344. Sima, Cholesky or QR fatorization for data ompression in subspae{based identiation?, Proeedings Seond NICONET Workshop, Paris{Versailles, [12℄ V. De. 3 (1999), pp.75{80. 14 Stewart and P. Van Dooren, Stability issues in the fatorization of strutured matries, SIAM J. Matrix Anal. Appl., Vol. 18 (1997), pp. 104{118. [13℄ M. Subspae Identiation for Linear Systems: Theory, Implementation, Appliations, Kluwer Aademi Publishers, [14℄ P. Van Overshee and B. De Moor, Dordreht, 1996. Van Overshee and B. De Moor, N4SID: Two subspae algorithms for the identiation of ombined deterministi-stohasti systems, Automatia, [15℄ P. 30(1) (1994), pp. 75{93. [16℄ M.Verhaegen, Subspae model identiation. Part 3: Analysis of the ordinary output{error state{spae model identiation algorithm, Int. J. Control, 58(3) (1993), pp. 555{586. 15

1/--страниц