A Recursive Algorithm for Linear System Identification - IEEE Xplore
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-34,NO.
492
3, JUNE 1986
A Recursive Algorithm for Linear System Identification
In Section I11 it will be shown that the least-squares minimization of the modified Kalman equation error based on the Wiener-Hopf technique [ 3 ] is equivalent to the Gram-Schmidt orthogonalization of intertwined Krylov [6] sequences and that this can be implemented, due to the skewness of K , by an order recursive procedure. This means that the algorithm also computes the lower order approximants compatible with the input-output record, I. INTRODUCTION Finally, some numerical examples are given in Section N recursive linear system identification [l], [ 2 ] ,one is IV . predominantly interested in obtaining a pole-zero model of an unknown system from measured records of the input 11. THE MODIFIED KALMANEQUATION ERROR and output. If u, ( t ) and u ( t ) are the time domain input to In what follows, wework in the real Hilbert space L,(O, and output of the system, respectively, then we are inter- T ) with scalar product ested in characterizing the impulse response h ( t ) by a sum T of complex exponentials, Le., (x,Y ) = x ( Q Y ( t ) dt, (3) n
Abstract-This paper deals with the pole-zero identification of a linear system from a measured input-output record. It is shown that the minimization of a modified version of the squared Kalman equation error can be implemented by an order recursive algorithm in the time domain. The algorithm is based on the Gram-Schmidt orthogonalization of intertwined Krylov sequences involving a skew self-adjoint linear operator.
I
5
0
h(t) =
C di exp (sit). i= 1
(1)
and the norm is Ix 1 = (x,x ) ~ ’ The ~ . adjoint of a linear K is denoted K t and defined as usual. The idenoperator In the sequel, n is referred to as the order of the system. I . tity operator is The quantities aiand diare the poles and the residues at If p ( z ) is the polynomial the poles, respectively. In the Laplace domain, the problem is to model the transfer function H(s)-which is the Laplace transform of h(t)-by a rational function as and K is a linear operator, then p ( K ) represents the operator Here V ( s ) and U,(s) are the Laplace transforms of the output and input, respectively. A(s) is a monic polynomial (leading coefficient = 1) of exact degree n and B ( s ) is a polynomial of degree n - 1. Equality in ( 2 ) is obtained when u ( t )and u, ( t )are noise free and the system order is exactly n. In this paper we will only consider real asymptotically stable systems, i.e., input, output, and impulse response are real functions and the poles are located in the left half-plane. The poles and residues are either real or occurring in complex conjugate pairs, and A(s), B ( s ) are polynomials with real coefficients. In Section I1 the identification problem is posed in terms of the least-squares minimization of a modified form of the Kalman equation error [ 3 ] - [ 5 ] .This will require the introduction of a skew self-adjoint linear operator K . Manuscript received May 23, 1985; revised October 23, 1985. The author is with the Laboratory of Electromagnetics and Acoustics, University of Ghent, Sint-Pietersnieuwstraat 41, B9000 Ghent, Belgium. IEEE Log Number 8407204.
n
Before proceeding, we need the following [ 5 ] proposition. Proposition I : If, in a linear vectorspace, n
u =
k= I
y(ak; ck),
(Tk,
ck
scalars,
(6)
and there is a linear operator K and a vector u such that Ky(a; c)
= oy(o; c)
+ cu
(7)
for all ak and ck under consideration, then there is a monic polynomial p (z) of exact degree n and a polynomial q (z) of degree n - 1 such that p(K)u
+ 4(K)u
=
0.
(8)
The proof is given in Appendix A. Proposition 1 is useful for linear system identification
0096-3518/86/0600-0492$01.OO O 1986 IEEE
purposes in the following way. Input u, (t)and outputu (t) of a linear system with impulse response h(t) are related by u(t)
=
Si
h(t -
7)
u,(r)
U(t) =
2
k= 1
dk exp ((Ykt)
1:
Putting
y ( a ; d ) = d exp ( a t )
s:
a IlYd7
it
(Y
+d
S0' U O d 7
+d
1'
y dr
T
(9)
u, d7.
(1 8) (19)
T
Adding both sides of (18) and (19) and dividing by 2 gives
y exp
- y(0) =
y - y(T) = d7
and, since h ( t ) is given by (l), n
y
(-0(k7)
U,(T)
dr.
(10)
-
$(Y(o) + y ( T ) ) = CVKY+ ~ K u ,
where the operator K is defined by
KX = 1 exp ( - a r ) u,(r)
d7,
(20)
2
(11)
[ 1' o
~ ( 7d7 )
+
d~].
(21)
We see thatK is a right-inverse of the operator d/dt (Le., d/dt K = I ) but it is a particular one, as is shown in the (10) can be written as next proposition. n Proposition 2: K is a skew self-adjoint operator (K = u = y ( a k ; dk) (12) -Kt). The proof is given in Appendix B. k= I Next we show that the term ( y (0)f y ( T ) ) in (20) where y (a;d ) satisfies vanishes for T + 03,provided u, (03) = lim, -, sU, (s) = 0. From (1 1) it follows that y(o) vanishes. This is equally d - y ( a ; d ) = ay(a;d ) duo. (13) true for y (03). Indeed, combining the requirement of dt asymptotic stability with the final value rule Equation (13), with K = d/dt, is of the canonical form y(03) = lim sY(s) (22) (7). This enables 'us to restate proposition 1 as follows. S-0 There are polynomials p ( z ) and q ( z ) such that gives [p(d/dt)]u + [q(d/dt)]u, = 0. (14) SUiJ0 ) y(03) = lim sY(s) = lim d This is not surprising, since (2) tells us that s-0 s-0 s -a
-
-4
+
4s) W ) -
B(s) UO(S)= 0,
(15)
-
a
which in the time domain translates to [A(d/dt)]u (t) - [B(d/dt)]U , ( t ) = 0.
=
A(s) V ( S )- B(s) U,(S)
(23)
S'O
(16) We may, therefore, if T = 03,rewrite (20) in'the canonical form (7) with
Hence, in this case, p ( z ) = A ( z ) , q (z) = - B (z). Formula (15) expresses the fact that, when no noise is present, the Kalman equation error [3j, [4]
E(s)
d lim sUO(s)= 0. --
(24)
c = -ad
(25)
Ku,.
(26)
u =
(17)
vanishes. When noise is present, most noniterative pole-zero modeling techniques utilize a minimum squared equation error criterion. This amounts to finding the coefficients of the polynomials A(s) and B ( s ) which minimize 1 E(s).(: where 1 I I is the equivalent norm in the Laplace domain. However, the fact that, in the time domain, one has to deal with thesuccessivederivatives of u ( t ) and u,(t), makes the operator d/dt numerically unattractive. A better way is to work with the successive integrals of input and output, as i s pointed out in [ 3 ] . This results in the minimization of a weighted version of the squared equation error., Here we follow a similar procedure. In order to obtain a skew self-adjoint integral operator K , we integrate (13 ) twice, oncewith lower limit zeroand once with lower limit T (reverse time integration). This gives, omitting the parameters a and d in y,
= l/a
0
Proposition 1 tells us that, in the absence of noise, the modified Kalman equation error e =P(K) u
+ q(K)
(27)
vanishes. If noise is present, we will minimize the squared modified equation error 1 e l2 with respect to the coefficients of the polynomials p ( z ) and q (z). 111. THE ALGORITHM Writing (27) explicitly gives n- I
e = K"u
n- 1
+ k = O pkKku + k = O qkKku.
(28)
The objective is to minimize I e 1' = ( e , e ) with respect to , n - 1. the coefficients P k , qk for k = 0 , 1, From elementary linearleast-squares approximation theory [ ' I ] ; we know that the minimum is obtained when
-
494
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. A S P - 3 4 , NO. 3, JUNE 1986
the following orthogonality conditions are satisfied:
The scalar a is such that (x;,
( e , Kku) = ( e , K%) = 0 for k = 0, 1,
,n
*
-
1.
(29)
The expressions (29) form the basis for the derivation of the algorithm. Definitions: A Krylov sequence [6] is a vector sequence of the form
-
( u , Kv, K’u, K3u,
*
(30)
j.
An intertwined Krylov sequence is a vector sequence of the form
( v , U , K v , Ku, K’u, K’u, K 3 u , K3u,
*
e).
x0 =
(xi,xj)
=
0
for i # j = 0 , 1,
( x n f 2 ,xj)
0
=
for j = 0, 1,
xjj
=
0
for j
(42)
=
0, 1,
,n - 3
(43)
and then determine the parameters by the remaining conditions for j = n - 2, n
xj) = 0
-
1, n, n
+ 1. (44)
Formulas (35) and (41) imply that + 2 , xj> = (Kxn> xj) = -
Kxj)
,n
for j = 0, 1,
u
-
3 (45)
since K is skew self-adjoint. Looking at (32) it is clear that applying the operator K on one of the xi augments its rank in the hierarchy by only 2, i.e., (32)
Kx~
xi+2
+ { x i + ] ,xi,
* *
. , X,)
(46)
and this is insufficient to produce a nonvanishing righthand side in (45). This proves (43). i f j, (33) Formulas (44) and (35) allow us to find expressions for the parameters ,8, y,and 8. We have , z } stands foralinear
and (xi, xj) = 0
+ 1.
. ,n
*
We first show that
(X,
*
, n + 1, (41)
*
then it is sufficient to prove that, for a particular choice of the parameters &, y n , and 8, in (35),
u
+ {u> x2 KU + (u, U > x3 = KU + (u,U , Ku} x4 = K’u + ( v , u , Kv, Ku) -XI =
(40)
0.
The proof is inductive. Suppose that
(31)
The optimal equation error e is therefore orthogonal to all vectors preceding the vector Knu in the intertwined Krylov sequence. If the vectors of the sequence (31) are linearly independent, then it is possible to orthogonalize them by a Gram-Schmidt orthogonalization (GSO) procedure [7], i.e., there are vectors xi such that
XI) =
for
-
where the expression {x, y , Combination of the vectors x, y, * ., 2 . (Kx,, xn-2) + 8,]x,-212 = 0 (47) Suppose the triangular GSO scheme (32) has been carried out. Because both xZnand e are of the form expressed (Kxn, X n - 1 ) + Y n l X n - l / 2 = 0. (48) by (28) and both are orthogonal to the same vectors K’u, K’u (i = 0 , 1, * * , n - l ) , andsince the GSO scheme The expression (Kx,, x,) vanishes identically because K is skew. Furthermore, is unique, it is clear that e
=
1
-
x2,.
The next proposition shows that the GSO scheme (32) can be implemented by a recursive procedure which takes advantage of the skewness of K . Proposition 3: The vectors xk satisfy the five-term recursive relationship xn+2 =
Kxn
+ Pnxn+1 n = 0,1,
f
7nXn-l
=
x,
eo = 81 = o
(36)
XI = u
+ au.
I2
= 0.
== - ( x n 7
Kxa - 2 )
(50)
+ /3,-2~,- I
= Kx,-~
(51)
gives (Kxn,xn-d
(37)
(49)
combined with
+ Yn-2Xn-3 + 8n-2Xn-4 (35)
= u
(Kxn, x , - 2)
X,
* * *
x-1 = x-2 = 0
PnlXn+l
Expressions (47) and (48) can be simplified, e.g., in (47)
+ enxn-2
with initial values Yo
(Kxn, x n + l ) +
(34)
= -(Xn,
Kxn-2)
-Ixn)2.
(52)
Similar considerations apply to (48). Defining the two intermediary quantities an and v,
(38)
un = Ixn12
(39)
vn =
(Kxn, xn+ 1 ) s
(53) (54)
the parameters
495
FOR LINEAR SYSTEM IDENTIFICATION
KNOCKAERT: ALGORITHM RECURSIVE
6, y,and 0 can be compactly written as
P,
-V,/W,+~
n = 0, 1,
y, =
V,,-~/W,,-~
n = 1 , 2,:
-
* *
(55)
-
(56)
en = W,/W,-2
*
n=2,3;**. Looking at (32), we see that x,, can be written as
x, = I I , ( K ) ~+ r, ( ~ ) u
TABLE Ia Preprocessing I.,
"I
(57) (58)
where l i l ,(z) and I ' ,( z ) are polynomials of degree [n/2], [(n - 2)/2], respectively ([n/2]stands for the integer part of d 2 ) . Inserting (58) in the recursion (35), and using the fact that this should be valid for all choices of u and u , leads to the following corollary. Corollary: The polynomials l il,( z ) and I?,(2) satisfy the five-term recursive relationship
X,+, = z x n + 6nXn-t with XP1 = X,, = 0 , and .
x,
=
1,
XI = a
(60)
for JJn ( 2 ) >
x, = 0 ,
"The index k runs through 0, 1,
for r n( 2 ) . Thispermitsa recursive construction of the relevant polynomials. The recursion (35), the initial values (36)i40), and the formulas (53)-(57).and (59)-(61) constitute the bulk of the algorithm. We now turn to the practical implementation of the algorithm. For data records terminated at time T, in.put and output are assumed to be identically zero from T to 03. This does introduce a truncation error which is small provided the data length T is such that u, (T) and u ( T ) are negligible. The integrations are hot performed by straightforwardRiemannsummation, but indirectly by first projecting input and output on the Legendre polynomials P , (a good reference concerning the use of Legendre polynomials and/or other orthogonal functions for identification purposes is [SI). This of course necessitates the linear transformation of the interval [0, TI into the interval [ - 1, 11 by the rule t = (7 1)T/2. The operator K , with f = Tl4, is then given by
. . . , m.
m
U(t)
k=O
1
1
dkpk(7)
dk = (k
0.5)
-1
U(t)
Pk(T)
d7,
(64) the algorithm can be run without numerical integration procedures, since (63) allows us to calculate Ku, or Ku by the rule rrn
where ck and
e,
ek
=
1
m
are related by
-fC1/1.5
em = fcm-j/(m - 0.5)
+
for k = 1, 2,
* * *
,m
- 1.
(66)
From Legendre polynomial theory [9], the scalar product (uo, is and has the following elegantproperty with respect to the Legendre polynomials (see Appendix C): The multiplicative factor 2fin (67) is redundant because formulas (55)-(57) show that only quotients of scalar products are of relevance. The algorithm used here for The advantage is that, after preprocessing input and out- finding the Legendre coefficients in (64) is based on trapezoidal rule integration and the properties of the associput as ated Legendre functions [9]. A Fortran listing is given in m 1 preprocessing u,(t) ckpk(7) ck = (k + 0.5) u o ( f ) pk(7) d7 Appendix D. Themainalgorithmandthe k=O -I phase are summarized in Table I.
c
/
496
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, SIGNAL AND
IV. NUMERICAL EXAMPLES Three examples are given to illustrate the technique. 1) Input: exp ( -24 Output: lO[exp (-t) - exp (- 2t)l Data: 201 uniformly spaced samples of a record of T = 8 3 1 Legendre coefficients calculated (m = 30). The true transfer function is
H(s) =
10
-
s
+ 1’
A first-order approximation yields 9.9982273
*(”
=s
+ 0.99992 14
’
2) Input: Output: Data:
exp ( - 3 t ) exp (- t) - exp ( -2t) 201 uniformly spaced samples of a record of T = 8 3 1 Legendre coefficients calculated ( m = 30). The true transfer function is
H(s)
=
PROCESSING, VOL. ASP-34, NO. 3, JUNE 1986
by orthogonalization than by passing to the normal equations and the Gram matrix, since the erroranalysis of the former reflects closer the perturbation theory of the leassquares problem itself [lZ]. A point much overlooked is that the Gram-Schmidt orthonormalization process is equivalent to the Cholesky decomposition of the inverse of the Gramian [ 5 ] ,without actually calculating the Gramian itself, and this means that, from a computational complexity point of view, orthogonalization procedures should be preferred. Moreover, the algorithm uses and is conditioned by a priori information about the skewness of K , information which would have partly been lost if one had simply constructed and processed the normal equations. Finally, preprocessing in terms of Legendre polynomials for system identification applications seems to be a promising approach, as indicated and advocated in 181. APPENDIX A If we define n
n
then p ( K ) can be factorized as
s + 3 s2 3s + 2 ’
n
+
.n( K - ail>
P(K) =
i = 1
A second-order approximation yields
+
0.9989215s 3.0652588 H(s) = s2 3.0008532s 2.0008076
+-
+
3) Input: exp (-2t) cos ( t ) Output:OS[exp ( - t ) - exp (-3t)l Data:201 equispaced data samples of a record of T = 7.1 3 1 Legendre coefficients calculated ( m = 30). The true transfer function is
H(s) =
s3
s2 + 4s + 5 + 6s2 + 11s + 6 ’
(A.2 )
since the identity operator commutes with any operator. Applying p ( K ) on both sides of n
u =
k= 1
Y(%;
Ck)
04.3)
yields n
Taking advantage of the commutativity of all factors in (A.2) and considering that ( K - o l ) y = cu is valid for all @k and ck under consideration, we may write n
A third-order approximatior, yields
+
1 . 0 0 0 6 0 5 3 ~+~ 3.9397009s 4.9362468 Finally, after defining the polynomial q ( z ) as H(s) = s3 + 5 . 9 9 9 2 6 9 9 ~+~ 10.9905812s 5.9921840’
+
This last example is borrowed from [ 101, There is close agreement between thetrueand approximate transfer functions, especially in the coefficients of the denominaaltor polynomial. In [lo] one sees a similar behavior, though the calculations there are based on the pencil-offunctions method.
n
it is clear that (A.5) is equivalent to
APPENDIX B V . CONCLUDINGREMARKS
Define two additional right-inverses K , and K2 of dldt A novel algorithm for continuous input-output pole- by zero modeling based on a least-squaresGram-Schmidt orthogonalization process has been presented.It is well Klx = X ( T ) dT (A4 known [ 111 that least-squares problems are better solved
s:
497
KNOCKAERT: RECURSIVE ALGORITHM FOR LINEAR SYSTEM IDENTIFICATION
c'
K2x =
-_
APPENDIX C X(T)
dT.
(A.9)
. T
+
4
Clearly, K = (K, K~). Integration per parts yields
sake the For of simplicity we take f = 0.5. From the theory of Legendre polynomials [9], we know that
d
PT
d
- P,+l(t) - - P,,-I(t) = (2n dt dt
(A. 10)
+ 1) P,(t).
(A.13)
Taking advantage of (A. 11)
1
T
(K,x, y)
= -
0
-
K2y xdt
+ [Klx -
(A. = 1 P,(1)
14)
and K2y]t.
(A.12)
Since [ K l x * K2y]: = 0, the conclusion is that K i - = -K2. Similarly, Ki. = -K1, and so K t = -K.
P,(-1) it is clear that P,+ I
- P,-1
= (-1)"
= (2n
+ 1) KP,.
APPENDIX D c This subroutine calculates the first m 1 Legendre coefficients 1 data-samples X ( j ) which are equispaced c D ( k ) from a record of n c in [ - 1, 11, endpoints included. It is based on trapezoidal intec gration and the properties of the associated Legendre functions. SUBROUTINE LEG(N, M,X,D) REAL*8 DR,.CS,S REAL*8 X ( 0 :500), QAl(0 :500), QA2(0 :500),CA(0 :500),D(0 : 100) QAl(N)=O.DO QA2(N) =O. DO DR =N/2. DO CS=l.DO DOKl=l,N-1
+
+
QAl(K1)=0.5DO*(-2.DO+Kl/DR)**2 QA2(Kl)=QAl(Kl) END DO
DO z=o, M
CA(N)=DR*QA2(N- 1) CA(0)= CS*CA(N)
cs=
-cs
IF(Z. EQ.0) QA2(0) =2. DO IF(Z. EQ.1) QA2(0)=2.M)13.DO IF(Z. NE.O.AND. I . NE. 1) QA2(0) =O. DO DOKl=l,N-l CA(Kl)=DR"(QA2(Kl+ 1) QA2(K1- 1) - 2. DO*QA2(K 1)) END DO S=O.DO DOKl=O,N s=s+ CA(K-l)*X(Kl) END DO D(Z)=S*(O.SW +Z) DOKl=l,N-1 CA(K1)=((2.D0*Z+1.D0)* 1(-l.DO+KlIDR)*QA2(Kl) - (1-2.DO)*QAl(K1))/(1+3.D0) END DO DOKl=l,N-1 QA1( K 1) = QA2(K 1) QA2(K 1) = CA(K 1) END DO END DO RETURN END
+
(A. 15) (A. 16)
498
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND
ACKNOWLEDGMENT Grateful acknowledgment ismadeto Dr. De Zutter for his encouragement and interest in this work. Appreciation goes to the anonymous reviewersfor their suggestions and
SIGNAL PROCESSING, VOL. ASP-34, NO. 3, JUNE 1986
[8] Paraskevopoulos, P. N. “Legendre identification approach series to and analysis of linear systems,” IEEE Trans. Automat. Contr., vol. AC-30, pp. 585-589,June1985.
L91 W. Magnus, F. Oberhettinger,and R. P. Soni, Formulas and Theorems for the Special Functions ojMathematica1Physics. New York: Springer-Verlag, 1966. COmnentS which allowed me to improve the Presentation [lo] T. K. Sarkar, 3. Nebat, D. D. Weiner, and v. K. Jain,“Suboptimal paper.of this approximation/identification electrofrom waveforms of transient magnetic systems by Pencil-of-Function method,” IEEE Trans. Antennas Propagat., vol. AP-28, pp. 928-933, Nov. 1980. REFERENCES [ l l ] G. H. Golub,“Numericalmethods for solvingleastsquaresproblems,” Numer.Math., vol. 7, pp.206-216,1965. C. R. Johnson, Jr,, “Adaptive IIR filtering:Currentresultsandopen Trans. Inform,Theory, vel. 1 ~ - 3 0pp, , 237-250, M ~ ~[I21 , G. W. Stewart, IntroductiontoMatrixComputations. New York: i s s u e s , ~IEEE ~ Academic,1973. 1984. L. Ljung and T. Soderstrom, Theory and Praciice of Recursive Identification. Cambridge, MA: MITPress,1983. T. K. Sarkar, S. A. Dianat, and D. D. Weiner,“A discussion of various approaches to the linear system identification problem,”IEEE Signal Processing, vol. ASSP-32, pp. 654Luc Knockaert (“81) received Physical the EnTrans. Acoust., Speech, gineer degree and Telecommunications Engineer 656, June 1984. degree from the State University of Ghent, BelR. E. Kalman, “Design of a self-optimizing control system,” Trans. gium, in 1974 and 1977, respectively. ASME, vol. 80, pp. 468-474, Feb. 1958. From 1979 to 1984 he was with the Universitk L. Knockaert, “Parametric modeling of electromagnetic waveNationale du Zaire. His current interests are in the forms,” Electromagnetics, vol. 4, no. 4, pp. 415-430, 1984. signal identification and electromagnetic problem A. S. Householder, The Theory of Matrices in Numerical Analysis. solving areas. New York: Dover, 1975. P. J . Davis, Interpolation and Approximation. New York: Dover, 1975.
492
3, JUNE 1986
A Recursive Algorithm for Linear System Identification
In Section I11 it will be shown that the least-squares minimization of the modified Kalman equation error based on the Wiener-Hopf technique [ 3 ] is equivalent to the Gram-Schmidt orthogonalization of intertwined Krylov [6] sequences and that this can be implemented, due to the skewness of K , by an order recursive procedure. This means that the algorithm also computes the lower order approximants compatible with the input-output record, I. INTRODUCTION Finally, some numerical examples are given in Section N recursive linear system identification [l], [ 2 ] ,one is IV . predominantly interested in obtaining a pole-zero model of an unknown system from measured records of the input 11. THE MODIFIED KALMANEQUATION ERROR and output. If u, ( t ) and u ( t ) are the time domain input to In what follows, wework in the real Hilbert space L,(O, and output of the system, respectively, then we are inter- T ) with scalar product ested in characterizing the impulse response h ( t ) by a sum T of complex exponentials, Le., (x,Y ) = x ( Q Y ( t ) dt, (3) n
Abstract-This paper deals with the pole-zero identification of a linear system from a measured input-output record. It is shown that the minimization of a modified version of the squared Kalman equation error can be implemented by an order recursive algorithm in the time domain. The algorithm is based on the Gram-Schmidt orthogonalization of intertwined Krylov sequences involving a skew self-adjoint linear operator.
I
5
0
h(t) =
C di exp (sit). i= 1
(1)
and the norm is Ix 1 = (x,x ) ~ ’ The ~ . adjoint of a linear K is denoted K t and defined as usual. The idenoperator In the sequel, n is referred to as the order of the system. I . tity operator is The quantities aiand diare the poles and the residues at If p ( z ) is the polynomial the poles, respectively. In the Laplace domain, the problem is to model the transfer function H(s)-which is the Laplace transform of h(t)-by a rational function as and K is a linear operator, then p ( K ) represents the operator Here V ( s ) and U,(s) are the Laplace transforms of the output and input, respectively. A(s) is a monic polynomial (leading coefficient = 1) of exact degree n and B ( s ) is a polynomial of degree n - 1. Equality in ( 2 ) is obtained when u ( t )and u, ( t )are noise free and the system order is exactly n. In this paper we will only consider real asymptotically stable systems, i.e., input, output, and impulse response are real functions and the poles are located in the left half-plane. The poles and residues are either real or occurring in complex conjugate pairs, and A(s), B ( s ) are polynomials with real coefficients. In Section I1 the identification problem is posed in terms of the least-squares minimization of a modified form of the Kalman equation error [ 3 ] - [ 5 ] .This will require the introduction of a skew self-adjoint linear operator K . Manuscript received May 23, 1985; revised October 23, 1985. The author is with the Laboratory of Electromagnetics and Acoustics, University of Ghent, Sint-Pietersnieuwstraat 41, B9000 Ghent, Belgium. IEEE Log Number 8407204.
n
Before proceeding, we need the following [ 5 ] proposition. Proposition I : If, in a linear vectorspace, n
u =
k= I
y(ak; ck),
(Tk,
ck
scalars,
(6)
and there is a linear operator K and a vector u such that Ky(a; c)
= oy(o; c)
+ cu
(7)
for all ak and ck under consideration, then there is a monic polynomial p (z) of exact degree n and a polynomial q (z) of degree n - 1 such that p(K)u
+ 4(K)u
=
0.
(8)
The proof is given in Appendix A. Proposition 1 is useful for linear system identification
0096-3518/86/0600-0492$01.OO O 1986 IEEE
purposes in the following way. Input u, (t)and outputu (t) of a linear system with impulse response h(t) are related by u(t)
=
Si
h(t -
7)
u,(r)
U(t) =
2
k= 1
dk exp ((Ykt)
1:
Putting
y ( a ; d ) = d exp ( a t )
s:
a IlYd7
it
(Y
+d
S0' U O d 7
+d
1'
y dr
T
(9)
u, d7.
(1 8) (19)
T
Adding both sides of (18) and (19) and dividing by 2 gives
y exp
- y(0) =
y - y(T) = d7
and, since h ( t ) is given by (l), n
y
(-0(k7)
U,(T)
dr.
(10)
-
$(Y(o) + y ( T ) ) = CVKY+ ~ K u ,
where the operator K is defined by
KX = 1 exp ( - a r ) u,(r)
d7,
(20)
2
(11)
[ 1' o
~ ( 7d7 )
+
d~].
(21)
We see thatK is a right-inverse of the operator d/dt (Le., d/dt K = I ) but it is a particular one, as is shown in the (10) can be written as next proposition. n Proposition 2: K is a skew self-adjoint operator (K = u = y ( a k ; dk) (12) -Kt). The proof is given in Appendix B. k= I Next we show that the term ( y (0)f y ( T ) ) in (20) where y (a;d ) satisfies vanishes for T + 03,provided u, (03) = lim, -, sU, (s) = 0. From (1 1) it follows that y(o) vanishes. This is equally d - y ( a ; d ) = ay(a;d ) duo. (13) true for y (03). Indeed, combining the requirement of dt asymptotic stability with the final value rule Equation (13), with K = d/dt, is of the canonical form y(03) = lim sY(s) (22) (7). This enables 'us to restate proposition 1 as follows. S-0 There are polynomials p ( z ) and q ( z ) such that gives [p(d/dt)]u + [q(d/dt)]u, = 0. (14) SUiJ0 ) y(03) = lim sY(s) = lim d This is not surprising, since (2) tells us that s-0 s-0 s -a
-
-4
+
4s) W ) -
B(s) UO(S)= 0,
(15)
-
a
which in the time domain translates to [A(d/dt)]u (t) - [B(d/dt)]U , ( t ) = 0.
=
A(s) V ( S )- B(s) U,(S)
(23)
S'O
(16) We may, therefore, if T = 03,rewrite (20) in'the canonical form (7) with
Hence, in this case, p ( z ) = A ( z ) , q (z) = - B (z). Formula (15) expresses the fact that, when no noise is present, the Kalman equation error [3j, [4]
E(s)
d lim sUO(s)= 0. --
(24)
c = -ad
(25)
Ku,.
(26)
u =
(17)
vanishes. When noise is present, most noniterative pole-zero modeling techniques utilize a minimum squared equation error criterion. This amounts to finding the coefficients of the polynomials A(s) and B ( s ) which minimize 1 E(s).(: where 1 I I is the equivalent norm in the Laplace domain. However, the fact that, in the time domain, one has to deal with thesuccessivederivatives of u ( t ) and u,(t), makes the operator d/dt numerically unattractive. A better way is to work with the successive integrals of input and output, as i s pointed out in [ 3 ] . This results in the minimization of a weighted version of the squared equation error., Here we follow a similar procedure. In order to obtain a skew self-adjoint integral operator K , we integrate (13 ) twice, oncewith lower limit zeroand once with lower limit T (reverse time integration). This gives, omitting the parameters a and d in y,
= l/a
0
Proposition 1 tells us that, in the absence of noise, the modified Kalman equation error e =P(K) u
+ q(K)
(27)
vanishes. If noise is present, we will minimize the squared modified equation error 1 e l2 with respect to the coefficients of the polynomials p ( z ) and q (z). 111. THE ALGORITHM Writing (27) explicitly gives n- I
e = K"u
n- 1
+ k = O pkKku + k = O qkKku.
(28)
The objective is to minimize I e 1' = ( e , e ) with respect to , n - 1. the coefficients P k , qk for k = 0 , 1, From elementary linearleast-squares approximation theory [ ' I ] ; we know that the minimum is obtained when
-
494
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. A S P - 3 4 , NO. 3, JUNE 1986
the following orthogonality conditions are satisfied:
The scalar a is such that (x;,
( e , Kku) = ( e , K%) = 0 for k = 0, 1,
,n
*
-
1.
(29)
The expressions (29) form the basis for the derivation of the algorithm. Definitions: A Krylov sequence [6] is a vector sequence of the form
-
( u , Kv, K’u, K3u,
*
(30)
j.
An intertwined Krylov sequence is a vector sequence of the form
( v , U , K v , Ku, K’u, K’u, K 3 u , K3u,
*
e).
x0 =
(xi,xj)
=
0
for i # j = 0 , 1,
( x n f 2 ,xj)
0
=
for j = 0, 1,
xjj
=
0
for j
(42)
=
0, 1,
,n - 3
(43)
and then determine the parameters by the remaining conditions for j = n - 2, n
xj) = 0
-
1, n, n
+ 1. (44)
Formulas (35) and (41) imply that + 2 , xj> = (Kxn> xj) = -
Kxj)
,n
for j = 0, 1,
u
-
3 (45)
since K is skew self-adjoint. Looking at (32) it is clear that applying the operator K on one of the xi augments its rank in the hierarchy by only 2, i.e., (32)
Kx~
xi+2
+ { x i + ] ,xi,
* *
. , X,)
(46)
and this is insufficient to produce a nonvanishing righthand side in (45). This proves (43). i f j, (33) Formulas (44) and (35) allow us to find expressions for the parameters ,8, y,and 8. We have , z } stands foralinear
and (xi, xj) = 0
+ 1.
. ,n
*
We first show that
(X,
*
, n + 1, (41)
*
then it is sufficient to prove that, for a particular choice of the parameters &, y n , and 8, in (35),
u
+ {u> x2 KU + (u, U > x3 = KU + (u,U , Ku} x4 = K’u + ( v , u , Kv, Ku) -XI =
(40)
0.
The proof is inductive. Suppose that
(31)
The optimal equation error e is therefore orthogonal to all vectors preceding the vector Knu in the intertwined Krylov sequence. If the vectors of the sequence (31) are linearly independent, then it is possible to orthogonalize them by a Gram-Schmidt orthogonalization (GSO) procedure [7], i.e., there are vectors xi such that
XI) =
for
-
where the expression {x, y , Combination of the vectors x, y, * ., 2 . (Kx,, xn-2) + 8,]x,-212 = 0 (47) Suppose the triangular GSO scheme (32) has been carried out. Because both xZnand e are of the form expressed (Kxn, X n - 1 ) + Y n l X n - l / 2 = 0. (48) by (28) and both are orthogonal to the same vectors K’u, K’u (i = 0 , 1, * * , n - l ) , andsince the GSO scheme The expression (Kx,, x,) vanishes identically because K is skew. Furthermore, is unique, it is clear that e
=
1
-
x2,.
The next proposition shows that the GSO scheme (32) can be implemented by a recursive procedure which takes advantage of the skewness of K . Proposition 3: The vectors xk satisfy the five-term recursive relationship xn+2 =
Kxn
+ Pnxn+1 n = 0,1,
f
7nXn-l
=
x,
eo = 81 = o
(36)
XI = u
+ au.
I2
= 0.
== - ( x n 7
Kxa - 2 )
(50)
+ /3,-2~,- I
= Kx,-~
(51)
gives (Kxn,xn-d
(37)
(49)
combined with
+ Yn-2Xn-3 + 8n-2Xn-4 (35)
= u
(Kxn, x , - 2)
X,
* * *
x-1 = x-2 = 0
PnlXn+l
Expressions (47) and (48) can be simplified, e.g., in (47)
+ enxn-2
with initial values Yo
(Kxn, x n + l ) +
(34)
= -(Xn,
Kxn-2)
-Ixn)2.
(52)
Similar considerations apply to (48). Defining the two intermediary quantities an and v,
(38)
un = Ixn12
(39)
vn =
(Kxn, xn+ 1 ) s
(53) (54)
the parameters
495
FOR LINEAR SYSTEM IDENTIFICATION
KNOCKAERT: ALGORITHM RECURSIVE
6, y,and 0 can be compactly written as
P,
-V,/W,+~
n = 0, 1,
y, =
V,,-~/W,,-~
n = 1 , 2,:
-
* *
(55)
-
(56)
en = W,/W,-2
*
n=2,3;**. Looking at (32), we see that x,, can be written as
x, = I I , ( K ) ~+ r, ( ~ ) u
TABLE Ia Preprocessing I.,
"I
(57) (58)
where l i l ,(z) and I ' ,( z ) are polynomials of degree [n/2], [(n - 2)/2], respectively ([n/2]stands for the integer part of d 2 ) . Inserting (58) in the recursion (35), and using the fact that this should be valid for all choices of u and u , leads to the following corollary. Corollary: The polynomials l il,( z ) and I?,(2) satisfy the five-term recursive relationship
X,+, = z x n + 6nXn-t with XP1 = X,, = 0 , and .
x,
=
1,
XI = a
(60)
for JJn ( 2 ) >
x, = 0 ,
"The index k runs through 0, 1,
for r n( 2 ) . Thispermitsa recursive construction of the relevant polynomials. The recursion (35), the initial values (36)i40), and the formulas (53)-(57).and (59)-(61) constitute the bulk of the algorithm. We now turn to the practical implementation of the algorithm. For data records terminated at time T, in.put and output are assumed to be identically zero from T to 03. This does introduce a truncation error which is small provided the data length T is such that u, (T) and u ( T ) are negligible. The integrations are hot performed by straightforwardRiemannsummation, but indirectly by first projecting input and output on the Legendre polynomials P , (a good reference concerning the use of Legendre polynomials and/or other orthogonal functions for identification purposes is [SI). This of course necessitates the linear transformation of the interval [0, TI into the interval [ - 1, 11 by the rule t = (7 1)T/2. The operator K , with f = Tl4, is then given by
. . . , m.
m
U(t)
k=O
1
1
dkpk(7)
dk = (k
0.5)
-1
U(t)
Pk(T)
d7,
(64) the algorithm can be run without numerical integration procedures, since (63) allows us to calculate Ku, or Ku by the rule rrn
where ck and
e,
ek
=
1
m
are related by
-fC1/1.5
em = fcm-j/(m - 0.5)
+
for k = 1, 2,
* * *
,m
- 1.
(66)
From Legendre polynomial theory [9], the scalar product (uo, is and has the following elegantproperty with respect to the Legendre polynomials (see Appendix C): The multiplicative factor 2fin (67) is redundant because formulas (55)-(57) show that only quotients of scalar products are of relevance. The algorithm used here for The advantage is that, after preprocessing input and out- finding the Legendre coefficients in (64) is based on trapezoidal rule integration and the properties of the associput as ated Legendre functions [9]. A Fortran listing is given in m 1 preprocessing u,(t) ckpk(7) ck = (k + 0.5) u o ( f ) pk(7) d7 Appendix D. Themainalgorithmandthe k=O -I phase are summarized in Table I.
c
/
496
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, SIGNAL AND
IV. NUMERICAL EXAMPLES Three examples are given to illustrate the technique. 1) Input: exp ( -24 Output: lO[exp (-t) - exp (- 2t)l Data: 201 uniformly spaced samples of a record of T = 8 3 1 Legendre coefficients calculated (m = 30). The true transfer function is
H(s) =
10
-
s
+ 1’
A first-order approximation yields 9.9982273
*(”
=s
+ 0.99992 14
’
2) Input: Output: Data:
exp ( - 3 t ) exp (- t) - exp ( -2t) 201 uniformly spaced samples of a record of T = 8 3 1 Legendre coefficients calculated ( m = 30). The true transfer function is
H(s)
=
PROCESSING, VOL. ASP-34, NO. 3, JUNE 1986
by orthogonalization than by passing to the normal equations and the Gram matrix, since the erroranalysis of the former reflects closer the perturbation theory of the leassquares problem itself [lZ]. A point much overlooked is that the Gram-Schmidt orthonormalization process is equivalent to the Cholesky decomposition of the inverse of the Gramian [ 5 ] ,without actually calculating the Gramian itself, and this means that, from a computational complexity point of view, orthogonalization procedures should be preferred. Moreover, the algorithm uses and is conditioned by a priori information about the skewness of K , information which would have partly been lost if one had simply constructed and processed the normal equations. Finally, preprocessing in terms of Legendre polynomials for system identification applications seems to be a promising approach, as indicated and advocated in 181. APPENDIX A If we define n
n
then p ( K ) can be factorized as
s + 3 s2 3s + 2 ’
n
+
.n( K - ail>
P(K) =
i = 1
A second-order approximation yields
+
0.9989215s 3.0652588 H(s) = s2 3.0008532s 2.0008076
+-
+
3) Input: exp (-2t) cos ( t ) Output:OS[exp ( - t ) - exp (-3t)l Data:201 equispaced data samples of a record of T = 7.1 3 1 Legendre coefficients calculated ( m = 30). The true transfer function is
H(s) =
s3
s2 + 4s + 5 + 6s2 + 11s + 6 ’
(A.2 )
since the identity operator commutes with any operator. Applying p ( K ) on both sides of n
u =
k= 1
Y(%;
Ck)
04.3)
yields n
Taking advantage of the commutativity of all factors in (A.2) and considering that ( K - o l ) y = cu is valid for all @k and ck under consideration, we may write n
A third-order approximatior, yields
+
1 . 0 0 0 6 0 5 3 ~+~ 3.9397009s 4.9362468 Finally, after defining the polynomial q ( z ) as H(s) = s3 + 5 . 9 9 9 2 6 9 9 ~+~ 10.9905812s 5.9921840’
+
This last example is borrowed from [ 101, There is close agreement between thetrueand approximate transfer functions, especially in the coefficients of the denominaaltor polynomial. In [lo] one sees a similar behavior, though the calculations there are based on the pencil-offunctions method.
n
it is clear that (A.5) is equivalent to
APPENDIX B V . CONCLUDINGREMARKS
Define two additional right-inverses K , and K2 of dldt A novel algorithm for continuous input-output pole- by zero modeling based on a least-squaresGram-Schmidt orthogonalization process has been presented.It is well Klx = X ( T ) dT (A4 known [ 111 that least-squares problems are better solved
s:
497
KNOCKAERT: RECURSIVE ALGORITHM FOR LINEAR SYSTEM IDENTIFICATION
c'
K2x =
-_
APPENDIX C X(T)
dT.
(A.9)
. T
+
4
Clearly, K = (K, K~). Integration per parts yields
sake the For of simplicity we take f = 0.5. From the theory of Legendre polynomials [9], we know that
d
PT
d
- P,+l(t) - - P,,-I(t) = (2n dt dt
(A. 10)
+ 1) P,(t).
(A.13)
Taking advantage of (A. 11)
1
T
(K,x, y)
= -
0
-
K2y xdt
+ [Klx -
(A. = 1 P,(1)
14)
and K2y]t.
(A.12)
Since [ K l x * K2y]: = 0, the conclusion is that K i - = -K2. Similarly, Ki. = -K1, and so K t = -K.
P,(-1) it is clear that P,+ I
- P,-1
= (-1)"
= (2n
+ 1) KP,.
APPENDIX D c This subroutine calculates the first m 1 Legendre coefficients 1 data-samples X ( j ) which are equispaced c D ( k ) from a record of n c in [ - 1, 11, endpoints included. It is based on trapezoidal intec gration and the properties of the associated Legendre functions. SUBROUTINE LEG(N, M,X,D) REAL*8 DR,.CS,S REAL*8 X ( 0 :500), QAl(0 :500), QA2(0 :500),CA(0 :500),D(0 : 100) QAl(N)=O.DO QA2(N) =O. DO DR =N/2. DO CS=l.DO DOKl=l,N-1
+
+
QAl(K1)=0.5DO*(-2.DO+Kl/DR)**2 QA2(Kl)=QAl(Kl) END DO
DO z=o, M
CA(N)=DR*QA2(N- 1) CA(0)= CS*CA(N)
cs=
-cs
IF(Z. EQ.0) QA2(0) =2. DO IF(Z. EQ.1) QA2(0)=2.M)13.DO IF(Z. NE.O.AND. I . NE. 1) QA2(0) =O. DO DOKl=l,N-l CA(Kl)=DR"(QA2(Kl+ 1) QA2(K1- 1) - 2. DO*QA2(K 1)) END DO S=O.DO DOKl=O,N s=s+ CA(K-l)*X(Kl) END DO D(Z)=S*(O.SW +Z) DOKl=l,N-1 CA(K1)=((2.D0*Z+1.D0)* 1(-l.DO+KlIDR)*QA2(Kl) - (1-2.DO)*QAl(K1))/(1+3.D0) END DO DOKl=l,N-1 QA1( K 1) = QA2(K 1) QA2(K 1) = CA(K 1) END DO END DO RETURN END
+
(A. 15) (A. 16)
498
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND
ACKNOWLEDGMENT Grateful acknowledgment ismadeto Dr. De Zutter for his encouragement and interest in this work. Appreciation goes to the anonymous reviewersfor their suggestions and
SIGNAL PROCESSING, VOL. ASP-34, NO. 3, JUNE 1986
[8] Paraskevopoulos, P. N. “Legendre identification approach series to and analysis of linear systems,” IEEE Trans. Automat. Contr., vol. AC-30, pp. 585-589,June1985.
L91 W. Magnus, F. Oberhettinger,and R. P. Soni, Formulas and Theorems for the Special Functions ojMathematica1Physics. New York: Springer-Verlag, 1966. COmnentS which allowed me to improve the Presentation [lo] T. K. Sarkar, 3. Nebat, D. D. Weiner, and v. K. Jain,“Suboptimal paper.of this approximation/identification electrofrom waveforms of transient magnetic systems by Pencil-of-Function method,” IEEE Trans. Antennas Propagat., vol. AP-28, pp. 928-933, Nov. 1980. REFERENCES [ l l ] G. H. Golub,“Numericalmethods for solvingleastsquaresproblems,” Numer.Math., vol. 7, pp.206-216,1965. C. R. Johnson, Jr,, “Adaptive IIR filtering:Currentresultsandopen Trans. Inform,Theory, vel. 1 ~ - 3 0pp, , 237-250, M ~ ~[I21 , G. W. Stewart, IntroductiontoMatrixComputations. New York: i s s u e s , ~IEEE ~ Academic,1973. 1984. L. Ljung and T. Soderstrom, Theory and Praciice of Recursive Identification. Cambridge, MA: MITPress,1983. T. K. Sarkar, S. A. Dianat, and D. D. Weiner,“A discussion of various approaches to the linear system identification problem,”IEEE Signal Processing, vol. ASSP-32, pp. 654Luc Knockaert (“81) received Physical the EnTrans. Acoust., Speech, gineer degree and Telecommunications Engineer 656, June 1984. degree from the State University of Ghent, BelR. E. Kalman, “Design of a self-optimizing control system,” Trans. gium, in 1974 and 1977, respectively. ASME, vol. 80, pp. 468-474, Feb. 1958. From 1979 to 1984 he was with the Universitk L. Knockaert, “Parametric modeling of electromagnetic waveNationale du Zaire. His current interests are in the forms,” Electromagnetics, vol. 4, no. 4, pp. 415-430, 1984. signal identification and electromagnetic problem A. S. Householder, The Theory of Matrices in Numerical Analysis. solving areas. New York: Dover, 1975. P. J . Davis, Interpolation and Approximation. New York: Dover, 1975.
