- Email: [email protected]
Approximation of multivariable linear systems with impulse response and autocorrelation sequences
0005-1098/83/030265-13503.00/0 PergamonPressLtd. 1983InternationalFederationof AutomaticControl
Automatica, Vol. 19, No. 3, pp. 265-277, 1983 Printedin Great Britain.
Approximation of Multivariable Linear Systems with Impulse Response and Autocorrelation Sequences* YUJIRO INOUYE? A least-squares approximation using afinite data of impulse response and autocorrelation sequences of a multivariable system provides a stable approximate system and it can be computed by an efficient recursive aloorithm. KeyWords---Modelling multivariable systems; correlation method; least-squares approximation; interpolation problem; fast recursive algorithm; sampled data systems. Abstract--This paper considers the construction of approximants of multi-input-multi-output, discrete-time linear systems from the finite data of the impulse response and autocorrelation sequences. In the approximation ofa multivariable linear system, it is common practice to use a finite portion of its impulse response sequence. This is formally equivalent to the Pad6 approximation technique, which may produce unstable approximants, even though the original system is stable. Mullis and Roberts proposed a new method, which yields stable approximants, in connection with approximation of digital filters. This is, however, restricted to the single-input-singie-output case. This paper extends their method to the multi-input-multi-output case and shows a fast recursive algorithm to construct stable approximants of linear systems.
method of approximation for the multivariable case, which contains a difficult test by decision methods for the existence of a solution to multivariable polynomial equations. Furthermore, it is practically impossible to implement their method for reasonable systems on a computer. On the other hand, Mullis and Roberts (1976) proposed a new method, which yields stable approximants, in connection with approximation of digital filters. This is, however, restricted to the single-input-single-output (or scalar) case. This paper extend the Mullis-Roberts method to the multivariable case. Since the approximation of multivariable systems can be regarded as the modelling of them, the extended method can be applied to the modelling of multivariable linear systems. This paper also shows a fast recursive algorithm to construct stable approximants of multivariable linear systems, which is closely related to a multivariate version of Levinson's algorithm in autoregression. For the sake of concise explanation of the results shown in the text, several proofs of the lemmas and theorems have been relegated to appendices. The following notation will be used in the present paper. The symbol tr H denotes the trace of a matrix H which is square. For a matrix H, rank H, H a', H* and H + denote, respectively, the rank, the transpose, the conjugate and the pseudo-inverse (that is, the generalized inverse defined by Penrose). IIn II denotes the Euclidean or Schur norm given by IIH II -- v/(tr HH*). For symmetric matrices K and L, K > L means K - L is positive semidefinite, and K > L means K - L is positive definite. For a positive semidefinite matrix L, ~/L denotes the symmetric square root of L. The symbols I and 0 denote, respectively, the identity matrix and the zero matrix.
1. INTRODUCTION IN THE approximation of a multivariable linear system, it is common practice to use a finite portion of its impulse response sequence. This is formally equivalent to the Pad6 approximation technique (Ledwich and Moore, 1976; Rissanen, 1971; Tether, 1970). It has been widely recognized, however, that the technique may produce unstable approximants, even though the original system is stable. Shamash (1974) and Hutton and Friedland (1975) proposed the methods for obtaining stable approximants, which are, however, limited to the case of singleinput and multiple-output, or multiple-input and single-output. The direct extension of their methods to the multivariable case may yield higher order approximants, even if the original system is of low order. Ledwich and Moore (1976) presented a * Received 13 November 1981; revised 24 May, 1982; revised 24 August 1982. The original version of this paper was presented at the 8th IFAC Congress on Control Science and Technology for the Progress of Society which was held in Kyoto, Japan during August 1981. The published proceedings of this IFAC meeting may be ordered from Pergamon Press Ltd, Headington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication in revised form by associate editor B. Francis under the direction of editor B. D. O. Anderson. t Department of Control Engineering Faculty of Engineering Science, Osaka University, Toyonaka, Osaka 60, Japan. 265
266
YUJIRO INOUYE
2. LEAST-SQUARES APPROXIMATION Consider a linear discrete system with q inputs and r outputs represented by its impulse response matrix sequence {Hk} ~ satisfying the following condition {Ho, H1, H2 . . . . }ff12
(1)
where 12 is the square summable space defined as follows:
system H(z) is represented by the following matrix fraction
Igl(z) = {I + Alz -1 +.., + A,z-"} -1 x{Qo+Qlz
Ye = - - A l y t - 1
k=0
(2)
which is square integrable on the unit circle Izl = 1. Let {ut} be an uncorrelated (wide-sense) stationary q-vector sequence (or a white random process) that is zero-mean and unit-variance, and let { y,} be the output r-vector sequence corresponding to the input {u,}, that is
Yt = ~ HkU,-k" k=O
(3)
The autocorrelation sequence {R/t }_~oo for the impulse response sequence {Hk} is defined as follows:
i=0
33t = -- AlYt-1 -- AEyt-2
-oo
"
'
"
-
-
A,,yt-,, (7)
The least-squares approximation problem considered in this paper is to find the coefficient matrices Ai and Qi that minimize E{ I[Yt -- 33t[[2} • The substitution of these obtained values Ai and Qi into (6) yields an approximant/4(z). Let et = Yt - 33,• It is notationally convenient in the following to define Ao = I and Hk = 0 for k < 0. Using (3) and (7) gives
k
j
. jHk-j-
(4)
where E{" } denotes expectation, C is the unit circle Izl = 1, and H/t = 0 for k < 0. It is clear from the definition (4) that R_/t = R~. The conditions of {u,} and the relation (3) give
H k=E{yt+kuTt},
Q,.ut-,..
+ Qout + Qxut- 1 + ... + Q,,ut-,,.
et =
-~
"'" A , , y t - .
On the other hand, if we made a particular guess for the Ai and Q~, we can predict the value 33t of the present output Yt from the values {Yt-1 . . . . . Yt-.; ut. . . . . ut-m} based the above relation as follows:
R k = E{yt+kyTt} = ~ H,+~H~
= ~j~j H(z)H(z) *?-ldz,
-- A 2 Y t - 2
+ Qout + QlUt-1 + ...
The z-transfer function matrix is defined by
(6)
If the approximant/4(z) is identical to the original system H(z), the inputs u, and outputs Yt are then related by a difference equation
12 = {{no, H i . . . . }l ~ I[nkl[ 2 < 00}. k=O
H(z) = ~ Hkz-/t
1+ ... + Qmz-,.}.
+
k Ut-k
k=m+l
j
)
AjHk-j ut-k.
(8)
Therefore, the following is obtained
E{l'et']2} = trlk~o(j~=oAjHk-~--Qk )
(5)
which implies H k is the cross correlation for timeshift k between the output {y,} and the input {u,}. In the approximation of a linear system {Hk}, the Pad6 approximation technique commonly uses finite matrices {Ho, H1 . . . . . H.,} of the impulse response sequence {Hk}~ as data. Throughout the paper, two kinds of data are given at the same time: the first ones are finite matrices {Ho. . . . . Hm} of the impulse response sequence {Hk} ~ and the second ones are finite matrices {Ro . . . . . R,} of the autocorrelation sequence {Rk}~-~. Assume an approximant H(z) of an original
×
Y
A~Hk_j
)wl •
(9)
Minimizing the above function (9) with respect to {Qo, Q1. . . . . e.,} yields
Qk= ~ AiHk-i j=O
for
O-
(10)
267
Approximation of multivariable linear systems
Definition. For the data {H o. . . . . H,~; Ro ..... Rn}, define a symmetric ( n + 1 ) x ( n + 1) [that is, (n + 1)r x (n + 1)r matrix] K(m, n) with the (i, j)th block element [K(m, n)]~j as follows:
autocorrelation sequence, K(m, n) is positive semidefinite. (2) Minimizing (9) with respect to {Qo. . . . . Qm} satisfies
m-j
[K(m, n)]ij = R j - i - ~ Hk+j-iH~.
(11)
Qk = £ AjHk-j for 0_< k_< m.
(14a)
j=o
k=O
where i, j = 0, 1, ..., n. When {Rk} is the autocorrelation sequence for the impulse response sequence {H~}, (11) is identical with
(3) If K(m, n) is positive semidefinite, then the minimization of (13) with respect to {AI . . . . , An} satisfies
[I, AI ..... An] K(m, n) = Fm.n[l, 0 . . . . . 03 (14b)
(12)
K(m, n)=
E
H'~-I
Lh,-,.J Lh,-,J which implies that K(m, n) is positive semidefinite. By means of (12), the substitution of (10) into (9) becomes
E{lle, ll2} =
tr [A o, A1 . . . . . An]
x K(m, n)[Ao, AI . . . . . An] T. (13) The minimization E{l]e, II2} with respect to {A1, A2 . . . . . An} requires the following lemma, which is verified by the Lagrange multiplier technique.
Lemma 1. Let K be a symmetric positive semidefinite n x n block (or nr x nr) matrix, and let be an n x 1 block (or nr x r) matrix. For the problem which is to minimize tr {XTKX} subject to OTX = I, a solution X exists and satisfies J(TK = ro T
where F is an r x r symmetric matrix uniquely determined by K, and satisfies
Furthermore, the solution ,Y satisfies
XTKX > XTK.~ for any X subject to OXX = I. Using the above arguments, the characterization of the least-squares approximation is obtained as follows.
Lemma 2. For the data {Ho . . . . . Hm; Ro . . . . . R~}, let K(ra, n) be the matrix defined by (11). (1) When the data {Ho . . . . , Hm; Ro, ..., Rn} are taken from an impulse response sequence and its
where Fm.n is a symmetric r x r matrix uniquely determined by the data. I shall call the problem of solving (14) the leastsquares approximation problem for the data {Ho, .... H~; Ro . . . . . Rn} with positive semidefiniteness of K(m, n). Some remarks will now be given here (14). Equation (14b) corresponds to the Yule-Walker equations in autoregression. All the properties derived in the sequel are based on (14). When matrix K(m, n) is nonsingular, the work of Mullis and Roberts in the scalar case could be extended to the multivariable case without essential difficulty. The case where K(m, n) is singular is rare in the scalar case but often occurs in the multivariable case. Equation (14b) may possess two solutions in this case. This fact should be taken into account in the singular case. Thus the questions are whether one can construct a stable approximant /~(z) which is uniquely determined from given data, and whether one can obtain a recursive algorithm for solving (14b) even when K (m, n) is singular. I shall show that these questions have been answered satisfactorily in this work. In the computation of (14b), however, the multiplications usually required are about proportional to n a, and the storage requirements are about proportional to n 2. Therefore, these become very large with increase of size n. Further reduction in computation time and storage is possible in solving (14b) because of the special form of K(m, n). A recursive algorithm in Section 3 will be obtained which is closely related to a multivariate version of Levinson's algorithm in autoregression. By substituting a solution {A 1..... An; Qo..... Qm} of(14) into (6), an approximant H(z) of the original system H(z) is obtained. The z-transfer function H(z) has an impulse response sequence {/tk} satisfying /~k= f 0
for k < 0
Qk -- £ AjI?-Ik-j j=l
for k _>-0.
(15)
268
YUJIRO INOUYE
where Qj = 0 for j > m. The approximant H(z) may be realized as a difference equation of order one in the following form
can be computed from the realization (16) as follows: Hk = ~gk for k < m - n ( CA k - l - ( m "~B f o r k > m - n + l
(16a)
-¢ct+1 = AYct + Bu,_(m_n)
ra--n
1~k = C A k I ( C T + ~ ISIi+kl?IT.
m--n
Ct = C2, + ~ HkU,-k
(18) (19)
i=0
(16b)
k=0
where ~, is an n-block vector (or nr-vector) defined by
-Pt+.-I- k=l~ H~-lU,÷.-k-"
Let / ~ ( m - 1, n - 1) be the matrix defined by (11) for the data {go,..., H,,- 1;/~o .... ,/~,- 1}. Equation (16c) and the same relation as (5) for g(z) then yield the fact that/~ is identical t o / ~ ( m - 1, n - 1). The following lemma, which is concerned with a Lyapunov equation, will be required later.
m
(16c)
Xt ~---
k=i
k=n
and the triplet {A, B, C} is given by
Lemma 3 (Appendix A). Let A and B be n x n and n x r matrices, respectively and let H ( z ) = (zI -- A ) - lB. (1) H(z) is stable if and only if there exist a symmetric positive semidefinite matrix K satisfying the Lyapunov equation K = A K A T + B B T.
(20)
-A
A=
I
(2) If all the eigenvalues of A lie in the open unit disk, then there exists at most one solution of (20). (3) Suppose H(z) is stable• Let
/
0
x
0
0
.•
. . . .
I
I( = ~ AkB(AkB) r.
(21)
k=p
.Bin - 1
B = L ]'Ira- n + 1
C =
[0 ... 0 1].
(16d)
In the case m > n, the above realization (16) of H(z) does not contain in the state space those states of pure delays which are serially connected to the input. Anoter realization can be presented which possesses all these states in its state space, but the realization (16) simplifies the arguments in the following. The stability of the approximant H(z) will be shown in Theorem 1. If the approximant g(z) is stable, then the state vector has a covariance matrix /~ satisfying /~ = E{:~t:~ T} = ~ AkB(AkB) T
(17a)
k=O
and the following Lyapunov equation I( = A I ( A T + B B T.
(17b)
Using the definition of the input sequence {Uk}, the impulse response and autocorrelation sequences
The matrix /( is then a solution of (20) which is symmetric and positive semidefinite. Furthermore, every other solution K of(20), that is symmetric and positive semidefinite, is ordered by K =>/(, namely, /~ is the minimal solution of (20) that is symmetric and positive definite. Some remarks will be given here on Lemma 3. The connections between stability and Lyapunov equations are discussed by Barnett (1971). Part 1 of Lemma 3 is a slight extension of Theorem 4.2 given by Barnett• Part 2 of Lemma 3 is similar to Barnett's Theorem 4.3. In the remainder of this section, the properties of the least-squares approximation problem will be given. The first important one is concerned with stability of the approximant H(z). Theorem 1. Assume the matrix K(m, n) defined by (11) is positive semidefinite. Then the approximant H(z) obtained by (6) and (14) is stable. Proof. Equation (14b) is equivalent to the following pair of equations K(m-
1, n - 1) = A K ( m -
1, n - 1)A x + BB r + BFrn.nBT
(22)
Approximation of multivariable linear systems
k=O
Ak[K(m,
n)]k,. = 0
H~=Hk
(23)
for
269 k = 0 , 1 . . . . . m.
(25)
where A and B are the matrices given by (16d), and /~ is the n x 1 block matrix given by
(2) The double sequence { F u } , whose element is uniquely determined by (14b) for (m, n), satisfies the following:
/~ = [I, 0 . . . O] T.
F..o = Ro - (HoH~ + . . . + H . H I )
(24)
This statement can be derived as follows. Equation (23) follows directly from equating the last block elements in both sides of (14b). From (11), one obtains [ K ( m - 1, n - 1)]i,j = [K(m, n)]i+lj+,
for i,j = 0 . . . . . n--1. This and the block companion form of A in (16d), together with (14b), yield [ A K ( m - 1, n - 1)]i a = [K(m, n)]ij+ 1
(26a)
F~ + ,,. < F.,.
(26b)
F.,. +, < F.,.
(26c)
rank Fro,. = rank K (m, n) - rank K (m - 1, n - 1), form>0, n>0.
(26d)
(3) Assume the data {Ho ..... Hm; Ro ..... R.} are taken from the impulse response and autocorrelation sequences of an original system H(z). Then the rank Fro,. is equal to r - p if and only if p is the largest integer such that for a p x r matrix Vofrank p and some matrices A1 ..... A., Qo..... Q., it holds that
for i,j = 0 . . . . . n - 1 V{y, + A l y , - x + ... + A,y,_,
[ A K ( m - 1, n - 1)AT]ij =
-
[K(m, n)]o,o - Fro,. for (i,j) = (0, O) [ K (m, n) ]ij otherwise. Also there holds from (11) [K(m, n)]ij = [K(m-
1, n -
- H . _ i H . _T j
1)]ij for i,j = 0 . . . . . n - 1 .
Then we obtain (22) from the last two equations. Conversely, it can be shown by an elementary calculation that equations (22) and (23) imply (14b). By (22) and the part 1 of Lemma 3, the z-transfer function ( z I - A)-*[B,/~x/Fm,.] is stable, which implies the stability of/-) (z) from the realization (16). This completes the proof of Theorem 1. A few remarks are given here. Note that the realization (6) may have unstable mode when it is not controllable (see the proof of Lemma 3 in Appendix A). This is not an issue in the scalar case, since the realization (6) is controllable if Fro,, # 0 or if Fro,. = 0 and Fro_ 1,.-1 # 0 (see Mullis and Roberts, 1976). Notice that, when K (m, n) is positive definite, the stability of H(z) could be also derived from the well-known minimum phase properties of predictors as in the work of Morf and co-workers (1977b). Theorem 2 (Appendix A). Presume the matrix K(m, n) is positive semidefinite. (1) The first m + 1 matrices of the impulse response sequence of the approximant H(z) are identical with those of the given data, that is
Qou, - . . . Q,,ut-,,} = o,
a.e.
for all t.
(27)
As a special case, the approximant/~(z) is identical with the original system H (z) if and only if F,,., = 0. Note that there may exist two different solutions of (14b) when the matrix K(m - 1, n - 1) is singular. The following theorem shows that these different solutions of (14b) produce the identical z-transfer function B(z). Theorem 3 (Appendix A). Presume the matrix K(m, n) is positive semidefinite. Then the z-transfer function ~(z) obtained by (6) and (14) is the same for all solutions of (14), that is, /~(z) is uniquely determined by the data {Ho . . . . . Hm; Ro . . . . . R.}. The approximant H(z) will now be examined to see how well matches the given data {Ho . . . . . Hm; Ro ..... R.}. From the part 1 of Theorem 2, the first m + 1 data are matched exactly. Some preparations will be made for the evaluation of differences between the data {Ro . . . . . R.} and the first n + 1 matrices of the autocorrelation sequence of the approximant H(z). Let
17(z) = A(z)- *z-m~/F.,.
(28a)
A(z) = I + AlZ- 1 + . . . + A.z-".
(28b)
where
We see from (22) and the part 1 of Lemma 3 that the system/-i7 (z) is stable. It will be said to be residual corresponding to the approximant/t(z). Let {/-iTk}
270
YUJIRO INOUYE
and {/~k} be respectively the impulse response and autocorrelation sequences of H(z). Note that ~k = 0 for k < m and ~,~ = ~/F,,,. Let ~(n - 1) be a symmetric n x n block matrix with the (i, j)th block element being/~j-i- Then iff(n - 1) is a block Toeplitz matrix whose first row block is [/~o, iffx. . . . . /~n-l]. It can be seen from the definitions that/~(n - 1) is identical to/~(m - 1, n - 1) defined by (11) for the data {0, ..., 0; /~o. . . . . /~n-1}. By means of the realization of H(z) [defined by (16) for the system H(z)], the part 3 of Lemma 3 and (17a), one obtains R ( n - 1) = / ~ ( m -
1, n - 1) = ~ AkBFm.nBT(AT) k. k=O
(29) The following lemma is also required. Lemma 4 (Appendix B). Let {x,} be a purely nondeterministic (wide-sense) stationary n-vector random process, and let K o = E { x , x [ } and K1 = E{xt+lx[} If the covariance Ko satisfies a Lyapunov equation
K o = AKo AT + BB r
m = n of the problem of computing a solution to (14b) may be handled by deriving an equivalent problem for which m' = n' = n. This technique was given by Mullis and Roberts (1976) in the scalar case and is given as follows. Denote the sequence {Ho, H~, H 2. . . . } by 9f ~, and define the operators O"L and aR by O'L~ = {HI, H2, H3,...} aR~,~ =
Let {Rk} be the autocorrelation sequence for 9f~, which is given by (4). For any positive integer k, the autocorrelation sequence for (as)kgf is the same as that of ~¢ itself. However, the autocorrelation sequence for (oL)k~ is given by k-1 j=0
Lemma 5 (Mullis and Roberts, 1976). For any 9f', let K(m, n) be defined by (11). Let
f a[m-n)gf~, if m>n
(30)
9f" =
and a condition
~,
if m = n
(33)
a~-m)gf ~, if m < n .
AKo=K 1
(31)
then Ko is the minimal solution of (30) that is symmetric and positive semidefinite. By means of these preparations, the evaluation of differences mentioned above is shown as follows. Theorem 4 (Appendix B). Presume the matrix K(m, n) is positive semidefinite. Let ~(z) be the approximant obtained by (6) and (14),/7(z) be the residual system defined by (28), and {/~k} and {/~k} be their respective autocorrelation sequences. (1) When the data {Ho . . . . . Hm; Ro . . . . . Rn} are taken from the impulse response and autocorrelation sequences of an original system H(z)
Rk--l~k=Rk
{0, Ho, H1 . . . . }.
for
k = 0 , 1 . . . . . n.
Define K'(n, n) for the sequence ~ ' by (11). Then K'(n, n) = K(m, n). The proof is made by an elementary calculation. We may compute a solution to the (m, n) problem by first deriving the data {H~ . . . . . H'n; R~ . . . . . R'n} and then applying the algorithm in the case m = n to get {AI, A2 . . . . . An}. {Qo, Q~ . . . . . Q=} is computed from (14a) by using the original data {no, ..., H,}. Now a recursive algorithm for computing a solution to (14b) for the case m = n is described. For notational simplicity, let K ( n ) = K ( n , n) and rn = r....
Algorithm. The data {Ho . . . . . Hu; Ro . . . . . Ru} are given. At the completion of step n, the quantities
(32)
A~(n), Bk(n), Ck(n): (0 <~k <%n), F,, f~, (2) Let A be the companion matrix (16d) for a solution {A 1. . . . . An} of (14b). If the eigenvalues of A all lie in the open unit disk (this assumption is attained if Fm,n > 0), the above relations (32) also hold.
have been calculated. To obtain step n + 1, the following calculations are carried out
A. = k=O ~ Ak(n)Rn+l-k
(34a)
An = - ~ Ak(n)Hn+ 1-k
(34b)
3. R E C U R S I V E A L G O R I T H M
This section will provide a re,curslve algorithm for computing a solution {A1 . . . . . An} to the equation (14b) for the case m = n. In the same way as the single input-output (or scalar) case, the case
k=O
O. = A.Bo(n)+A.Co(n )
(34c)
Approximation of multivariable linear systems (I)n = A.Co(n)T + Anf~n rn+,
(34d) (34e)
= rn - anO T - A.e T
fln + I -_ f~. + ~ . TF+. + I ~
(34t)
Ak(n + 1) = A~(n) -- A.B.+ t -k(n) -- A.C.+ 1-k(n) ] B~(n+ l ) - B k ( n ) - - ® . TF . ++ l A. + l - k ( n + l ) -
~
-
Ck(n + 1) = Ck(n) - ~ . TF . ++ 1 A . + , -k(n + 1) 0~
(34g)
where A.+,(n) = B.+,(n) = 0 and C,+t(n) = 0. As initial values for n --- 0, set A0(0) = I, r o = R o - H o HT, Bo(0) = r ~ Co(0 )
= H oT F o+, ~o = I +HToF~Ho .
As far as the matrices {Fo, F, . . . . . FM} are positive semidefinite (see Theorem 5 below), this procedure is carried out until step M is obtained. When K(M) is positive semidefine, this provides a solution {AI(M) . . . . . AM(M)} of (14b) for m = n --- M, and the approximant H(z) is given by (M
a(z)= tk=o 2
)-IM/M )
\
2 ! 2 A,(U)U-A/ k=Ok/=O (35)
If F. = 0 for some step n, the approximant ~(z) for step n is identical with the original system. If F~ is not positive semidefinite for step n, then K(n) is not positive semidefinite (see Theorem 5 below), which implies there is no system which matches the data {Ho . . . . . H.; Ro . . . . . R.}, that is, the data are never identical with those taken from the impulse response and autocorrelation sequences for any system with its impulse response sequence in space !s (see Theorem 6). The proof of validity of the algorithm is found in Appendix C. If one sets H~ = 0 for 0 ~< i ~
271
pseudo-inverse (Golub and Kahan, 1965; Dongarra and co-workers, 1979). The matrices A. and A. are the errors in the predicted values of R.+I and Hn+,, respectively, based on the nth approximation. If these both vanish, then the n + 1st approximation is equivalent to the nth. The Levinson algorithm in autoregression was also extended to non-Toeplitz cases (however closeto-Toeplitz cases) by Morf and co-workers (1977a) and Friedlander and co-workers (1979), where they employed the displacement rank properties. As shown below, the above algorithm can be used to examine whether the matrix K(m, n) for some data {Ho . . . . . Hm; Ro . . . . . R.} is positive semidefinite or not. Theorem 5 (Appendix C). Let K(M) be the matrix defined by (11) for some data {H o. . . . . HM; Ro. . . . . RM}, and {Ao(j) . . . . . A~(j); F j} be the quantities generated from the above algorithm (34) at step j. Then, K(M) is positive semidefinite if and only if Fo . . . . . F u are all positive semidefinite and {Ao(j) . . . . . Aft.j); F~} satisfies (14b) for m = n = j , where 0 ~
Yt = ~ HkUt-k + v,
(36)
k=0
where {Hk} •12,{ut} is the input sequence possessing the properties explained in Section 2, and {v,} is a stationary noise random process uncorrelated to {u,} with mean-zero. When some data {Ho . . . . . Hm; Ro . . . . . R.} are given, the problem considered here is to construct a system with noise (36) such that
E{yt+ku T } = H k for 0~
for
O<~k<~n.
(37a) (37b)
This is called the interpolation problem for the data {Ho, ..., Hm; R o, ..., R,}. The data is said to be consistent if and only if there exists a solution to *The least-squares approximation problem can also be formulatedfor the system(36) instead of (3), and obtain that it is characterizedby (14), where the data are computed by (37).
272
YUJIRO INOUYE
the problem. When there exists a solution to the problem, it is required that a system is constructed with noise such that the covariance E{v:~} of noise is as small as possible. When the covariance E{vt vXt} may be reduced to zero, the system obtained is said to match the data exactly. If a system H(z) can be represented in the form (6), it will be said to be of order at most (m, n). A solvability condition of the interpolation problem and a procedure for obtaining a system with noise that matches the data will be given. Theorem 6 (Appendix D). The data {H o. . . . . H~,; R0 . . . . . R.} for the interpolation problem is consistent if and only if the matrix K(m, n) defined by (11) for the data is positive semidefinite. Given the data {H 0. . . . . Hm; Ro .... , R.} for which K(m, n) >=O, there are two cases where Fro,. = 0 or Fro,. # 0. In the former case, there almost always exists a system/4 (z) which matches exactly the data. This is shown below. In the latter case the problem of artificially adding H= + 1 (or R. + 1) in a way that K(m + 1, n) > 0 and F r o + l , n = 0 [or K(m,n + 1) =>0 and Fro,. +1 = O] will be considered. Theorem 7 (Appendix D). Assume the positive semidefiniteness of K(m, n)for the data {Ho. . . . . Hm; R o. . . . . R.}. (1) There exists precisely one system H (z) of order at most (m, n) which matches exactly the data if and only if Fro,. = 0 [or equivalently rank K(m, n) < nr] (38a)
K ( m - 1, n - 1) = ~ AkBBT(AT) k
(38b)
k=O
where A is the companion matrix given by (16d) for a solution {At . . . . . A,} of the (m, n) least-squares approximation problem, and B is defined by B T = [-H T' Hm_ T t . . . . . H sT_ . + t ] .
(38c)
(2) When F r o . . # 0 [or equivalently rank K(m, n) > nr], there are at least two values of Hm+ 1 for which K(m + 1, n) >- 0 and Fm + t,. = 0 in the case rank Ro < q, and there are at least two values of R.+I for which K(m, n + 1) =>0 and Fro,.+ 1 = 0. Two examples are presented to illustrate the above results. These become counterexamples to the statements b and c of Proposition 1 provided by Mullis and Roberts (1976). Example 1. Consider the data { 1, 0; 2, 1} in the scalar case. It follows from the definition (11) that K(1, 1)=
LI 1]>0 1 =
which implies rank K(1, 1) = 1, and so F 1,1 = 0. The solution of(14b) for these data is given by A ~ = - 1. We obtain B = 0 in (38c), which implies the value of the right hand side of (38b) vanishes. Therefore, the condition (38b) does not hold for the data { 1, 0; 2, 1}. Part 1 of Theorem 7 shows that there is no system that matches exactly the data. This indicates that the statement b of Proposition 1 presented by Mullis and Roberts (1976) is not correct. It states that if K(m, n) is positive semidefinite but det K(m, n) = O, then there is precisely one system which matches exactly the data. The approximant/4(z) = 1 can be obtained for the data {1, 0; 2, 1}. By a simple calculation, it is seen that the following system with noise Vo matches the data Yt = Ut + V0 where mean E{vo} = 0 and variance E{v 2} = 1. Example 2. Consider the data { 1, 0; 2} in the scalar case. K(1, 0) becomes 1, and so the two values of H 2 are equal to + 1 for which F2. o = 0. It can be seen that H(z) = 1 _+ z -2 for the data {1, 0, + 1; 2}. The approximant H(z) matches exactly the data {1, 0, + 1; 2}. On the other hand, the two values ofR 1 become _+1 for which F 1A = 0. It follows that the condition (38b) does not hold for the data {1, 0; 2, +1} and thus there is no system that matches exactly the data { 1, 0; 2, 1}, or { 1, 0; 2, - 1}. This indicates that the statement c of Proposition 1 provided by Mullis and Roberts (1976) is not true. The approximants for the data {1, 0; 2, + 1} are same and are given by H(z) = 1. By a simple calculation, it follows that the following system with noise (_+ 1)tVo matches the data y, = u, + (+ 1)% where mean E{vo} = 0 and variance E{v2o} = 1. 5. CONCLUSIONS In this paper, the least-squares approximation problem was reduced to that of solving a linear (matrix-valued) equation. This equation can be solved recursively by the algorithm, which was closely related to a multivariate version of Levinson's algorithm in autoregression. A procedure for obtaining the solutions of the interpolation problem is presented through the solution of the least-squares approximation problem. A problem left unsolved however is that of finding the system of minimal order that matches exactly the data. Acknowledgements--The author
would like to thank Professor Y. Sakawa for constant encouragement and Dr H. Kimura for useful discussions.
273
Approximation of multivariable linear systems REFERENCES Barnett, S. (1971). Matrices in Control Theory. Von Nonstrand Reinhold. Boullion, T. L. and P. L. Odell (1971). Generalized Inverse Matrices. John Wiley. Dongarra, J. J., C. B. Moler, J. R. Bunch and G. W. Stewart (1979). L I N P A K Users" Guide. SIAM Press. Friedlander, B., M. Moff, T. Kailath and L. Ljung (1979). New inversion formulas for matrices classified in terms of their distance from Toeplitz matrices. Linear Algebra & Applic., 27, 31. Golub, G. and W. Kahan (1965). Calculating the singular values and pseudo-inverse of a matrix. J. S I A M Numer. Anal., Ser. B, 2, 205. Hutton, M. F. and 13. Friedland (1975). Routh approximations for reducing order of linear, time-invariant systems. IEEE Trans Aut. Control, AC-20, 329. Inouye, Y. (1981). Approximation of multivariabte linear systems with impulse response and autocorrelation sequences. Technical report no. 81-05, Department of Control Engineering, Osaka University, Japan. Ledwich, G. and J. B. Moore (1976). Minimal stable partial realization. Automatica, 12, 497. Markel, J. D. and A. H. Gray, Jr (1973). On autocorrelation equations as applied to speech analysis. IEEE Trans Audio Electroacoust., AE-21, 69. Morf, M., B. Dickinson, T. Kailath and A. Vieira (1977a). Efficient solution of covariance equations for l!near prediction. IEEE Trans Acoust., Speech, Signal Processing, ASSP-25, 429. Morf, M., D. T. Lee, J. R. Nickolls and A. Vieira (1977b). A classification of algorithm for ARMA models and ladder realization. Proc. 1977 IEEE Conf. on Acoust., Speech, Signal Processing. Mullis, C. T. and R. A. Roberts (1976). The use of second-order information in the approximation of discrete-time linear systems. IEEE Trans Acoust., Speech, Signal Processing, ASSP-24, 226. Rissanen, J. (1971). Recursive identification of linear system. S I A M J. Control, 9, 420. Rosenbrock, H. H. (1970). State-space and Multivariable Theory. Nelson. Shamash, Y. (1974). Stable reduced-order models using Pad& type approximations. IEEE Trans Aut. Control, AC-12, 615. Thether, A. J. (1970). Construction of minimal linear state variable models from finite input-output data. IEEE Trans Aut. Control, AC-15, 427. Whittle, P. (1963). On the fitting of multivariate autoregressions, and the approximate factorization of a spectral density matrix. Biometrika, 50, 129.
(1 - 12[2)v*Kv = IIv*e II2 >~0.
(A5)
If 2 is an unstable eigenvalue of A, that is, 12l > 1, then (A5) implies v*B = O. Therefore v*[a - 21, B] = 0
(A6)
rank [A - 21, B] < n.
(A7)
which yields
This fact means that 2 is an input-decoupling zero, namely, the mode associated with 2 is uncontrollable (Rosenbroek, 1970, p. 64). Hence this unstable mode is absent in the output behavior of H (z). Proof of(2) of Lemma 3 The uniqueness of the solution of (20) comes from the fact that if all the eigenvalues of A lie in the open unit disk, then the matrix equation K = A K A r has the unique trivial solution K=0. Proof of(3) of Lemma 3 The former part obviously follows from (AI) and (A2) above. The latter part comes as follows. Let
A(n) =
~
AkB(A~B) r.
(A8)
kfn+l
Then from (A8), (AI) becomes
g = ~ AkB(A~B)T + A(n).
(A9)
k=O
Since the sequence {A(n)} is monotonically decreasing a n d / ( exists, it follows lim A(n) = O.
(AI0)
n~ov
From (20), one obtains K = A K A T + BB T A K A • = A2K(A2) T + AB(AB) T A"K(A~)T = A" +l K ( A ~+t)T + A.B(AnB)T.
(All)
Hence
APPENDIX A
n
Proof of(I) of Lemma 3 If H(z) is stable, then there exists the limit ~ A~B(AkB)r as shown in (17a), which is denoted by k=o I( = ~ AkB(AkB) x.
(A1)
K = An+IK(A~+I)T+ ~ AkB(A~B)T.
(AI2)
kffiO
Combining (A9) and (A12) gives K - I?, = A n+ IK(An+ 1)r _ A(n).
(AI3)
k=O
Since K is positive semidefinite, (AI3) implies
Then (A1) becomes
g - / ~ _>_-A(n).
Letting n--* ~ in the above inequality together with (A10) yields K>=K.
kffi
Conversely, assume that (20) holds. Let 2 be an eigenvalue of A and let v* be a left eigenvector associated with 4. Then v*A = 2v*,
v* :~ 0.
(A3)
From (20) and (A3), it follows v*Kv = o*AKArv + v*BBrv = 1212v*Kv + v*BBrv
which means
(A14)
AJ'B(AtB) T A r = AI~A T + BB r. (A2)
l( = BB T + A
(A4)
Proof of Theorem 2
The relation (25) follows from (14a) and (15). Equation (26a) comes directly from the definition (11) of K(m, n) and the fact that Fm.o = K(m, 0). The relations (26b) and (26c) follow from Lemma 1, the definition (11) of K(m, n), and the fact that the solution to the (m, n) problem satisfies the constraints for the (m + 1, n) and (m, n + I) problems. Equation (26d) comes as follows: let {A t, A 2. . . . . A,; Fro.n} be a solution of (14b). Let T
274
YUJIRO INOUYE
be the (n + 1) x (n + 1) block matrix, whose diagonal elements are unity, defined by
T=[; A,.;.A.] lblOCknblocks.
The above two equations give
/~;-,]*
I
F/t',-1]
(A15)
..... / , J
[D1,
D,]
!
Then from (14b) and the definition (11) of K ( m , n) TK(m, n)T r =
It
K(m-
0
]
Since T is nonsingular, the above equation implies (26d). From (8) and (10), the following can be obtained in an argument similar to the derivation of (13) E { e , eT} = [A o. . . . . A . ] K ( m ,
, [ D 1. . . . .
Dn] T = 0
for t > m +
(A16)
1) "
1, n -
!
/~;-.J
n)[A o. . . . . A.] r.
1
which implies
• Dj/~',_j = 0
for t _->m + 1.
(A25)
./=1
(A17)
Assume/t~, = / t ~ for k = 0, 1. . . . . t - 1, where t ->_m + 1. Then from (15) and (A25), it follows that
This and (14b) give E { e , e~} = F.,..
(A18)
On the other hand, (7) yields
%
H t - H
""t =
_
~, (A'jI$1't_jA']H;'_j) --
i=i
e, = Yt - Y', = Y, + A I Y , + . . . + A,,yf_,, -- Qou, - . . .
- Q,,,u,_,..
It is clear that rank F.,. < r - p if and only if there exists a p x r matrix of rank p such that VFm.V r = 0. Therefore, rank F.,. < r - p if and only if V E { q e'~}vr = 0, or equivalently Ve, = 0 a.e. for all t, which implies (27). This completes the proof of Theorem 2.
=-
~Dj/t;_j=0
APPENDIX B
{At. . . . . A;} and {A~'. . . . . A~'}be two
solutions of (14b).
L e t D~ = A'~ - A':, where i = 1, 2 . . . . . n. From the two equations
obtained by substituting {A'1. . . . . A~} and {A~'. . . . . A~'} into {A~ . . . . . A,} of (14b), respectively, one obtains [D l .....
(A26)
This means the sequences { ~ } and { / ~ } are identical, which completes the proof of Theorem 3.
Proof of Theorem 3
Let
fort=>m+l.
j=l
D , , ] K ( m - 1, n - 1) = 0.
for
All the properties mentioned in Lemma 4 are not changed by a coordinate transformation, we can assume without loss of generality that the pair {A, B} is written in the following reachable canonical form
(A19)
Let /t'(z) and /t"(z) be respectively the z-transfer functions [obtained from (6) and (14)] corresponding to {A'l . . . . . A',,} and {A'~ . . . . . A'.'}, let {H~} and {/-/~'} be respectively their impulse response sequences, and let ~ ' ( m - 1, n - 1) be the covariance matrix given by (17a) for the system/~'(z). Part 1 of Theorem 2 gives
H;=H~=H;'
P r o o f o f L e m m a 4.
k = 0 , 1 . . . . . m.
(A20)
A=[A1
;:],
B = [ B1]
(A27)
where the pair {A 1, BI } is reachable (or controllable from origin). According to the above partition, let x, =
F
x~ '
[Kp
r 0 = LK2
K022_], K, -- LK**'
r~2J"
(A28)
From (30), (31), (A27), and (A28)
From (22)
Ko22 = A 2 K 022A2r A2 K 2 2 = K122.
K ( m - 1, n - 1) = A ' K ( m - 1, n - 1)A 'r + B ' B 'T + ~ F . , . ~ r.
(A21) where A' and B' are the matrices given by (16d) for the system /t'(z). Part 3 of Lemma 3 gives K ( m - 1, n - 1) > ~ A ' k ( B ' B 'T + BF,.,./~X)(A'r) k
(A29) (A30)
From (A29) and (A30), it follows K 22 = A2(K22) T.
(A31)
Using (A30) and (A31) yields
k=O
E{ (x2+1 - A2x2)(x2+ , -
> ~ A ' k B ' ( A ' k B ' ) "r = K ' ( m -
1, n -
1).
k=O
(A22)
D,,]K'(m-
1, n - 1)[Dr . . . . . D.] r = 0.
(A23)
= Ko - K1 A2T -- A2(K22) T -t- A 2 K 022 A T2 = 0
(A32)
~
I~;-1] fq;--1]T •
',-, L m-.l
x2+ 1 = A2 x2,
a.e. - oo < t < oo.
(A33)
Since {xf} is purely non-deterministic, this means x 2 = 0 a.e. for all t. Thus
On the other hand from (16c) [or 12] one obtains
g ' l m - 1, ~ - 1) =
A2x2)T}
22
which implies
This and (A19) yield [DI .....
22
.
K 2 = o, K o~' (A24)
=
o, K 2 = o.
(A34)
This yields, together with (30) and (A28) K ~ I = A I K o11 AIT + BxB~"
(A35)
Approximation of multivariable linear systems Since the pair {At, Bt} is reachable, by parts 1 and 2 of Lemma 3, it follows
275
Equation (23) which was derived from (14b) gives [K(m, n)]o,. = - ~'. A~[K(m, n)]k,.
K ~ I = ~ A~B~(A~BOT
(A47a)
kffiI
(A36)
kffiO
Similarly one obtains which means, together with (A27), (A28), and (A34) JR(m, n)]o.. = Ko = ~ AkB(AkB) r.
~ Ak[g(m, n)]k.~
(A47b)
kffil
(A37)
kffiO
~. = - ~ AtR.-k. By part 3 of Lemma 3, this implies that Ko is the minimal solution of the Lyapunov equation (30) that is symmetric and positive semidefinite. This completes the proof of Lemma 4. Proof of Theorem 4 Let xt be an n-block vector #oven by
(A47c)
k=l
From the definition (11), the last block columns of K(m, n) are respectively #oven by m-n
[ K ( m , n)]e,. = R . - e - Y, Hk÷.-,H~ for i = 0 . . . . . n
(A48a)
k=O m-n
~ nk-tUt+n-k
YI+n-I-
[/~(m, n)]~,. = R . _ / -
kffil
~ /-~k+._d-irk for i = 0 . . . . . n.
(A48b)
k=0
xt =
Yt+.-i - ~. Hk-eut+.-k
(A38)
Since H e = / l e for 0 -< i -< m, (A46)-(A48) yield
kffit
R. = / ~ . +/~..
Yt
Hk-nUt+s-k kffin
-- ~
Since the input sequence {ut} is stationary and purely nondeterministic, so the output sequence {Yt} and therefore the sequence {xt} defined by (A38) become. Let Ko = E{x, x,r} and KI = {xt÷ 1 xtr} • Then from an elementary calculation, (A38) and (11) give Ko = K(m
-
[K1]ea = [K(m, n)],a+~
1, n -
1)
(A39)
for i,j = 0, 1. . . . . n - 1. (A40)
For the #oven data, (22) and (23) are obtained which are equivalent to (14b), These, together with (A39) and (A40), ensure the relations Ko = AKo Ar + BB T + BFm,n/~T
(A41)
A K o = Kt.
(A42)
By Lemma 4, Ko is the minimal solution of (A41) that is symmetric and positive semidefinite. This yields by the part 3 of Lemma 3 K(m-
1,
n-
This implies the part 1. If the eigenvalues of the companion matrix A all lie in the open unit disk, the relation (A43) also follows from (22) and the part 2 of Lemma 3, which yields the relations (32) in the same argument as above. This completes the proof of Theorem 4. APPENDIX C Proof of validity of the algorithm (34) The following notation will be used. For a matrix K, Im K and Ker K denote, in the same order, the image of K and the kernel of K. K" denotes a reflexive generalized inverse of K, that is, a matrix that satisfies K K ' K = K and K ' K K " = K" (Bouillon and Odell, 1971). Using the definition (11) of K(n), it can be shown that with the proper definitions of the matrices E(n) and D(n)
/ffiO
kffiO
(A43) which implies from (17a) and (29)
I-"'7
V
K(n + l) = /
K(n)
I R-n- 1
1)= kffiO ~. Ak[B,g~(F.,.)][B, B~(F.,,.)]T(A~ = ~ AtBBT(AT~ + ~ AkBF.,.IIT(A~
(A49)
...
x [n.r+l ..... Hor]
=l=E(n) I D(n)r]. L D(n)
K(n)J
(A50)
Let A(n) = I-Ao(n). . . . . A.(n)]
(ASia)
B(n) = [Bn(n). . . . . Bo(n)]
(A51b)
C(n) = [C.(n) . . . . . Co(n)]
(A51c)
K(m - 1, n - 1) = / ( ( m - 1, n - 1) + / ( ( m - 1, n - 1). (A44) where A(n) satisfies (14b) for the case m = n, and F r o m t h e definition (11), the first block rows of K ( m - 1, n - 1) and K ( m - 1, n - 1) are respectively #oven by [K(m-
1, n - 1)]o.1 = Rj-
=-l-j kffiO m-l-j
~
(A52a)
C(n) = [ n l . . . . . nro]K(n) ".
(A52b)
Hk+1HT
for j = 0, 1. . . . . n - 1. (A45a) EK(m-l,n-l)]oa=l~1-
B(n) = [0 .... O, lJK(n)"
~
H~+~HT
Here K(n)" is a reflexive generalized inverse of K(n) recursively given by Lemma C2 below. The following two lemmas are required: /.emma C 1. Suppose K(n + 1) is positive semidefinite. Then
k=0
for j = 0, 1. . . . . n - 1. (A45b)
Im D(n) = Im K(n) Ker {[0 . . . . . 0, A.] + A . [ H .r . . . . . H~]} = Ker K(n)
(A53) (A54)
Since H~ = / ~ l for 0 ~ i ~ m, (A44), (A45), and (29) yield R~ = i~ + ~,j for j = 0, 1. . . . . n - 1.
(A46)
where A. and A. #oven by (34a) and (34b), respectively. Proof: It is easy to show that (A53) follows from (A50) and
276
YUJIRO INOUYE
the positive semidefiniteness of K(n + 1), From (A53) and the symmetry of K(n), it follows
which implies
I Bn(n)T ] Ker D(n) r ~ Ker K(n).
(A55)
:
=
[Adn ) ... An(n)]
LB,in:5
On the other hand, it holds [F., 0 . . . . . 0, A.] + An[H.X+ 1. . . . . HoT] = [A(n), 0]K(n + 1) = [A(n), 0][E(~2)
O(n)r] (A56) K(n) A"
It is easily shown that (A55) and (A56) imply (m4). Lemma C.2. Suppose A(n) and r . satisfy (14h) for the case m = n, and K(n) is positive semidefinite. Then a reflexive
generalized inverse of K(n) that is symmetric is given by K(n)" = A(n)TF~*A(n) +
D ( n - 1)T
.
(A66)
Equation (A59a) therefore becomes, from (A63), (A66), and (34a)(34c) O . = D(n)T B(n) T
A.Bo(n) + A.Co(n) - L A~(n)F.+ 1 -k
=
k=l
[01 0 ? 0
{ [i:l -,-- I -,l = ®. - L Ak(n)Fn+l-k"
K(n - 1)'
(A67)
k=l
By means of (A53) for n - 1, it is easy to show that K(n)" given by (A57) satisfies the definition of the reflexive generalized inverse. From (A50) and the above two lemmas, it is easy to derive the following:
A(n + 1) = [A.(n), 0] - A.[0, B(n)] - A.[0, C(n)]
(A58a)
B(n + 1) = I-0, B(n)] -OrF~++~A(n + 1)
(A58b)
C ( n + l ) = [0, C(n)] - @. "r F.+ + 1A(n + 1)
(A58c)
rn+, = r n - A.~).x - A.~.r
Here use is made of the fact Bo(n) = Bo(n)r. It follows from (A62b), (A50) and (A61) that [G.X ... G~x] = [C.(n) ... C,(n)]r(n - 1)+ Co(n)[R_.... R_ 1] -
O ( n - 1)r
where ~)., ~., ~). and ~ . are given as follows: (A59a)
~n = D(n)rC(n) T - H.+ ~
(A59b)
~,, = - A ( n + 1)[0, .... 0, I] 7
(A59c)
@n = - A ( n + 1)[H T"+~ . . . . . n ~ ] r.
(A59d)
(A68)
which means, together with (A64)
(A58d)
~ . = O(n)r B(n) r
fJ.[H .T. . . H r]
x
• Lc,in:J
= [Al(n ) .., A~(n)]
{?] R1
C ° ( n ) r - L_HI_J
:'1} ~1
. (A69)
Equation (A59b) therefore becomes, from (A63), (A69), (A61), and (34a), (34b), and (34d) n
Here the following equation is used which is derived from (A52) and (A54):
~) n = D(n)rC(n) T
H.+ 1 =
-
AnC0(n)
r
+ A.fl. - k~l Ak(n)Gn+ 1 -~
n
= On -- k~=tAk(n)Gn+ 1 -k.
[ 0 . . . 0 A.] + A.[nr~... H T] - {A.B(n) + A.C(n)}K(n) = { [ 0 . . . 0 A , ] + A n[H r . . . H ro] } {I --
K(n)'K(n) }
= 0.
(A60)
(A70)
Here use is made of the fact fl n = ll~. On the other hand, since
Define
A ( n + 1)K(n + 1) = [F.+I, 0 ... 0] n . = t + c(n)EtC .., HOr] ~.
(A71)
(A61) it follows from (A50), (A58), (A59a), and (A59b) that
Let [ F , + l , 0 ... 03 = [F,, 0 ... 0, A,] + A.[0, H~ ... H~]
[ F ~ . . . FOr] = B(n)K(n)
(A62a)
[Gr~ ... GOr] = C ( n ) r ( n ) - EH~ ... HT].
(A62b)
From the definition of D(n), it follows D(n)T = [ D ( n - 1)r, R , + , ] - H . + , [ H ]
- [A.O.r + A.t~T, A,B(n)K(n) + A.C(n)K(n)] which implies, together with (A62) IF,+ 1, 0 . . . 03 = IF,, 0 . . . 0, An]
... HTo]. (A63)
--T T - EA.e, + a . ¢- T, , AnEF.... For ] + A.[G~ ... GOr]].
(A73)
Therefore, it holds
Equations (14b) and (A50) give
D(n - 1)7 = - [Adn) ... A.(n)]K(n - 1).
(A64)
AnEF~ ...
F ,r ] + A,[G,T ... GI] = 0.
(A74)
Using (A67), (A70), and (A74), one obtains
It follows from (A62), (A50) and (A52) that
A . ( 6 r _ ®T) + An(~x _ @1-) = 0
[F~r
(A72)
FIr] = B(n) L k - . . . . R_,
(A75)
which yields, together with (A58d), Fn+ 1 = F. - A.@T - A.@T.
= [B.(n) ... n l ( n ) ] K ( n - 1)+ Bo(n)[R_,,... R_ 1]
- C o ( n ) r [ H ] . . . HI]
(A65)
(A76)
It follows from (A58a), (34b), (A52a), (A52b), and (A61) that
277
Approximation of multivariable linear systems T A(n + I)[H.+,, .... H ~ ] T = [A(n), 0 ] [ H . T+ I ..... HoT] T
where i, j, = 0, 1..... M. Then, we obtain from (14b)
- A.E 0, B(n)] [H.T+, ..... H ~ ] T -- A.[0, C(n)] [H.T+, ..... HoT] T
T(M)K(M)T(M) T = block diag {Fu, F u - , . . . . . Fo}. (A81)
= A(n)[HT.+, ..... H [ ] T - A.B(n)[H T ..... s T ] T
-- A.C(n)[HT. . . . . . HTo] T = -- A. - A.[0 ... 0 I]K(n)" [HT. ..... H T ] T -
A.C(n)[HT.. . . . .
ST] T
Since T(M) is nonsingular, the positive semidefiniteness of K(M) implies that of all Fo, F1 ..... Fu. Conversely, if {Ao(j).... , A~); F~} satisfies (14b) for m = n = j (where 0 < j ~ M), then (A81) holds. Therefore, if Fo, F1. . . . . F u are all positive semidefinite, then K(M) is also positive semidefinite.
= -- A.[0 ... 0 I]C(n) T - A.{I + C(n)[HT. ..... Ho~] ~} = - A.Co(n) T - A.f2.
which, together with (A59d) and (34d), implies ~. = @.. It follows from (A58a), (A59c) and (34c) that 0 n =
APPENDIX D
(A77)
- - A ( n -]- I)[0 ... 0
~/]
= A,[0, B ( n ) ] [ 0 ...
I] T
Proof of Theorem 6 If the interpolation problem for the data {Ho..... H,; Ro..... R,} is solvable, then there exists a system with noise (36). Since the random sequences {u,} and {v,} are uncorrelated, (36) and the definition (11) of K(m, n) give
+ A.[0, C(,)] [0... o z]~ = A.Bo(n) + A.Co(n) = O..
(A78)
Kim. ,) =
~
.n~_~
~-~
+ ~(,) (AS2)
k>ra+l
Using (A61), (A58c), (A59d), and the fact that ~. = (~., one obtains n.+, -f~. = {C(n + I ) - [0, C(n)]} In.T+, ..... HTo] T = - ~ .^r F . + t+A ( n + ~T
+
^
__
I)[HT+, . . . . . S ~ ] ~ T
+
= (b. F.+ ,~. - ~. F.+ ~ .
(A79)
which becomes (34f). This completes the derivation of the algorithm (34).
-
.Hk_..
u_ilk_..
where ~(n) is the positive semidefinite (n + 1) × (n + 1) block matrix with the (i, j)th block element being E{v,_, v~_~}. This implies that K(ra, n) is positive semidefinite. Conversely, assume K(m, n) is positive semidefinite. Let yf = ~ /t~u,_, + v, vt = ~, ~k~t-k + 6
Proof of Theorem 5 AS far as K(M) is positive semidefinite, {Ao(j). . . . . Aj(j); F j} satisfies (14b) for m = n =j, where O~_j<_M (see the above proof). Let T(M) be the upper triangular (M + 1)x (M + 1) block matrix defined by .... f0, for i > j [ltMLh'j = ~Aj_,(M - j ) ,
(A80) for i < j
(A83a)
k=O
(A83b)
k=0
where {/~,[ and {~,} are the impulse response sequences for H(z) and H(z), respectively, {6,} is a white noise process uncorrelated to {u,}, and {~,} is a deterministic process uncorrelated to {u,} and {at}. Then it can be shown that the system (A83) matches the given data. The details of the proof is found in the work of Inouye (1981). The proof of Theorem 7 is also found in the work of Inouye (1981).
Automatica, Vol. 19, No. 3, pp. 265-277, 1983 Printedin Great Britain.
Approximation of Multivariable Linear Systems with Impulse Response and Autocorrelation Sequences* YUJIRO INOUYE? A least-squares approximation using afinite data of impulse response and autocorrelation sequences of a multivariable system provides a stable approximate system and it can be computed by an efficient recursive aloorithm. KeyWords---Modelling multivariable systems; correlation method; least-squares approximation; interpolation problem; fast recursive algorithm; sampled data systems. Abstract--This paper considers the construction of approximants of multi-input-multi-output, discrete-time linear systems from the finite data of the impulse response and autocorrelation sequences. In the approximation ofa multivariable linear system, it is common practice to use a finite portion of its impulse response sequence. This is formally equivalent to the Pad6 approximation technique, which may produce unstable approximants, even though the original system is stable. Mullis and Roberts proposed a new method, which yields stable approximants, in connection with approximation of digital filters. This is, however, restricted to the single-input-singie-output case. This paper extends their method to the multi-input-multi-output case and shows a fast recursive algorithm to construct stable approximants of linear systems.
method of approximation for the multivariable case, which contains a difficult test by decision methods for the existence of a solution to multivariable polynomial equations. Furthermore, it is practically impossible to implement their method for reasonable systems on a computer. On the other hand, Mullis and Roberts (1976) proposed a new method, which yields stable approximants, in connection with approximation of digital filters. This is, however, restricted to the single-input-single-output (or scalar) case. This paper extend the Mullis-Roberts method to the multivariable case. Since the approximation of multivariable systems can be regarded as the modelling of them, the extended method can be applied to the modelling of multivariable linear systems. This paper also shows a fast recursive algorithm to construct stable approximants of multivariable linear systems, which is closely related to a multivariate version of Levinson's algorithm in autoregression. For the sake of concise explanation of the results shown in the text, several proofs of the lemmas and theorems have been relegated to appendices. The following notation will be used in the present paper. The symbol tr H denotes the trace of a matrix H which is square. For a matrix H, rank H, H a', H* and H + denote, respectively, the rank, the transpose, the conjugate and the pseudo-inverse (that is, the generalized inverse defined by Penrose). IIn II denotes the Euclidean or Schur norm given by IIH II -- v/(tr HH*). For symmetric matrices K and L, K > L means K - L is positive semidefinite, and K > L means K - L is positive definite. For a positive semidefinite matrix L, ~/L denotes the symmetric square root of L. The symbols I and 0 denote, respectively, the identity matrix and the zero matrix.
1. INTRODUCTION IN THE approximation of a multivariable linear system, it is common practice to use a finite portion of its impulse response sequence. This is formally equivalent to the Pad6 approximation technique (Ledwich and Moore, 1976; Rissanen, 1971; Tether, 1970). It has been widely recognized, however, that the technique may produce unstable approximants, even though the original system is stable. Shamash (1974) and Hutton and Friedland (1975) proposed the methods for obtaining stable approximants, which are, however, limited to the case of singleinput and multiple-output, or multiple-input and single-output. The direct extension of their methods to the multivariable case may yield higher order approximants, even if the original system is of low order. Ledwich and Moore (1976) presented a * Received 13 November 1981; revised 24 May, 1982; revised 24 August 1982. The original version of this paper was presented at the 8th IFAC Congress on Control Science and Technology for the Progress of Society which was held in Kyoto, Japan during August 1981. The published proceedings of this IFAC meeting may be ordered from Pergamon Press Ltd, Headington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication in revised form by associate editor B. Francis under the direction of editor B. D. O. Anderson. t Department of Control Engineering Faculty of Engineering Science, Osaka University, Toyonaka, Osaka 60, Japan. 265
266
YUJIRO INOUYE
2. LEAST-SQUARES APPROXIMATION Consider a linear discrete system with q inputs and r outputs represented by its impulse response matrix sequence {Hk} ~ satisfying the following condition {Ho, H1, H2 . . . . }ff12
(1)
where 12 is the square summable space defined as follows:
system H(z) is represented by the following matrix fraction
Igl(z) = {I + Alz -1 +.., + A,z-"} -1 x{Qo+Qlz
Ye = - - A l y t - 1
k=0
(2)
which is square integrable on the unit circle Izl = 1. Let {ut} be an uncorrelated (wide-sense) stationary q-vector sequence (or a white random process) that is zero-mean and unit-variance, and let { y,} be the output r-vector sequence corresponding to the input {u,}, that is
Yt = ~ HkU,-k" k=O
(3)
The autocorrelation sequence {R/t }_~oo for the impulse response sequence {Hk} is defined as follows:
i=0
33t = -- AlYt-1 -- AEyt-2
-oo
"
'
"
-
-
A,,yt-,, (7)
The least-squares approximation problem considered in this paper is to find the coefficient matrices Ai and Qi that minimize E{ I[Yt -- 33t[[2} • The substitution of these obtained values Ai and Qi into (6) yields an approximant/4(z). Let et = Yt - 33,• It is notationally convenient in the following to define Ao = I and Hk = 0 for k < 0. Using (3) and (7) gives
k
j
. jHk-j-
(4)
where E{" } denotes expectation, C is the unit circle Izl = 1, and H/t = 0 for k < 0. It is clear from the definition (4) that R_/t = R~. The conditions of {u,} and the relation (3) give
H k=E{yt+kuTt},
Q,.ut-,..
+ Qout + Qxut- 1 + ... + Q,,ut-,,.
et =
-~
"'" A , , y t - .
On the other hand, if we made a particular guess for the Ai and Q~, we can predict the value 33t of the present output Yt from the values {Yt-1 . . . . . Yt-.; ut. . . . . ut-m} based the above relation as follows:
R k = E{yt+kyTt} = ~ H,+~H~
= ~j~j H(z)H(z) *?-ldz,
-- A 2 Y t - 2
+ Qout + QlUt-1 + ...
The z-transfer function matrix is defined by
(6)
If the approximant/4(z) is identical to the original system H(z), the inputs u, and outputs Yt are then related by a difference equation
12 = {{no, H i . . . . }l ~ I[nkl[ 2 < 00}. k=O
H(z) = ~ Hkz-/t
1+ ... + Qmz-,.}.
+
k Ut-k
k=m+l
j
)
AjHk-j ut-k.
(8)
Therefore, the following is obtained
E{l'et']2} = trlk~o(j~=oAjHk-~--Qk )
(5)
which implies H k is the cross correlation for timeshift k between the output {y,} and the input {u,}. In the approximation of a linear system {Hk}, the Pad6 approximation technique commonly uses finite matrices {Ho, H1 . . . . . H.,} of the impulse response sequence {Hk}~ as data. Throughout the paper, two kinds of data are given at the same time: the first ones are finite matrices {Ho. . . . . Hm} of the impulse response sequence {Hk} ~ and the second ones are finite matrices {Ro . . . . . R,} of the autocorrelation sequence {Rk}~-~. Assume an approximant H(z) of an original
×
Y
A~Hk_j
)wl •
(9)
Minimizing the above function (9) with respect to {Qo, Q1. . . . . e.,} yields
Qk= ~ AiHk-i j=O
for
O-
(10)
267
Approximation of multivariable linear systems
Definition. For the data {H o. . . . . H,~; Ro ..... Rn}, define a symmetric ( n + 1 ) x ( n + 1) [that is, (n + 1)r x (n + 1)r matrix] K(m, n) with the (i, j)th block element [K(m, n)]~j as follows:
autocorrelation sequence, K(m, n) is positive semidefinite. (2) Minimizing (9) with respect to {Qo. . . . . Qm} satisfies
m-j
[K(m, n)]ij = R j - i - ~ Hk+j-iH~.
(11)
Qk = £ AjHk-j for 0_< k_< m.
(14a)
j=o
k=O
where i, j = 0, 1, ..., n. When {Rk} is the autocorrelation sequence for the impulse response sequence {H~}, (11) is identical with
(3) If K(m, n) is positive semidefinite, then the minimization of (13) with respect to {AI . . . . , An} satisfies
[I, AI ..... An] K(m, n) = Fm.n[l, 0 . . . . . 03 (14b)
(12)
K(m, n)=
E
H'~-I
Lh,-,.J Lh,-,J which implies that K(m, n) is positive semidefinite. By means of (12), the substitution of (10) into (9) becomes
E{lle, ll2} =
tr [A o, A1 . . . . . An]
x K(m, n)[Ao, AI . . . . . An] T. (13) The minimization E{l]e, II2} with respect to {A1, A2 . . . . . An} requires the following lemma, which is verified by the Lagrange multiplier technique.
Lemma 1. Let K be a symmetric positive semidefinite n x n block (or nr x nr) matrix, and let be an n x 1 block (or nr x r) matrix. For the problem which is to minimize tr {XTKX} subject to OTX = I, a solution X exists and satisfies J(TK = ro T
where F is an r x r symmetric matrix uniquely determined by K, and satisfies
Furthermore, the solution ,Y satisfies
XTKX > XTK.~ for any X subject to OXX = I. Using the above arguments, the characterization of the least-squares approximation is obtained as follows.
Lemma 2. For the data {Ho . . . . . Hm; Ro . . . . . R~}, let K(ra, n) be the matrix defined by (11). (1) When the data {Ho . . . . , Hm; Ro, ..., Rn} are taken from an impulse response sequence and its
where Fm.n is a symmetric r x r matrix uniquely determined by the data. I shall call the problem of solving (14) the leastsquares approximation problem for the data {Ho, .... H~; Ro . . . . . Rn} with positive semidefiniteness of K(m, n). Some remarks will now be given here (14). Equation (14b) corresponds to the Yule-Walker equations in autoregression. All the properties derived in the sequel are based on (14). When matrix K(m, n) is nonsingular, the work of Mullis and Roberts in the scalar case could be extended to the multivariable case without essential difficulty. The case where K(m, n) is singular is rare in the scalar case but often occurs in the multivariable case. Equation (14b) may possess two solutions in this case. This fact should be taken into account in the singular case. Thus the questions are whether one can construct a stable approximant /~(z) which is uniquely determined from given data, and whether one can obtain a recursive algorithm for solving (14b) even when K (m, n) is singular. I shall show that these questions have been answered satisfactorily in this work. In the computation of (14b), however, the multiplications usually required are about proportional to n a, and the storage requirements are about proportional to n 2. Therefore, these become very large with increase of size n. Further reduction in computation time and storage is possible in solving (14b) because of the special form of K(m, n). A recursive algorithm in Section 3 will be obtained which is closely related to a multivariate version of Levinson's algorithm in autoregression. By substituting a solution {A 1..... An; Qo..... Qm} of(14) into (6), an approximant H(z) of the original system H(z) is obtained. The z-transfer function H(z) has an impulse response sequence {/tk} satisfying /~k= f 0
for k < 0
Qk -- £ AjI?-Ik-j j=l
for k _>-0.
(15)
268
YUJIRO INOUYE
where Qj = 0 for j > m. The approximant H(z) may be realized as a difference equation of order one in the following form
can be computed from the realization (16) as follows: Hk = ~gk for k < m - n ( CA k - l - ( m "~B f o r k > m - n + l
(16a)
-¢ct+1 = AYct + Bu,_(m_n)
ra--n
1~k = C A k I ( C T + ~ ISIi+kl?IT.
m--n
Ct = C2, + ~ HkU,-k
(18) (19)
i=0
(16b)
k=0
where ~, is an n-block vector (or nr-vector) defined by
-Pt+.-I- k=l~ H~-lU,÷.-k-"
Let / ~ ( m - 1, n - 1) be the matrix defined by (11) for the data {go,..., H,,- 1;/~o .... ,/~,- 1}. Equation (16c) and the same relation as (5) for g(z) then yield the fact that/~ is identical t o / ~ ( m - 1, n - 1). The following lemma, which is concerned with a Lyapunov equation, will be required later.
m
(16c)
Xt ~---
k=i
k=n
and the triplet {A, B, C} is given by
Lemma 3 (Appendix A). Let A and B be n x n and n x r matrices, respectively and let H ( z ) = (zI -- A ) - lB. (1) H(z) is stable if and only if there exist a symmetric positive semidefinite matrix K satisfying the Lyapunov equation K = A K A T + B B T.
(20)
-A
A=
I
(2) If all the eigenvalues of A lie in the open unit disk, then there exists at most one solution of (20). (3) Suppose H(z) is stable• Let
/
0
x
0
0
.•
. . . .
I
I( = ~ AkB(AkB) r.
(21)
k=p
.Bin - 1
B = L ]'Ira- n + 1
C =
[0 ... 0 1].
(16d)
In the case m > n, the above realization (16) of H(z) does not contain in the state space those states of pure delays which are serially connected to the input. Anoter realization can be presented which possesses all these states in its state space, but the realization (16) simplifies the arguments in the following. The stability of the approximant H(z) will be shown in Theorem 1. If the approximant g(z) is stable, then the state vector has a covariance matrix /~ satisfying /~ = E{:~t:~ T} = ~ AkB(AkB) T
(17a)
k=O
and the following Lyapunov equation I( = A I ( A T + B B T.
(17b)
Using the definition of the input sequence {Uk}, the impulse response and autocorrelation sequences
The matrix /( is then a solution of (20) which is symmetric and positive semidefinite. Furthermore, every other solution K of(20), that is symmetric and positive semidefinite, is ordered by K =>/(, namely, /~ is the minimal solution of (20) that is symmetric and positive definite. Some remarks will be given here on Lemma 3. The connections between stability and Lyapunov equations are discussed by Barnett (1971). Part 1 of Lemma 3 is a slight extension of Theorem 4.2 given by Barnett• Part 2 of Lemma 3 is similar to Barnett's Theorem 4.3. In the remainder of this section, the properties of the least-squares approximation problem will be given. The first important one is concerned with stability of the approximant H(z). Theorem 1. Assume the matrix K(m, n) defined by (11) is positive semidefinite. Then the approximant H(z) obtained by (6) and (14) is stable. Proof. Equation (14b) is equivalent to the following pair of equations K(m-
1, n - 1) = A K ( m -
1, n - 1)A x + BB r + BFrn.nBT
(22)
Approximation of multivariable linear systems
k=O
Ak[K(m,
n)]k,. = 0
H~=Hk
(23)
for
269 k = 0 , 1 . . . . . m.
(25)
where A and B are the matrices given by (16d), and /~ is the n x 1 block matrix given by
(2) The double sequence { F u } , whose element is uniquely determined by (14b) for (m, n), satisfies the following:
/~ = [I, 0 . . . O] T.
F..o = Ro - (HoH~ + . . . + H . H I )
(24)
This statement can be derived as follows. Equation (23) follows directly from equating the last block elements in both sides of (14b). From (11), one obtains [ K ( m - 1, n - 1)]i,j = [K(m, n)]i+lj+,
for i,j = 0 . . . . . n--1. This and the block companion form of A in (16d), together with (14b), yield [ A K ( m - 1, n - 1)]i a = [K(m, n)]ij+ 1
(26a)
F~ + ,,. < F.,.
(26b)
F.,. +, < F.,.
(26c)
rank Fro,. = rank K (m, n) - rank K (m - 1, n - 1), form>0, n>0.
(26d)
(3) Assume the data {Ho ..... Hm; Ro ..... R.} are taken from the impulse response and autocorrelation sequences of an original system H(z). Then the rank Fro,. is equal to r - p if and only if p is the largest integer such that for a p x r matrix Vofrank p and some matrices A1 ..... A., Qo..... Q., it holds that
for i,j = 0 . . . . . n - 1 V{y, + A l y , - x + ... + A,y,_,
[ A K ( m - 1, n - 1)AT]ij =
-
[K(m, n)]o,o - Fro,. for (i,j) = (0, O) [ K (m, n) ]ij otherwise. Also there holds from (11) [K(m, n)]ij = [K(m-
1, n -
- H . _ i H . _T j
1)]ij for i,j = 0 . . . . . n - 1 .
Then we obtain (22) from the last two equations. Conversely, it can be shown by an elementary calculation that equations (22) and (23) imply (14b). By (22) and the part 1 of Lemma 3, the z-transfer function ( z I - A)-*[B,/~x/Fm,.] is stable, which implies the stability of/-) (z) from the realization (16). This completes the proof of Theorem 1. A few remarks are given here. Note that the realization (6) may have unstable mode when it is not controllable (see the proof of Lemma 3 in Appendix A). This is not an issue in the scalar case, since the realization (6) is controllable if Fro,, # 0 or if Fro,. = 0 and Fro_ 1,.-1 # 0 (see Mullis and Roberts, 1976). Notice that, when K (m, n) is positive definite, the stability of H(z) could be also derived from the well-known minimum phase properties of predictors as in the work of Morf and co-workers (1977b). Theorem 2 (Appendix A). Presume the matrix K(m, n) is positive semidefinite. (1) The first m + 1 matrices of the impulse response sequence of the approximant H(z) are identical with those of the given data, that is
Qou, - . . . Q,,ut-,,} = o,
a.e.
for all t.
(27)
As a special case, the approximant/~(z) is identical with the original system H (z) if and only if F,,., = 0. Note that there may exist two different solutions of (14b) when the matrix K(m - 1, n - 1) is singular. The following theorem shows that these different solutions of (14b) produce the identical z-transfer function B(z). Theorem 3 (Appendix A). Presume the matrix K(m, n) is positive semidefinite. Then the z-transfer function ~(z) obtained by (6) and (14) is the same for all solutions of (14), that is, /~(z) is uniquely determined by the data {Ho . . . . . Hm; Ro . . . . . R.}. The approximant H(z) will now be examined to see how well matches the given data {Ho . . . . . Hm; Ro ..... R.}. From the part 1 of Theorem 2, the first m + 1 data are matched exactly. Some preparations will be made for the evaluation of differences between the data {Ro . . . . . R.} and the first n + 1 matrices of the autocorrelation sequence of the approximant H(z). Let
17(z) = A(z)- *z-m~/F.,.
(28a)
A(z) = I + AlZ- 1 + . . . + A.z-".
(28b)
where
We see from (22) and the part 1 of Lemma 3 that the system/-i7 (z) is stable. It will be said to be residual corresponding to the approximant/t(z). Let {/-iTk}
270
YUJIRO INOUYE
and {/~k} be respectively the impulse response and autocorrelation sequences of H(z). Note that ~k = 0 for k < m and ~,~ = ~/F,,,. Let ~(n - 1) be a symmetric n x n block matrix with the (i, j)th block element being/~j-i- Then iff(n - 1) is a block Toeplitz matrix whose first row block is [/~o, iffx. . . . . /~n-l]. It can be seen from the definitions that/~(n - 1) is identical to/~(m - 1, n - 1) defined by (11) for the data {0, ..., 0; /~o. . . . . /~n-1}. By means of the realization of H(z) [defined by (16) for the system H(z)], the part 3 of Lemma 3 and (17a), one obtains R ( n - 1) = / ~ ( m -
1, n - 1) = ~ AkBFm.nBT(AT) k. k=O
(29) The following lemma is also required. Lemma 4 (Appendix B). Let {x,} be a purely nondeterministic (wide-sense) stationary n-vector random process, and let K o = E { x , x [ } and K1 = E{xt+lx[} If the covariance Ko satisfies a Lyapunov equation
K o = AKo AT + BB r
m = n of the problem of computing a solution to (14b) may be handled by deriving an equivalent problem for which m' = n' = n. This technique was given by Mullis and Roberts (1976) in the scalar case and is given as follows. Denote the sequence {Ho, H~, H 2. . . . } by 9f ~, and define the operators O"L and aR by O'L~ = {HI, H2, H3,...} aR~,~ =
Let {Rk} be the autocorrelation sequence for 9f~, which is given by (4). For any positive integer k, the autocorrelation sequence for (as)kgf is the same as that of ~¢ itself. However, the autocorrelation sequence for (oL)k~ is given by k-1 j=0
Lemma 5 (Mullis and Roberts, 1976). For any 9f', let K(m, n) be defined by (11). Let
f a[m-n)gf~, if m>n
(30)
9f" =
and a condition
~,
if m = n
(33)
a~-m)gf ~, if m < n .
AKo=K 1
(31)
then Ko is the minimal solution of (30) that is symmetric and positive semidefinite. By means of these preparations, the evaluation of differences mentioned above is shown as follows. Theorem 4 (Appendix B). Presume the matrix K(m, n) is positive semidefinite. Let ~(z) be the approximant obtained by (6) and (14),/7(z) be the residual system defined by (28), and {/~k} and {/~k} be their respective autocorrelation sequences. (1) When the data {Ho . . . . . Hm; Ro . . . . . Rn} are taken from the impulse response and autocorrelation sequences of an original system H(z)
Rk--l~k=Rk
{0, Ho, H1 . . . . }.
for
k = 0 , 1 . . . . . n.
Define K'(n, n) for the sequence ~ ' by (11). Then K'(n, n) = K(m, n). The proof is made by an elementary calculation. We may compute a solution to the (m, n) problem by first deriving the data {H~ . . . . . H'n; R~ . . . . . R'n} and then applying the algorithm in the case m = n to get {AI, A2 . . . . . An}. {Qo, Q~ . . . . . Q=} is computed from (14a) by using the original data {no, ..., H,}. Now a recursive algorithm for computing a solution to (14b) for the case m = n is described. For notational simplicity, let K ( n ) = K ( n , n) and rn = r....
Algorithm. The data {Ho . . . . . Hu; Ro . . . . . Ru} are given. At the completion of step n, the quantities
(32)
A~(n), Bk(n), Ck(n): (0 <~k <%n), F,, f~, (2) Let A be the companion matrix (16d) for a solution {A 1. . . . . An} of (14b). If the eigenvalues of A all lie in the open unit disk (this assumption is attained if Fm,n > 0), the above relations (32) also hold.
have been calculated. To obtain step n + 1, the following calculations are carried out
A. = k=O ~ Ak(n)Rn+l-k
(34a)
An = - ~ Ak(n)Hn+ 1-k
(34b)
3. R E C U R S I V E A L G O R I T H M
This section will provide a re,curslve algorithm for computing a solution {A1 . . . . . An} to the equation (14b) for the case m = n. In the same way as the single input-output (or scalar) case, the case
k=O
O. = A.Bo(n)+A.Co(n )
(34c)
Approximation of multivariable linear systems (I)n = A.Co(n)T + Anf~n rn+,
(34d) (34e)
= rn - anO T - A.e T
fln + I -_ f~. + ~ . TF+. + I ~
(34t)
Ak(n + 1) = A~(n) -- A.B.+ t -k(n) -- A.C.+ 1-k(n) ] B~(n+ l ) - B k ( n ) - - ® . TF . ++ l A. + l - k ( n + l ) -
~
-
Ck(n + 1) = Ck(n) - ~ . TF . ++ 1 A . + , -k(n + 1) 0~
(34g)
where A.+,(n) = B.+,(n) = 0 and C,+t(n) = 0. As initial values for n --- 0, set A0(0) = I, r o = R o - H o HT, Bo(0) = r ~ Co(0 )
= H oT F o+, ~o = I +HToF~Ho .
As far as the matrices {Fo, F, . . . . . FM} are positive semidefinite (see Theorem 5 below), this procedure is carried out until step M is obtained. When K(M) is positive semidefine, this provides a solution {AI(M) . . . . . AM(M)} of (14b) for m = n --- M, and the approximant H(z) is given by (M
a(z)= tk=o 2
)-IM/M )
\
2 ! 2 A,(U)U-A/ k=Ok/=O (35)
If F. = 0 for some step n, the approximant ~(z) for step n is identical with the original system. If F~ is not positive semidefinite for step n, then K(n) is not positive semidefinite (see Theorem 5 below), which implies there is no system which matches the data {Ho . . . . . H.; Ro . . . . . R.}, that is, the data are never identical with those taken from the impulse response and autocorrelation sequences for any system with its impulse response sequence in space !s (see Theorem 6). The proof of validity of the algorithm is found in Appendix C. If one sets H~ = 0 for 0 ~< i ~
271
pseudo-inverse (Golub and Kahan, 1965; Dongarra and co-workers, 1979). The matrices A. and A. are the errors in the predicted values of R.+I and Hn+,, respectively, based on the nth approximation. If these both vanish, then the n + 1st approximation is equivalent to the nth. The Levinson algorithm in autoregression was also extended to non-Toeplitz cases (however closeto-Toeplitz cases) by Morf and co-workers (1977a) and Friedlander and co-workers (1979), where they employed the displacement rank properties. As shown below, the above algorithm can be used to examine whether the matrix K(m, n) for some data {Ho . . . . . Hm; Ro . . . . . R.} is positive semidefinite or not. Theorem 5 (Appendix C). Let K(M) be the matrix defined by (11) for some data {H o. . . . . HM; Ro. . . . . RM}, and {Ao(j) . . . . . A~(j); F j} be the quantities generated from the above algorithm (34) at step j. Then, K(M) is positive semidefinite if and only if Fo . . . . . F u are all positive semidefinite and {Ao(j) . . . . . Aft.j); F~} satisfies (14b) for m = n = j , where 0 ~
Yt = ~ HkUt-k + v,
(36)
k=0
where {Hk} •12,{ut} is the input sequence possessing the properties explained in Section 2, and {v,} is a stationary noise random process uncorrelated to {u,} with mean-zero. When some data {Ho . . . . . Hm; Ro . . . . . R.} are given, the problem considered here is to construct a system with noise (36) such that
E{yt+ku T } = H k for 0~
for
O<~k<~n.
(37a) (37b)
This is called the interpolation problem for the data {Ho, ..., Hm; R o, ..., R,}. The data is said to be consistent if and only if there exists a solution to *The least-squares approximation problem can also be formulatedfor the system(36) instead of (3), and obtain that it is characterizedby (14), where the data are computed by (37).
272
YUJIRO INOUYE
the problem. When there exists a solution to the problem, it is required that a system is constructed with noise such that the covariance E{v:~} of noise is as small as possible. When the covariance E{vt vXt} may be reduced to zero, the system obtained is said to match the data exactly. If a system H(z) can be represented in the form (6), it will be said to be of order at most (m, n). A solvability condition of the interpolation problem and a procedure for obtaining a system with noise that matches the data will be given. Theorem 6 (Appendix D). The data {H o. . . . . H~,; R0 . . . . . R.} for the interpolation problem is consistent if and only if the matrix K(m, n) defined by (11) for the data is positive semidefinite. Given the data {H 0. . . . . Hm; Ro .... , R.} for which K(m, n) >=O, there are two cases where Fro,. = 0 or Fro,. # 0. In the former case, there almost always exists a system/4 (z) which matches exactly the data. This is shown below. In the latter case the problem of artificially adding H= + 1 (or R. + 1) in a way that K(m + 1, n) > 0 and F r o + l , n = 0 [or K(m,n + 1) =>0 and Fro,. +1 = O] will be considered. Theorem 7 (Appendix D). Assume the positive semidefiniteness of K(m, n)for the data {Ho. . . . . Hm; R o. . . . . R.}. (1) There exists precisely one system H (z) of order at most (m, n) which matches exactly the data if and only if Fro,. = 0 [or equivalently rank K(m, n) < nr] (38a)
K ( m - 1, n - 1) = ~ AkBBT(AT) k
(38b)
k=O
where A is the companion matrix given by (16d) for a solution {At . . . . . A,} of the (m, n) least-squares approximation problem, and B is defined by B T = [-H T' Hm_ T t . . . . . H sT_ . + t ] .
(38c)
(2) When F r o . . # 0 [or equivalently rank K(m, n) > nr], there are at least two values of Hm+ 1 for which K(m + 1, n) >- 0 and Fm + t,. = 0 in the case rank Ro < q, and there are at least two values of R.+I for which K(m, n + 1) =>0 and Fro,.+ 1 = 0. Two examples are presented to illustrate the above results. These become counterexamples to the statements b and c of Proposition 1 provided by Mullis and Roberts (1976). Example 1. Consider the data { 1, 0; 2, 1} in the scalar case. It follows from the definition (11) that K(1, 1)=
LI 1]>0 1 =
which implies rank K(1, 1) = 1, and so F 1,1 = 0. The solution of(14b) for these data is given by A ~ = - 1. We obtain B = 0 in (38c), which implies the value of the right hand side of (38b) vanishes. Therefore, the condition (38b) does not hold for the data { 1, 0; 2, 1}. Part 1 of Theorem 7 shows that there is no system that matches exactly the data. This indicates that the statement b of Proposition 1 presented by Mullis and Roberts (1976) is not correct. It states that if K(m, n) is positive semidefinite but det K(m, n) = O, then there is precisely one system which matches exactly the data. The approximant/4(z) = 1 can be obtained for the data {1, 0; 2, 1}. By a simple calculation, it is seen that the following system with noise Vo matches the data Yt = Ut + V0 where mean E{vo} = 0 and variance E{v 2} = 1. Example 2. Consider the data { 1, 0; 2} in the scalar case. K(1, 0) becomes 1, and so the two values of H 2 are equal to + 1 for which F2. o = 0. It can be seen that H(z) = 1 _+ z -2 for the data {1, 0, + 1; 2}. The approximant H(z) matches exactly the data {1, 0, + 1; 2}. On the other hand, the two values ofR 1 become _+1 for which F 1A = 0. It follows that the condition (38b) does not hold for the data {1, 0; 2, +1} and thus there is no system that matches exactly the data { 1, 0; 2, 1}, or { 1, 0; 2, - 1}. This indicates that the statement c of Proposition 1 provided by Mullis and Roberts (1976) is not true. The approximants for the data {1, 0; 2, + 1} are same and are given by H(z) = 1. By a simple calculation, it follows that the following system with noise (_+ 1)tVo matches the data y, = u, + (+ 1)% where mean E{vo} = 0 and variance E{v2o} = 1. 5. CONCLUSIONS In this paper, the least-squares approximation problem was reduced to that of solving a linear (matrix-valued) equation. This equation can be solved recursively by the algorithm, which was closely related to a multivariate version of Levinson's algorithm in autoregression. A procedure for obtaining the solutions of the interpolation problem is presented through the solution of the least-squares approximation problem. A problem left unsolved however is that of finding the system of minimal order that matches exactly the data. Acknowledgements--The author
would like to thank Professor Y. Sakawa for constant encouragement and Dr H. Kimura for useful discussions.
273
Approximation of multivariable linear systems REFERENCES Barnett, S. (1971). Matrices in Control Theory. Von Nonstrand Reinhold. Boullion, T. L. and P. L. Odell (1971). Generalized Inverse Matrices. John Wiley. Dongarra, J. J., C. B. Moler, J. R. Bunch and G. W. Stewart (1979). L I N P A K Users" Guide. SIAM Press. Friedlander, B., M. Moff, T. Kailath and L. Ljung (1979). New inversion formulas for matrices classified in terms of their distance from Toeplitz matrices. Linear Algebra & Applic., 27, 31. Golub, G. and W. Kahan (1965). Calculating the singular values and pseudo-inverse of a matrix. J. S I A M Numer. Anal., Ser. B, 2, 205. Hutton, M. F. and 13. Friedland (1975). Routh approximations for reducing order of linear, time-invariant systems. IEEE Trans Aut. Control, AC-20, 329. Inouye, Y. (1981). Approximation of multivariabte linear systems with impulse response and autocorrelation sequences. Technical report no. 81-05, Department of Control Engineering, Osaka University, Japan. Ledwich, G. and J. B. Moore (1976). Minimal stable partial realization. Automatica, 12, 497. Markel, J. D. and A. H. Gray, Jr (1973). On autocorrelation equations as applied to speech analysis. IEEE Trans Audio Electroacoust., AE-21, 69. Morf, M., B. Dickinson, T. Kailath and A. Vieira (1977a). Efficient solution of covariance equations for l!near prediction. IEEE Trans Acoust., Speech, Signal Processing, ASSP-25, 429. Morf, M., D. T. Lee, J. R. Nickolls and A. Vieira (1977b). A classification of algorithm for ARMA models and ladder realization. Proc. 1977 IEEE Conf. on Acoust., Speech, Signal Processing. Mullis, C. T. and R. A. Roberts (1976). The use of second-order information in the approximation of discrete-time linear systems. IEEE Trans Acoust., Speech, Signal Processing, ASSP-24, 226. Rissanen, J. (1971). Recursive identification of linear system. S I A M J. Control, 9, 420. Rosenbrock, H. H. (1970). State-space and Multivariable Theory. Nelson. Shamash, Y. (1974). Stable reduced-order models using Pad& type approximations. IEEE Trans Aut. Control, AC-12, 615. Thether, A. J. (1970). Construction of minimal linear state variable models from finite input-output data. IEEE Trans Aut. Control, AC-15, 427. Whittle, P. (1963). On the fitting of multivariate autoregressions, and the approximate factorization of a spectral density matrix. Biometrika, 50, 129.
(1 - 12[2)v*Kv = IIv*e II2 >~0.
(A5)
If 2 is an unstable eigenvalue of A, that is, 12l > 1, then (A5) implies v*B = O. Therefore v*[a - 21, B] = 0
(A6)
rank [A - 21, B] < n.
(A7)
which yields
This fact means that 2 is an input-decoupling zero, namely, the mode associated with 2 is uncontrollable (Rosenbroek, 1970, p. 64). Hence this unstable mode is absent in the output behavior of H (z). Proof of(2) of Lemma 3 The uniqueness of the solution of (20) comes from the fact that if all the eigenvalues of A lie in the open unit disk, then the matrix equation K = A K A r has the unique trivial solution K=0. Proof of(3) of Lemma 3 The former part obviously follows from (AI) and (A2) above. The latter part comes as follows. Let
A(n) =
~
AkB(A~B) r.
(A8)
kfn+l
Then from (A8), (AI) becomes
g = ~ AkB(A~B)T + A(n).
(A9)
k=O
Since the sequence {A(n)} is monotonically decreasing a n d / ( exists, it follows lim A(n) = O.
(AI0)
n~ov
From (20), one obtains K = A K A T + BB T A K A • = A2K(A2) T + AB(AB) T A"K(A~)T = A" +l K ( A ~+t)T + A.B(AnB)T.
(All)
Hence
APPENDIX A
n
Proof of(I) of Lemma 3 If H(z) is stable, then there exists the limit ~ A~B(AkB)r as shown in (17a), which is denoted by k=o I( = ~ AkB(AkB) x.
(A1)
K = An+IK(A~+I)T+ ~ AkB(A~B)T.
(AI2)
kffiO
Combining (A9) and (A12) gives K - I?, = A n+ IK(An+ 1)r _ A(n).
(AI3)
k=O
Since K is positive semidefinite, (AI3) implies
Then (A1) becomes
g - / ~ _>_-A(n).
Letting n--* ~ in the above inequality together with (A10) yields K>=K.
kffi
Conversely, assume that (20) holds. Let 2 be an eigenvalue of A and let v* be a left eigenvector associated with 4. Then v*A = 2v*,
v* :~ 0.
(A3)
From (20) and (A3), it follows v*Kv = o*AKArv + v*BBrv = 1212v*Kv + v*BBrv
which means
(A14)
AJ'B(AtB) T A r = AI~A T + BB r. (A2)
l( = BB T + A
(A4)
Proof of Theorem 2
The relation (25) follows from (14a) and (15). Equation (26a) comes directly from the definition (11) of K(m, n) and the fact that Fm.o = K(m, 0). The relations (26b) and (26c) follow from Lemma 1, the definition (11) of K(m, n), and the fact that the solution to the (m, n) problem satisfies the constraints for the (m + 1, n) and (m, n + I) problems. Equation (26d) comes as follows: let {A t, A 2. . . . . A,; Fro.n} be a solution of (14b). Let T
274
YUJIRO INOUYE
be the (n + 1) x (n + 1) block matrix, whose diagonal elements are unity, defined by
T=[; A,.;.A.] lblOCknblocks.
The above two equations give
/~;-,]*
I
F/t',-1]
(A15)
..... / , J
[D1,
D,]
!
Then from (14b) and the definition (11) of K ( m , n) TK(m, n)T r =
It
K(m-
0
]
Since T is nonsingular, the above equation implies (26d). From (8) and (10), the following can be obtained in an argument similar to the derivation of (13) E { e , eT} = [A o. . . . . A . ] K ( m ,
, [ D 1. . . . .
Dn] T = 0
for t > m +
(A16)
1) "
1, n -
!
/~;-.J
n)[A o. . . . . A.] r.
1
which implies
• Dj/~',_j = 0
for t _->m + 1.
(A25)
./=1
(A17)
Assume/t~, = / t ~ for k = 0, 1. . . . . t - 1, where t ->_m + 1. Then from (15) and (A25), it follows that
This and (14b) give E { e , e~} = F.,..
(A18)
On the other hand, (7) yields
%
H t - H
""t =
_
~, (A'jI$1't_jA']H;'_j) --
i=i
e, = Yt - Y', = Y, + A I Y , + . . . + A,,yf_,, -- Qou, - . . .
- Q,,,u,_,..
It is clear that rank F.,. < r - p if and only if there exists a p x r matrix of rank p such that VFm.V r = 0. Therefore, rank F.,. < r - p if and only if V E { q e'~}vr = 0, or equivalently Ve, = 0 a.e. for all t, which implies (27). This completes the proof of Theorem 2.
=-
~Dj/t;_j=0
APPENDIX B
{At. . . . . A;} and {A~'. . . . . A~'}be two
solutions of (14b).
L e t D~ = A'~ - A':, where i = 1, 2 . . . . . n. From the two equations
obtained by substituting {A'1. . . . . A~} and {A~'. . . . . A~'} into {A~ . . . . . A,} of (14b), respectively, one obtains [D l .....
(A26)
This means the sequences { ~ } and { / ~ } are identical, which completes the proof of Theorem 3.
Proof of Theorem 3
Let
fort=>m+l.
j=l
D , , ] K ( m - 1, n - 1) = 0.
for
All the properties mentioned in Lemma 4 are not changed by a coordinate transformation, we can assume without loss of generality that the pair {A, B} is written in the following reachable canonical form
(A19)
Let /t'(z) and /t"(z) be respectively the z-transfer functions [obtained from (6) and (14)] corresponding to {A'l . . . . . A',,} and {A'~ . . . . . A'.'}, let {H~} and {/-/~'} be respectively their impulse response sequences, and let ~ ' ( m - 1, n - 1) be the covariance matrix given by (17a) for the system/~'(z). Part 1 of Theorem 2 gives
H;=H~=H;'
P r o o f o f L e m m a 4.
k = 0 , 1 . . . . . m.
(A20)
A=[A1
;:],
B = [ B1]
(A27)
where the pair {A 1, BI } is reachable (or controllable from origin). According to the above partition, let x, =
F
x~ '
[Kp
r 0 = LK2
K022_], K, -- LK**'
r~2J"
(A28)
From (30), (31), (A27), and (A28)
From (22)
Ko22 = A 2 K 022A2r A2 K 2 2 = K122.
K ( m - 1, n - 1) = A ' K ( m - 1, n - 1)A 'r + B ' B 'T + ~ F . , . ~ r.
(A21) where A' and B' are the matrices given by (16d) for the system /t'(z). Part 3 of Lemma 3 gives K ( m - 1, n - 1) > ~ A ' k ( B ' B 'T + BF,.,./~X)(A'r) k
(A29) (A30)
From (A29) and (A30), it follows K 22 = A2(K22) T.
(A31)
Using (A30) and (A31) yields
k=O
E{ (x2+1 - A2x2)(x2+ , -
> ~ A ' k B ' ( A ' k B ' ) "r = K ' ( m -
1, n -
1).
k=O
(A22)
D,,]K'(m-
1, n - 1)[Dr . . . . . D.] r = 0.
(A23)
= Ko - K1 A2T -- A2(K22) T -t- A 2 K 022 A T2 = 0
(A32)
~
I~;-1] fq;--1]T •
',-, L m-.l
x2+ 1 = A2 x2,
a.e. - oo < t < oo.
(A33)
Since {xf} is purely non-deterministic, this means x 2 = 0 a.e. for all t. Thus
On the other hand from (16c) [or 12] one obtains
g ' l m - 1, ~ - 1) =
A2x2)T}
22
which implies
This and (A19) yield [DI .....
22
.
K 2 = o, K o~' (A24)
=
o, K 2 = o.
(A34)
This yields, together with (30) and (A28) K ~ I = A I K o11 AIT + BxB~"
(A35)
Approximation of multivariable linear systems Since the pair {At, Bt} is reachable, by parts 1 and 2 of Lemma 3, it follows
275
Equation (23) which was derived from (14b) gives [K(m, n)]o,. = - ~'. A~[K(m, n)]k,.
K ~ I = ~ A~B~(A~BOT
(A47a)
kffiI
(A36)
kffiO
Similarly one obtains which means, together with (A27), (A28), and (A34) JR(m, n)]o.. = Ko = ~ AkB(AkB) r.
~ Ak[g(m, n)]k.~
(A47b)
kffil
(A37)
kffiO
~. = - ~ AtR.-k. By part 3 of Lemma 3, this implies that Ko is the minimal solution of the Lyapunov equation (30) that is symmetric and positive semidefinite. This completes the proof of Lemma 4. Proof of Theorem 4 Let xt be an n-block vector #oven by
(A47c)
k=l
From the definition (11), the last block columns of K(m, n) are respectively #oven by m-n
[ K ( m , n)]e,. = R . - e - Y, Hk÷.-,H~ for i = 0 . . . . . n
(A48a)
k=O m-n
~ nk-tUt+n-k
YI+n-I-
[/~(m, n)]~,. = R . _ / -
kffil
~ /-~k+._d-irk for i = 0 . . . . . n.
(A48b)
k=0
xt =
Yt+.-i - ~. Hk-eut+.-k
(A38)
Since H e = / l e for 0 -< i -< m, (A46)-(A48) yield
kffit
R. = / ~ . +/~..
Yt
Hk-nUt+s-k kffin
-- ~
Since the input sequence {ut} is stationary and purely nondeterministic, so the output sequence {Yt} and therefore the sequence {xt} defined by (A38) become. Let Ko = E{x, x,r} and KI = {xt÷ 1 xtr} • Then from an elementary calculation, (A38) and (11) give Ko = K(m
-
[K1]ea = [K(m, n)],a+~
1, n -
1)
(A39)
for i,j = 0, 1. . . . . n - 1. (A40)
For the #oven data, (22) and (23) are obtained which are equivalent to (14b), These, together with (A39) and (A40), ensure the relations Ko = AKo Ar + BB T + BFm,n/~T
(A41)
A K o = Kt.
(A42)
By Lemma 4, Ko is the minimal solution of (A41) that is symmetric and positive semidefinite. This yields by the part 3 of Lemma 3 K(m-
1,
n-
This implies the part 1. If the eigenvalues of the companion matrix A all lie in the open unit disk, the relation (A43) also follows from (22) and the part 2 of Lemma 3, which yields the relations (32) in the same argument as above. This completes the proof of Theorem 4. APPENDIX C Proof of validity of the algorithm (34) The following notation will be used. For a matrix K, Im K and Ker K denote, in the same order, the image of K and the kernel of K. K" denotes a reflexive generalized inverse of K, that is, a matrix that satisfies K K ' K = K and K ' K K " = K" (Bouillon and Odell, 1971). Using the definition (11) of K(n), it can be shown that with the proper definitions of the matrices E(n) and D(n)
/ffiO
kffiO
(A43) which implies from (17a) and (29)
I-"'7
V
K(n + l) = /
K(n)
I R-n- 1
1)= kffiO ~. Ak[B,g~(F.,.)][B, B~(F.,,.)]T(A~ = ~ AtBBT(AT~ + ~ AkBF.,.IIT(A~
(A49)
...
x [n.r+l ..... Hor]
=l=E(n) I D(n)r]. L D(n)
K(n)J
(A50)
Let A(n) = I-Ao(n). . . . . A.(n)]
(ASia)
B(n) = [Bn(n). . . . . Bo(n)]
(A51b)
C(n) = [C.(n) . . . . . Co(n)]
(A51c)
K(m - 1, n - 1) = / ( ( m - 1, n - 1) + / ( ( m - 1, n - 1). (A44) where A(n) satisfies (14b) for the case m = n, and F r o m t h e definition (11), the first block rows of K ( m - 1, n - 1) and K ( m - 1, n - 1) are respectively #oven by [K(m-
1, n - 1)]o.1 = Rj-
=-l-j kffiO m-l-j
~
(A52a)
C(n) = [ n l . . . . . nro]K(n) ".
(A52b)
Hk+1HT
for j = 0, 1. . . . . n - 1. (A45a) EK(m-l,n-l)]oa=l~1-
B(n) = [0 .... O, lJK(n)"
~
H~+~HT
Here K(n)" is a reflexive generalized inverse of K(n) recursively given by Lemma C2 below. The following two lemmas are required: /.emma C 1. Suppose K(n + 1) is positive semidefinite. Then
k=0
for j = 0, 1. . . . . n - 1. (A45b)
Im D(n) = Im K(n) Ker {[0 . . . . . 0, A.] + A . [ H .r . . . . . H~]} = Ker K(n)
(A53) (A54)
Since H~ = / ~ l for 0 ~ i ~ m, (A44), (A45), and (29) yield R~ = i~ + ~,j for j = 0, 1. . . . . n - 1.
(A46)
where A. and A. #oven by (34a) and (34b), respectively. Proof: It is easy to show that (A53) follows from (A50) and
276
YUJIRO INOUYE
the positive semidefiniteness of K(n + 1), From (A53) and the symmetry of K(n), it follows
which implies
I Bn(n)T ] Ker D(n) r ~ Ker K(n).
(A55)
:
=
[Adn ) ... An(n)]
LB,in:5
On the other hand, it holds [F., 0 . . . . . 0, A.] + An[H.X+ 1. . . . . HoT] = [A(n), 0]K(n + 1) = [A(n), 0][E(~2)
O(n)r] (A56) K(n) A"
It is easily shown that (A55) and (A56) imply (m4). Lemma C.2. Suppose A(n) and r . satisfy (14h) for the case m = n, and K(n) is positive semidefinite. Then a reflexive
generalized inverse of K(n) that is symmetric is given by K(n)" = A(n)TF~*A(n) +
D ( n - 1)T
.
(A66)
Equation (A59a) therefore becomes, from (A63), (A66), and (34a)(34c) O . = D(n)T B(n) T
A.Bo(n) + A.Co(n) - L A~(n)F.+ 1 -k
=
k=l
[01 0 ? 0
{ [i:l -,-- I -,l = ®. - L Ak(n)Fn+l-k"
K(n - 1)'
(A67)
k=l
By means of (A53) for n - 1, it is easy to show that K(n)" given by (A57) satisfies the definition of the reflexive generalized inverse. From (A50) and the above two lemmas, it is easy to derive the following:
A(n + 1) = [A.(n), 0] - A.[0, B(n)] - A.[0, C(n)]
(A58a)
B(n + 1) = I-0, B(n)] -OrF~++~A(n + 1)
(A58b)
C ( n + l ) = [0, C(n)] - @. "r F.+ + 1A(n + 1)
(A58c)
rn+, = r n - A.~).x - A.~.r
Here use is made of the fact Bo(n) = Bo(n)r. It follows from (A62b), (A50) and (A61) that [G.X ... G~x] = [C.(n) ... C,(n)]r(n - 1)+ Co(n)[R_.... R_ 1] -
O ( n - 1)r
where ~)., ~., ~). and ~ . are given as follows: (A59a)
~n = D(n)rC(n) T - H.+ ~
(A59b)
~,, = - A ( n + 1)[0, .... 0, I] 7
(A59c)
@n = - A ( n + 1)[H T"+~ . . . . . n ~ ] r.
(A59d)
(A68)
which means, together with (A64)
(A58d)
~ . = O(n)r B(n) r
fJ.[H .T. . . H r]
x
• Lc,in:J
= [Al(n ) .., A~(n)]
{?] R1
C ° ( n ) r - L_HI_J
:'1} ~1
. (A69)
Equation (A59b) therefore becomes, from (A63), (A69), (A61), and (34a), (34b), and (34d) n
Here the following equation is used which is derived from (A52) and (A54):
~) n = D(n)rC(n) T
H.+ 1 =
-
AnC0(n)
r
+ A.fl. - k~l Ak(n)Gn+ 1 -~
n
= On -- k~=tAk(n)Gn+ 1 -k.
[ 0 . . . 0 A.] + A.[nr~... H T] - {A.B(n) + A.C(n)}K(n) = { [ 0 . . . 0 A , ] + A n[H r . . . H ro] } {I --
K(n)'K(n) }
= 0.
(A60)
(A70)
Here use is made of the fact fl n = ll~. On the other hand, since
Define
A ( n + 1)K(n + 1) = [F.+I, 0 ... 0] n . = t + c(n)EtC .., HOr] ~.
(A71)
(A61) it follows from (A50), (A58), (A59a), and (A59b) that
Let [ F , + l , 0 ... 03 = [F,, 0 ... 0, A,] + A.[0, H~ ... H~]
[ F ~ . . . FOr] = B(n)K(n)
(A62a)
[Gr~ ... GOr] = C ( n ) r ( n ) - EH~ ... HT].
(A62b)
From the definition of D(n), it follows D(n)T = [ D ( n - 1)r, R , + , ] - H . + , [ H ]
- [A.O.r + A.t~T, A,B(n)K(n) + A.C(n)K(n)] which implies, together with (A62) IF,+ 1, 0 . . . 03 = IF,, 0 . . . 0, An]
... HTo]. (A63)
--T T - EA.e, + a . ¢- T, , AnEF.... For ] + A.[G~ ... GOr]].
(A73)
Therefore, it holds
Equations (14b) and (A50) give
D(n - 1)7 = - [Adn) ... A.(n)]K(n - 1).
(A64)
AnEF~ ...
F ,r ] + A,[G,T ... GI] = 0.
(A74)
Using (A67), (A70), and (A74), one obtains
It follows from (A62), (A50) and (A52) that
A . ( 6 r _ ®T) + An(~x _ @1-) = 0
[F~r
(A72)
FIr] = B(n) L k - . . . . R_,
(A75)
which yields, together with (A58d), Fn+ 1 = F. - A.@T - A.@T.
= [B.(n) ... n l ( n ) ] K ( n - 1)+ Bo(n)[R_,,... R_ 1]
- C o ( n ) r [ H ] . . . HI]
(A65)
(A76)
It follows from (A58a), (34b), (A52a), (A52b), and (A61) that
277
Approximation of multivariable linear systems T A(n + I)[H.+,, .... H ~ ] T = [A(n), 0 ] [ H . T+ I ..... HoT] T
where i, j, = 0, 1..... M. Then, we obtain from (14b)
- A.E 0, B(n)] [H.T+, ..... H ~ ] T -- A.[0, C(n)] [H.T+, ..... HoT] T
T(M)K(M)T(M) T = block diag {Fu, F u - , . . . . . Fo}. (A81)
= A(n)[HT.+, ..... H [ ] T - A.B(n)[H T ..... s T ] T
-- A.C(n)[HT. . . . . . HTo] T = -- A. - A.[0 ... 0 I]K(n)" [HT. ..... H T ] T -
A.C(n)[HT.. . . . .
ST] T
Since T(M) is nonsingular, the positive semidefiniteness of K(M) implies that of all Fo, F1 ..... Fu. Conversely, if {Ao(j).... , A~); F~} satisfies (14b) for m = n = j (where 0 < j ~ M), then (A81) holds. Therefore, if Fo, F1. . . . . F u are all positive semidefinite, then K(M) is also positive semidefinite.
= -- A.[0 ... 0 I]C(n) T - A.{I + C(n)[HT. ..... Ho~] ~} = - A.Co(n) T - A.f2.
which, together with (A59d) and (34d), implies ~. = @.. It follows from (A58a), (A59c) and (34c) that 0 n =
APPENDIX D
(A77)
- - A ( n -]- I)[0 ... 0
~/]
= A,[0, B ( n ) ] [ 0 ...
I] T
Proof of Theorem 6 If the interpolation problem for the data {Ho..... H,; Ro..... R,} is solvable, then there exists a system with noise (36). Since the random sequences {u,} and {v,} are uncorrelated, (36) and the definition (11) of K(m, n) give
+ A.[0, C(,)] [0... o z]~ = A.Bo(n) + A.Co(n) = O..
(A78)
Kim. ,) =
~
.n~_~
~-~
+ ~(,) (AS2)
k>ra+l
Using (A61), (A58c), (A59d), and the fact that ~. = (~., one obtains n.+, -f~. = {C(n + I ) - [0, C(n)]} In.T+, ..... HTo] T = - ~ .^r F . + t+A ( n + ~T
+
^
__
I)[HT+, . . . . . S ~ ] ~ T
+
= (b. F.+ ,~. - ~. F.+ ~ .
(A79)
which becomes (34f). This completes the derivation of the algorithm (34).
-
.Hk_..
u_ilk_..
where ~(n) is the positive semidefinite (n + 1) × (n + 1) block matrix with the (i, j)th block element being E{v,_, v~_~}. This implies that K(ra, n) is positive semidefinite. Conversely, assume K(m, n) is positive semidefinite. Let yf = ~ /t~u,_, + v, vt = ~, ~k~t-k + 6
Proof of Theorem 5 AS far as K(M) is positive semidefinite, {Ao(j). . . . . Aj(j); F j} satisfies (14b) for m = n =j, where O~_j<_M (see the above proof). Let T(M) be the upper triangular (M + 1)x (M + 1) block matrix defined by .... f0, for i > j [ltMLh'j = ~Aj_,(M - j ) ,
(A80) for i < j
(A83a)
k=O
(A83b)
k=0
where {/~,[ and {~,} are the impulse response sequences for H(z) and H(z), respectively, {6,} is a white noise process uncorrelated to {u,}, and {~,} is a deterministic process uncorrelated to {u,} and {at}. Then it can be shown that the system (A83) matches the given data. The details of the proof is found in the work of Inouye (1981). The proof of Theorem 7 is also found in the work of Inouye (1981).