- Email: [email protected]
Statistical Identification of Parallel Cascades of Linear and Nonlinear Systems
Copyrigh t © IFAC Ide ntificatio n a nd System Pa rame te r Estimation 1982. Washington D.C" USA 1982
STATISTICAL IDENTIFICATION OF PARALLEL CASCADES OF LINEAR AND NON LINEAR SYSTEMS M. Korenberg R ogers, Bereskin and Parr, P. O. Box 313, Commerce Court Wes t, Suit e 3710, T oronto, On tario, Canada
Abstract. An orthogonal statistical method is proposed for identifying a nonlinear s\'stem by successively adding a basic model paralle l to si~i l ar such models. The basic model selected must not be closed under acldition. A first such model is obtained to approximate the output of the non l inear s\'stem . The difference between the system output and the first mode l output is treated as a new output, and a second model is obtained to approxirate the l atter output . From this output the second model ou t put is subtracted to produce a third output, and a third model is obtained to approximate t'lis output, and so on. The method is illustrated when the nonlinear system is assumed to have a \\'iener functional expansion, and the basic model is a cas cade of a linear system and a static nonlinearity. It is shown that the nonlinear system may be approximated to an arbitrary accuracy by the sum of the models obtained . and that a functiona l expansion or a differentia l c~uation expansion may then be wri t ten down for the nonlinear system . An explicit time domain so l ution is a l so proposed for identifying a cascade of a linear. a static non linear and a linear system, which gives the impulse rpsronses direct l y by specia l first order cross corre l ations. Keywords. Nonlinear systems; identification; models; correlation net~ods; discrete systems; parameter estimat ion; time-domain analysis; cascacles. I\TRODUCTIO\ A typical approach to identifying a nonlinear system is to approximate its outpu t by means of some basic model, for e xample a Vo l terra functional expansion. Once the best such expansion is obtained, in the least-square sense, the repre sentation cannot be improved by adding a second such expan s ion paral l e l to the first, because the sum of an\' two Vo l terra series is a Volterra series. The situation is different, however, when the basic model chosen is not closed under addition. for example, suppose t hat the basic mode l -:omDrises a linear s;'stem follOl,ed by a static nonlinearit\' in the form of a pol:~omial or ;-'Oher series, and that the "best" impulse response and coefficients def i ning the nonlinearit\· are found to approximate the nonlinear system . The residual, i.e. the dif ference betl,een the svstel'1 output and the model output, can itse l f be treated as a neh output, and a second mode l of a l inear system followed b\' a static nonlinearity can be ob tained to approximate the neh output. The sum of the tlW model O:ltputS I,ill be c l oser to the system output, in the mean square sense, than is the first mode l ou t put a l one. The difference betheen the system output and the sum of the thO model outputs may be treated as a third output, and a third such
669
model obtained to approximate this output. The procedure may be repeated indefinitely to approximate the nonlinear system to an arbitrary accuracy in the mean square sense by a sum of the mode l s obtained separately. This procedure is predicated on the fact that the sum of two such models is not equi\'alent to a sing l e mode l of the same type. In other hords, it is necessary (but not sufficient) that the basic model selected not be closed under addition. An advantage of the proce dure is that, while the final representation comprises a parallel array of the models, at each stage of t he identification on l y a sing l e mode l is computed. ~loreover the procedure is orthogonal in tha t an additional mode l may be identified to improve the appro ximation hithout recalcu l ating the previous models obtained. Although there are many possible choices for the basic mode l , the procedure wi l l be illustrated here I,hen the mode l is a cascade of a l inear system fol l ohed by a static nonl·inearity. BOI,'ever it I,'ill be appreciated that the me th od proposed is not limited to successive approximat ions by the cascade model, and that other models not closed under addi t ion may be substi t uted.
M. Korenberg
670
Before illustrating the proposed method, however, the identification of a cascade of a linear, a static nonlinear, and a linear system (Spekreijse, 1969; Spekreijseand Oosting, 1970 ; Korenberg, 1973a, 1973b) is re-examined. It is shown that the i~pulse response of each linear system can be obt~ined directly as a special first order cross correlation. Next the proposed method is illustrated when the basic model comprises a linear system followed by a static nonlinearity, the extension to a basic model having a second linear system following the nonlinearity being trivial. Finally, it is shown that a functional expansion or R differential equation expansion can be obtained immediately from t~e parallel casc ade r epresentation identified.
A similar formula for IK(s) I then follows from (3). Analogous formulas exist for the phase of G(s) and K(s). Thus the linear systems may Qe identified from the first and second order cross correlations (and the static nonlinearity then determined). Alternately, setting sl=-s2=s in (4) and neglecting the constant IG(s) I =
2[Y le~::_S) 'j ~ xx~(O)__
However the latter formula for IG(s) I is less general than (6) and cannot be used if K(O)=O, which could occur, for example, if the second linear system in Fig. 1 involves a differentiator. y i e !~s
Setting ° =° =0 in (2) 1 2 4>XXY(O,O) = A2
CASCADE OF A LINEAR, A STATIC NONLINEAR, AND A LINEAR SYSTEM Spekreijse (1969) suggested an interesting nethod for identifying a cascade of a linear, a static nonlinear, and a linear system (Fig. 1) using special test inputs and re~eated experiments. Korenberg (1973a, 1973b) nroposed a correlation technique for explicitly identifying each corronent in the cascade following a single "aussian white-noise experiment. The latter technique involves evaluating the time-averaged first order cross correlation y(t)x(t-o) Alfoo k(T)g(o- T)dT
o
(1)
and the time-averaged ' second order cross correlation, shown to equal 4>XX/ O l, 0 2 ) = A2 foo k( T) g(01-T)g(02-T)dT ,
o
(2)
"'here x(t) is a \yhite "aussian input, get) and k(t) are the impulse responses of the first and second linear systems r espectivel y , Al and A2 are constants, yet) is the cascade output, and the mean of the output is subtracted before the second order cross correlat ion is comput ed (c.f . l1armarelis and Harmarelis, 1978). Taking the Laplace trans:Oorm of (1) and the two-dimensional Laplace transform of (2) yields respectively ~
and Hence
xy (s) = AlK(s )r;(s)
(3)
xXy(S l ,s2) = A2 K(sl+s2)G(sl)G(s2) (4) G(sl)G(s2) G(s 1+s2)
Al xxv(s l,s2 ) =:\2
(5)
and. letting sl=-s/2, s2=s, and neglecting arbltrarv constants . 1<1> (-s/2,s)1 IG(s) I = xxv ( 6)
1~y(s/2)I
(7)
J:
l:(T)i(o-T) dT
(8)
Billings and Fakhouri (lQ7U) have shown elegantly that by expressinp, (1) and (8) in sampled data form and assuming the sampled correlation functions am', the linear systems can be fit by rational pu1,sP' transfer functions . the linear systems CRn be identified. The cascade of Fig . 1 has indeed been extensively studied, e.g. Lammers and De Boer (1979), Bi 11 ings (1980). When no assumptions are made about the form of the linear systems, usin p 4> x (0,0) in place of 4> (0 ,0 2 ) may lose f06 much in1 formation f3 Yunlquely identify the linear systems. For example, if pet) and k(t) each assume only the values 0 and I but are different functions, then using (1) and (8) it will not be possible to det ermi ne the order in which the linear systems appear in the cascade. This is because g2(t)=g(t) and k 2 (t)=k(t), so that g anc, I: may be interchanged in (1) and in (8). This problem never arises using (1) and (2), or e quivalentl y (3) and ( 4), so that get) and k(t) can then be uniquely identified up to constants. It is now shown that get) <'.nd k(t) may be directl y determined in the time domain from 4> (0) and a special cross correlation. Let
4>~r(t) be the inverse of 4> (t), so that the c3Xvolution of these two fGrictions is the unit impulse of Dirac. Form the linear system L whose impulse response is the conv01u-at -1 tion of e and 4> (t) (for some a>O), i.e. xy e- at * 4> -l(t) xy and feed the output y of the cascade of Fig. 1 through L. Let pet) be the corresponding output of L. It will be proved that g(o) is given by the time average p(t)x(t)'x(t- o )
,°
(9)
Let 4> xp (0) and 4> xxp (° 1 2 ) be respectively the first and second order cross correlations of x and p defined in the usual manner. Note that the linear system L will have transfer
Parallel Cascades of Linear and Nonlinear Systems
671
GENERAL NONLINEAR IDENTIFICATION
function (s+a)1>
(s) xy Therefore, analogous to (3), neglecting arbitrary constants.
For convenience, the identification method will be demonstrated for a discrete nonlinear system N having finite memory of length M and a Wiener functional expansion. The extension to continuous systems is straightforward. Let x(n), n=0,1,2, ... be a white Gaussian input with mean zero and variance unity. Hen ce the time average autocorrelation is
(10)
(s+a)
x(n) ·x(n-k ) = o (k)
Analogous to (5), neglecting arbitrary constant s, G(sl)G(s:z) G(sl+s2 ) (11 ) (sl+s2+a)~xxp(Sl,S2)
by (10). We now take the limit of each side of (11) as s2 .... 00
where the delta function o( k)=O for k#O, and 0 (0)=l. Let yen), n=0,1,2, .. . be the system output, and y the corres ponding time average . We first find a cascade of a linear system with delta response hl (n), and a static nonlinea:it y f , whose output zl approximates l y(n)-y. Assume that the first order cross correlation ~xy
(k) = y(n)x( n-k)
is not identicallY zero. hl(n) = s2 • lim s +s s2-+OO 1 2
~
Then set (13 )
xy (n)
Define
Yl (n) = yen) - )'
lim s 2G(s2) s2 .... 00
-
Um (sl+s2)G(s l+ s2) s 2-+00
lim get)
-z
L""
(14 )
-2
Clearly y =0 and Yl=y -y. We wish to represent tAe static nonlincarity f by a weighted sum of Hermite polynomials H.{ · }. To this end, choose a , i=1,2, ... tolminimize li
t .... O
G(s l) lim get) t .... O (12)
The first fel" Hermite polyr!Omials are ~1
It was assumed here that get) contains no impulses at t=O and t'13t get) is not identical l y zero immediately to the ri ght of t=O . However, even in these circumstances, G(s l ) is the limit of the left side of (11). Hen ce, taking the l im it of the right side of
Hr. f l'
.
~
(16)
hI (k)x(n-k)
k=O
l~
( 1-1
H/· 1 =
'2 M hl(k)X(n-k) \ - L hiCk) k=O ) k=O ( 17)
( 11) , 1 im (s l+ s~+a)'l> _ xxp (sl , s2) 52 -+oc
(18 ) ~I
3
. t
xxp
(s
Therefore
s)
l' 2
..,
)' hi Ck ) k=O
~1
L
hl(k)x(n-k)
k=O
\ote that since Yl=O, the zero order Hermite polynomial HO ~as omitted from expression ( 15 ). If only the first P terms are significant in the Wiener expansion for the nonlinear system \, the upper limit for i in (15) may be set equal to P. Since the a 1 . minimize (15) and the H.{· } are orthogonal t6 one another, 1
p(t)x(t) ·x (t-O l )
M. Korenberg
672
Yl u
li
(n)Hi{k~O hI (k)X(n-k)} (19)
H~{k~O
In particular, since
~
(k)=~
1
(k)=hl(k), it
Once the al' are determined, the static nonlinearity 1fl is identified. The output of the cascade is
i~l aliHi{k~n hl(k)X(n-k)} (20) and the resiCl.ual nnt anproxil'lated by the cascade is (21 )
Y-;z(n) = yl(n) - zl(n)
equal to the mean-square of the first term on the right side of (20). From this it follows eas i ly that (22) Finally, (21) implies ~xy
As noted,
2
(k) = ~xy (k) - ~ 1
~
xz 1
2
l k=O
h2 (k)X(n-k)}J
The H.{o} are as illustrated in (17) and (18) except that hI in !hese eouations is replaced now by h 2 · Slnce Y2=0 and ~xy =0, the zero
and first order Hermite POlyno~ials HO' Hl{o} are omitted in (24). Again, if only tfie first P terms are significant in the Wiener expansion for the nonlinear system N, the upper limit for i in (2 d ) may be set eoual to P. The formula for u · is analogous 2 to (19), with Yl and hI replaced respectively by Y and h . 2 2
The output of the second cascade is z2(n)
""2
is easy to show Y2 = Yl - zl' so that the meansquare of the residual is decreased the greater the mean-souare of the cascade output. Since the H {·} in (20) are orthogonal to one i another, therefore z~ is p:reater than or
1
L
(24 )
Usinp: (19), (20) and (21), it ""2""2
n
( 11
a2 · H . i
hI (k)X(n-k)}
xYl xY follows from (19) that ull=l.
Clearly )-2=0.
To represent f2 as a sum of Hermite polynomials, we find coefficients a 2i to minimize
L
i=2
f 11L k=O
a 2 .Il. l 1
1
h 2 ( k)X(n-k)}
(25) and the residual not approximated by the first two cascades is (26) -
""2
-2-
""2
Clearly Y3=0, and Y3 = Y2 - z7.' It follows from (25) that ~ =0; therefore ~ =0. xZ 2 xY3 We now find a third cascade of a linear system with delta response h (n), and a static nonlinearity f3' whose ou~put z3 approximates the residual Y3' Set (27)
(k)
(k)=hl(k), and because x(n-k)
xYl is orthogonal to all terms on the right side of(20) except the first,
ensuring that h ~O (or we would change the sign of 6 in (~7) and begin again, reversing the sign of 6 in the formula for h 4 (n) and so on).
-~-,-- --.
a ll
.L
hl(j)x(n-j)·x(n-k)
J=O
= hI (k)
Note that the sign of 6 in (23) and (27) is reversed. In general, the sign of 6 alternates according to the formula
,
h . (n) = Hence t!J
J
=0 .
xY2 I';e next find another cascade of a linear system with delta response h (n), and a static 2 nonlinearity f2' whose output approximates the residual Y . Since ~ (n)=O, we must 2 xY2 now define h (n) from higher order cross 2 correlations. Assume that for some m, the second order cross correlation
Zz
~xxy
(k,m) = Y2(n)x(n-k)x(n-m) 2
is not zero for at least one k. h2 (n) =
~
xxY,}
Then set
(;t ,I'l) + 6 (n-m)
,
(23)
ensuring that h 2 (n) is not identically zero. (If hZ =O, we would change the sign of 6 in (23).)
~xxy. (n,m l J
+
(-1)j 6 (n-m), (28)
ensuring, as above, that hi(n) is not identically zero. To find f3' we choose coefficients u 3 · to minimize an expression analogous to l~4) with Y2' u 2i ' and h2 replaced respectivel y with Y3' u 3i ' and h 3 · The formula for u 3i is analogous to (19), with 1'1 and hI replaced respectively by y and h~. The output of the third cascade, Z3tn), is" defined analogously to (25), with u 2 · and h2 replaced respectively with a and fi . 3i 3 Let E be an arbitrary small positive number. We now show that if (29)
673
Parallel Cascades of Linear and Nonlinear Systems then for all k (k,T1) 1O. F lrst, Z2 IS rreater than or equal to the mean-square of the First term on the right side of (25), and after substituting for Y22 in this ter~, we obtain 1<1>
H;tk~n
~1
xxy (kl,m)
L k=(J
is not zero for at least one k. hK(n) = cP
cP xxy (m,m) 2
cP~Xy 2 (k ,m) I
~1
2
~1
) 1~= 0
M \' L
].;2=0
<
El
(32)
xxy (kl,m)cP xxv (k 2 ,m) 3 '3
•xxv (k l ,k 2 ) . 3
xxy (m,m)
+
xxxYK
+
where El depends on E ~nd tends to zero as E->O. Slmi larl y, from (27) and (29) ,
Ille
IfXXYK is ne~ligible, hK(n) should be defined instead from a hi~ler order cross correlation. Assume that for S0r'e ml and m2 , the third order cross correlation xxxy K(k ,m l ,m 2 ) = YK(n)x(n-k)x(n-m l )x(n-m 2 )
L
+
K
for all k and m. Therefore we can reduce the second order cross correlation of the residue arbitrarilY close to zern.
h 2 (k)X(n-k) }
+
"2
lxxy (k,m)1 < E' (31)
It is easy to show that the denominator is not zero and is les s than a bound which is independent of h, and depends solely on yz. Therefore, using~(73) and (29), it follows from the numerator of (31) that
xXY 2 (k 1 ,k Z)
2
zK + zK +1 < E Hence, in view of the above theorem,
~2HztL h 2 (k)x(n-k) ~2
-~ ~ ~=~;=======':'"" ( M ,
h (n) from (28), and so on. It is clear that 6 each ~ calculated reduces the mean-square of the re~idue, which can be reduced at most to zero. Therefore, we will eventually reach a cascade, say the Kth cascade where, no matter what value of m is used in (28) to define hK'
3
- 2
cP k=O
xXY3
(k,m) 1 <
E2
(33)
wher e E, depends on E and tends to zero as [ ->0. It folloh's from (26) that -cP = xXY xxy 2
1
XX\' .2
(k
'
m) \ I
I,h ere [ _ depends on E and tends to zero as [ ->0. Tfiis completes the nroof. The theorem has the folloh'ing use. Assume that we have adopted a very small positive value for E. If (29) is satisfied, then when we calculate the fourth cascade, we change the value of m in defining h 4 (n) from the gene ral formula (28). Again, if - 2
=~ +
z5
< E
we would change the va lue of m in defining
(34)
o(n-m?)
(r.,m l ,m 2 )+ 0 (n-m l ) xxxYK+I - 0 (n-m?) (35)
hK+2(n)=cPxXXYK+2 (n, ml,m2)-0(n-ml ) - 0( n-m 2 ;
(36)
hK 3(n)= (n,m l ,m 2 )- 0 (n-m l ) + xxxYK+3 +0( n-rn ) (37) 2
The signs of the 0 functions in the formula for ~K+4 are as in (34), and then the sequence of SIgn changes is repeated as just set out. It can be shown that if +Z2
K+2
+ --;:'K+3 2
(38)
where the delta responses are defined in (34) to (37), then for all k (39) (k,m 1 ,m 2 )1 < E" XXXYK ",here E" depends on E ancl. tends to zero as E+(J. 1<1>
-2
(n,r'l,r1 2 )+0(n-m l )
ensuring that hK(n) is not identically zero (by changing the sign of one of the 0 function s in (34) if necessary). In representing the static nonlinearity fK' the zero, first and second order Hermite polynomials are omitted since YK=(J, =0, cP ~O . xYK xXYK The procedure is continued analogously to that set out above for th~ second and third cascades. In defining t~le delta responses for the linear systems in cascades K+l, K+2, and K+3, the sign of one of t~le 0 functions in (34) is reversed each time as follows: hK+l(n)=cP
2
Then set
M. Korenberg
674 FUNCTIONAL EXPANSIONS AND DIFFERENTIAL EQUATION EXPANSIONS
Once the above identification is complete, the system output (actually y-y) can be expressed using the sum of the cascade outputs, and by collecting Hermite polynomials of like degree, we immediately obtain the Wiener functional expansion for y in terms of orthogonal functionals. Alternately; each of the static nonlinearities can be rearranged as an ordinary polynomial or power series, so that the output of the jth
casca:~(:~ :quyva:~~t[l~ J
i=O
J 1 k=O
h.(k)X(n_k)]i J
( 40) Usinr this representation for z., summing over all the cascades and collecting like powers, ~e immediately obtain the Volterra functional expansion for y. A method of developing a differential equation expansion for the nonlinear system will now be outlined. Although a discrete representation was obtained for the nonlinear system, it will be assumed below that the real system has an output yet) defined in continuous time. For simplicity, it will be supposed that y=O; otherwise the mean would be subtracted from the output. Let the input to the nonlinear system be x(t) = A sin wt
(41 )
From the above identified representation, we can immediately write down the corresponding output as the sum of the individual cascade outputs: 00
yet)
L
Bli(A'H l Sin( wt+ 8 1 ))i
L
B2i (A'H 2 sin(wt+8 2 ))i+ ...
\" L
A si(t)
i=O +
i=O
i
(42)
i=O
output derivatives respectively. Substituting (41) and (42) into (44) and equating the coefficients of the linear terms in A yields:
I
cisii)(t) =
i=O
~
i=O
b. sin(i)wt (45) 1
Since sl(t) is known from (43), the parameters c. and b. (and the order of the differenfial eqGation) may be estimated from Bode plots, or from linear regression on (45) using various values of w. Similarly, substituting (41) and (42) i~to (44) and equating the coefficients of A- would yield an equation enabling us to estimate the parameters c .. and b... The re~aining parameters a}J estimited rulalogously so that expansion (44) may be developed to an arbitrary extent. CONCLUSION The general identification procedure proposed here was illustrated using successive approximations of a simple cascade model. In a future paper, the proc~dure will be implemented when the basic model is a simple non linear differential equation rather than the cascade model employed above. REFERENCES Billings, S .A. (1980). Proc. lEE, llZ., Pt. D. 272 -28 5. Billings, S.A., and S.Y. Fakhouri (1978). Proc. lEE, 125, 69l-6~7. Korenberg, H. (1973a). Proc. 16th Hidwest Symp. Circuit Theory, 18.2, 1-9. Korenberg, ~1. (1973b) . Proc. 10th Ann. Rocky Hountain Bioeng. S~P .. 47 -52. Lamme~H.C . and E. De"Boer (1979). Proc. IEEE, 67, 432-434. Harmarelis,P.Z., and V. Z. ~larJTlarelis (1978). Analysis of Ph ys iological Systems - The White-Noise Approach . Plenum Press, New York. Spekreijse, H . (1969). Vision Res. ~, 14611472 . Spekreijse, H., and H. Oosting (1970). Kybernetik, 7 , 22 - 31.
where 00
Si (t) =
L
j=O
S·· (H. Sin( wt+ 8 .))i (43) J Jl J
Here the B. . , H. and 8. are known, the coefficients B~~ de~ining the static nonlin earitv in the j-t~lcascade, and the H. and 8. de- . fining the gain and phase char~cteristics of the linear system in this cascade. The nonlinear system may be approximated to an arbitrary accuracy using a differential equation expansion of the form:
I
i=O
c.y(i) + 1
I I
c .. y(i)/j) + ... i=O j=O IJ
~
b.x(i)+
~
~
b .. x(i)x(j)+ .. . IJ . (44) where xCi) and y(l) are the i-th input and i=O
1
i=O j=O
k(~)l-~ ~-~t~~ ~l ~t~:_Plt2. ~ '--.J
1_ _ _
:.r__
I
g (C1) =p=-CCr::tT)X""'C~tT).,...x'C7t -: _oC"C)
Fig. 1: Cascade of a linear, a static nonlinear, and a linear system. The dotted portion illustrates feeding the output of the cascade through an additional linear filter in order to identify get).
STATISTICAL IDENTIFICATION OF PARALLEL CASCADES OF LINEAR AND NON LINEAR SYSTEMS M. Korenberg R ogers, Bereskin and Parr, P. O. Box 313, Commerce Court Wes t, Suit e 3710, T oronto, On tario, Canada
Abstract. An orthogonal statistical method is proposed for identifying a nonlinear s\'stem by successively adding a basic model paralle l to si~i l ar such models. The basic model selected must not be closed under acldition. A first such model is obtained to approximate the output of the non l inear s\'stem . The difference between the system output and the first mode l output is treated as a new output, and a second model is obtained to approxirate the l atter output . From this output the second model ou t put is subtracted to produce a third output, and a third model is obtained to approximate t'lis output, and so on. The method is illustrated when the nonlinear system is assumed to have a \\'iener functional expansion, and the basic model is a cas cade of a linear system and a static nonlinearity. It is shown that the nonlinear system may be approximated to an arbitrary accuracy by the sum of the models obtained . and that a functiona l expansion or a differentia l c~uation expansion may then be wri t ten down for the nonlinear system . An explicit time domain so l ution is a l so proposed for identifying a cascade of a linear. a static non linear and a linear system, which gives the impulse rpsronses direct l y by specia l first order cross corre l ations. Keywords. Nonlinear systems; identification; models; correlation net~ods; discrete systems; parameter estimat ion; time-domain analysis; cascacles. I\TRODUCTIO\ A typical approach to identifying a nonlinear system is to approximate its outpu t by means of some basic model, for e xample a Vo l terra functional expansion. Once the best such expansion is obtained, in the least-square sense, the repre sentation cannot be improved by adding a second such expan s ion paral l e l to the first, because the sum of an\' two Vo l terra series is a Volterra series. The situation is different, however, when the basic model chosen is not closed under addition. for example, suppose t hat the basic mode l -:omDrises a linear s;'stem follOl,ed by a static nonlinearit\' in the form of a pol:~omial or ;-'Oher series, and that the "best" impulse response and coefficients def i ning the nonlinearit\· are found to approximate the nonlinear system . The residual, i.e. the dif ference betl,een the svstel'1 output and the model output, can itse l f be treated as a neh output, and a second mode l of a l inear system followed b\' a static nonlinearity can be ob tained to approximate the neh output. The sum of the tlW model O:ltputS I,ill be c l oser to the system output, in the mean square sense, than is the first mode l ou t put a l one. The difference betheen the system output and the sum of the thO model outputs may be treated as a third output, and a third such
669
model obtained to approximate this output. The procedure may be repeated indefinitely to approximate the nonlinear system to an arbitrary accuracy in the mean square sense by a sum of the mode l s obtained separately. This procedure is predicated on the fact that the sum of two such models is not equi\'alent to a sing l e mode l of the same type. In other hords, it is necessary (but not sufficient) that the basic model selected not be closed under addition. An advantage of the proce dure is that, while the final representation comprises a parallel array of the models, at each stage of t he identification on l y a sing l e mode l is computed. ~loreover the procedure is orthogonal in tha t an additional mode l may be identified to improve the appro ximation hithout recalcu l ating the previous models obtained. Although there are many possible choices for the basic mode l , the procedure wi l l be illustrated here I,hen the mode l is a cascade of a l inear system fol l ohed by a static nonl·inearity. BOI,'ever it I,'ill be appreciated that the me th od proposed is not limited to successive approximat ions by the cascade model, and that other models not closed under addi t ion may be substi t uted.
M. Korenberg
670
Before illustrating the proposed method, however, the identification of a cascade of a linear, a static nonlinear, and a linear system (Spekreijse, 1969; Spekreijseand Oosting, 1970 ; Korenberg, 1973a, 1973b) is re-examined. It is shown that the i~pulse response of each linear system can be obt~ined directly as a special first order cross correlation. Next the proposed method is illustrated when the basic model comprises a linear system followed by a static nonlinearity, the extension to a basic model having a second linear system following the nonlinearity being trivial. Finally, it is shown that a functional expansion or R differential equation expansion can be obtained immediately from t~e parallel casc ade r epresentation identified.
A similar formula for IK(s) I then follows from (3). Analogous formulas exist for the phase of G(s) and K(s). Thus the linear systems may Qe identified from the first and second order cross correlations (and the static nonlinearity then determined). Alternately, setting sl=-s2=s in (4) and neglecting the constant IG(s) I =
2[Y le~::_S) 'j ~ xx~(O)__
However the latter formula for IG(s) I is less general than (6) and cannot be used if K(O)=O, which could occur, for example, if the second linear system in Fig. 1 involves a differentiator. y i e !~s
Setting ° =° =0 in (2) 1 2 4>XXY(O,O) = A2
CASCADE OF A LINEAR, A STATIC NONLINEAR, AND A LINEAR SYSTEM Spekreijse (1969) suggested an interesting nethod for identifying a cascade of a linear, a static nonlinear, and a linear system (Fig. 1) using special test inputs and re~eated experiments. Korenberg (1973a, 1973b) nroposed a correlation technique for explicitly identifying each corronent in the cascade following a single "aussian white-noise experiment. The latter technique involves evaluating the time-averaged first order cross correlation y(t)x(t-o) Alfoo k(T)g(o- T)dT
o
(1)
and the time-averaged ' second order cross correlation, shown to equal 4>XX/ O l, 0 2 ) = A2 foo k( T) g(01-T)g(02-T)dT ,
o
(2)
"'here x(t) is a \yhite "aussian input, get) and k(t) are the impulse responses of the first and second linear systems r espectivel y , Al and A2 are constants, yet) is the cascade output, and the mean of the output is subtracted before the second order cross correlat ion is comput ed (c.f . l1armarelis and Harmarelis, 1978). Taking the Laplace trans:Oorm of (1) and the two-dimensional Laplace transform of (2) yields respectively ~
and Hence
xy (s) = AlK(s )r;(s)
(3)
Al xxv(s l,s2 ) =:\2
(5)
and. letting sl=-s/2, s2=s, and neglecting arbltrarv constants . 1<1> (-s/2,s)1 IG(s) I = xxv ( 6)
1~y(s/2)I
(7)
J:
l:(T)i(o-T) dT
(8)
Billings and Fakhouri (lQ7U) have shown elegantly that by expressinp, (1) and (8) in sampled data form and assuming the sampled correlation functions am', the linear systems can be fit by rational pu1,sP' transfer functions . the linear systems CRn be identified. The cascade of Fig . 1 has indeed been extensively studied, e.g. Lammers and De Boer (1979), Bi 11 ings (1980). When no assumptions are made about the form of the linear systems, usin p 4> x (0,0) in place of 4> (0 ,0 2 ) may lose f06 much in1 formation f3 Yunlquely identify the linear systems. For example, if pet) and k(t) each assume only the values 0 and I but are different functions, then using (1) and (8) it will not be possible to det ermi ne the order in which the linear systems appear in the cascade. This is because g2(t)=g(t) and k 2 (t)=k(t), so that g anc, I: may be interchanged in (1) and in (8). This problem never arises using (1) and (2), or e quivalentl y (3) and ( 4), so that get) and k(t) can then be uniquely identified up to constants. It is now shown that get) <'.nd k(t) may be directl y determined in the time domain from 4> (0) and a special cross correlation. Let
4>~r(t) be the inverse of 4> (t), so that the c3Xvolution of these two fGrictions is the unit impulse of Dirac. Form the linear system L whose impulse response is the conv01u-at -1 tion of e and 4> (t) (for some a>O), i.e. xy e- at * 4> -l(t) xy and feed the output y of the cascade of Fig. 1 through L. Let pet) be the corresponding output of L. It will be proved that g(o) is given by the time average p(t)x(t)'x(t- o )
,°
(9)
Let 4> xp (0) and 4> xxp (° 1 2 ) be respectively the first and second order cross correlations of x and p defined in the usual manner. Note that the linear system L will have transfer
Parallel Cascades of Linear and Nonlinear Systems
671
GENERAL NONLINEAR IDENTIFICATION
function (s+a)1>
(s) xy Therefore, analogous to (3), neglecting arbitrary constants.
For convenience, the identification method will be demonstrated for a discrete nonlinear system N having finite memory of length M and a Wiener functional expansion. The extension to continuous systems is straightforward. Let x(n), n=0,1,2, ... be a white Gaussian input with mean zero and variance unity. Hen ce the time average autocorrelation is
(10)
(s+a)
x(n) ·x(n-k ) = o (k)
Analogous to (5), neglecting arbitrary constant s, G(sl)G(s:z) G(sl+s2 ) (11 ) (sl+s2+a)~xxp(Sl,S2)
by (10). We now take the limit of each side of (11) as s2 .... 00
where the delta function o( k)=O for k#O, and 0 (0)=l. Let yen), n=0,1,2, .. . be the system output, and y the corres ponding time average . We first find a cascade of a linear system with delta response hl (n), and a static nonlinea:it y f , whose output zl approximates l y(n)-y. Assume that the first order cross correlation ~xy
(k) = y(n)x( n-k)
is not identicallY zero. hl(n) = s2 • lim s +s s2-+OO 1 2
~
Then set (13 )
xy (n)
Define
Yl (n) = yen) - )'
lim s 2G(s2) s2 .... 00
-
Um (sl+s2)G(s l+ s2) s 2-+00
lim get)
-z
L""
(14 )
-2
Clearly y =0 and Yl=y -y. We wish to represent tAe static nonlincarity f by a weighted sum of Hermite polynomials H.{ · }. To this end, choose a , i=1,2, ... tolminimize li
t .... O
G(s l) lim get) t .... O (12)
The first fel" Hermite polyr!Omials are ~1
It was assumed here that get) contains no impulses at t=O and t'13t get) is not identical l y zero immediately to the ri ght of t=O . However, even in these circumstances, G(s l ) is the limit of the left side of (11). Hen ce, taking the l im it of the right side of
Hr. f l'
.
~
(16)
hI (k)x(n-k)
k=O
l~
( 1-1
H/· 1 =
'2 M hl(k)X(n-k) \ - L hiCk) k=O ) k=O ( 17)
( 11) , 1 im (s l+ s~+a)'l> _ xxp (sl , s2) 52 -+oc
(18 ) ~I
3
. t
xxp
(s
Therefore
s)
l' 2
..,
)' hi Ck ) k=O
~1
L
hl(k)x(n-k)
k=O
\ote that since Yl=O, the zero order Hermite polynomial HO ~as omitted from expression ( 15 ). If only the first P terms are significant in the Wiener expansion for the nonlinear system \, the upper limit for i in (15) may be set equal to P. Since the a 1 . minimize (15) and the H.{· } are orthogonal t6 one another, 1
p(t)x(t) ·x (t-O l )
M. Korenberg
672
Yl u
li
(n)Hi{k~O hI (k)X(n-k)} (19)
H~{k~O
In particular, since
~
(k)=~
1
(k)=hl(k), it
Once the al' are determined, the static nonlinearity 1fl is identified. The output of the cascade is
i~l aliHi{k~n hl(k)X(n-k)} (20) and the resiCl.ual nnt anproxil'lated by the cascade is (21 )
Y-;z(n) = yl(n) - zl(n)
equal to the mean-square of the first term on the right side of (20). From this it follows eas i ly that (22) Finally, (21) implies ~xy
As noted,
2
(k) = ~xy (k) - ~ 1
~
xz 1
2
l k=O
h2 (k)X(n-k)}J
The H.{o} are as illustrated in (17) and (18) except that hI in !hese eouations is replaced now by h 2 · Slnce Y2=0 and ~xy =0, the zero
and first order Hermite POlyno~ials HO' Hl{o} are omitted in (24). Again, if only tfie first P terms are significant in the Wiener expansion for the nonlinear system N, the upper limit for i in (2 d ) may be set eoual to P. The formula for u · is analogous 2 to (19), with Yl and hI replaced respectively by Y and h . 2 2
The output of the second cascade is z2(n)
""2
is easy to show Y2 = Yl - zl' so that the meansquare of the residual is decreased the greater the mean-souare of the cascade output. Since the H {·} in (20) are orthogonal to one i another, therefore z~ is p:reater than or
1
L
(24 )
Usinp: (19), (20) and (21), it ""2""2
n
( 11
a2 · H . i
hI (k)X(n-k)}
xYl xY follows from (19) that ull=l.
Clearly )-2=0.
To represent f2 as a sum of Hermite polynomials, we find coefficients a 2i to minimize
L
i=2
f 11L k=O
a 2 .Il. l 1
1
h 2 ( k)X(n-k)}
(25) and the residual not approximated by the first two cascades is (26) -
""2
-2-
""2
Clearly Y3=0, and Y3 = Y2 - z7.' It follows from (25) that ~ =0; therefore ~ =0. xZ 2 xY3 We now find a third cascade of a linear system with delta response h (n), and a static nonlinearity f3' whose ou~put z3 approximates the residual Y3' Set (27)
(k)
(k)=hl(k), and because x(n-k)
xYl is orthogonal to all terms on the right side of(20) except the first,
ensuring that h ~O (or we would change the sign of 6 in (~7) and begin again, reversing the sign of 6 in the formula for h 4 (n) and so on).
-~-,-- --.
a ll
.L
hl(j)x(n-j)·x(n-k)
J=O
= hI (k)
Note that the sign of 6 in (23) and (27) is reversed. In general, the sign of 6 alternates according to the formula
,
h . (n) = Hence t!J
J
=0 .
xY2 I';e next find another cascade of a linear system with delta response h (n), and a static 2 nonlinearity f2' whose output approximates the residual Y . Since ~ (n)=O, we must 2 xY2 now define h (n) from higher order cross 2 correlations. Assume that for some m, the second order cross correlation
Zz
~xxy
(k,m) = Y2(n)x(n-k)x(n-m) 2
is not zero for at least one k. h2 (n) =
~
xxY,}
Then set
(;t ,I'l) + 6 (n-m)
,
(23)
ensuring that h 2 (n) is not identically zero. (If hZ =O, we would change the sign of 6 in (23).)
~xxy. (n,m l J
+
(-1)j 6 (n-m), (28)
ensuring, as above, that hi(n) is not identically zero. To find f3' we choose coefficients u 3 · to minimize an expression analogous to l~4) with Y2' u 2i ' and h2 replaced respectivel y with Y3' u 3i ' and h 3 · The formula for u 3i is analogous to (19), with 1'1 and hI replaced respectively by y and h~. The output of the third cascade, Z3tn), is" defined analogously to (25), with u 2 · and h2 replaced respectively with a and fi . 3i 3 Let E be an arbitrary small positive number. We now show that if (29)
673
Parallel Cascades of Linear and Nonlinear Systems then for all k (k,T1) 1
H;tk~n
~1
L k=(J
is not zero for at least one k. hK(n) = cP
cP xxy (m,m) 2
cP~Xy 2 (k ,m) I
~1
2
~1
) 1~= 0
M \' L
].;2=0
<
El
(32)
•
+
xxxYK
+
where El depends on E ~nd tends to zero as E->O. Slmi larl y, from (27) and (29) ,
Ille
If
L
+
K
for all k and m. Therefore we can reduce the second order cross correlation of the residue arbitrarilY close to zern.
h 2 (k)X(n-k) }
+
"2
lxxy (k,m)1 < E' (31)
It is easy to show that the denominator is not zero and is les s than a bound which is independent of h, and depends solely on yz. Therefore, using~(73) and (29), it follows from the numerator of (31) that
2
zK + zK +1 < E Hence, in view of the above theorem,
~2HztL h 2 (k)x(n-k) ~2
-~ ~ ~=~;=======':'"" ( M ,
h (n) from (28), and so on. It is clear that 6 each ~ calculated reduces the mean-square of the re~idue, which can be reduced at most to zero. Therefore, we will eventually reach a cascade, say the Kth cascade where, no matter what value of m is used in (28) to define hK'
3
- 2
cP k=O
xXY3
(k,m) 1 <
E2
(33)
wher e E, depends on E and tends to zero as [ ->0. It folloh's from (26) that -cP = xXY xxy 2
O in view o~ xx. -, (29).
Therefor e (32) and (33) imply
I ~1,
\
I, k=O
1
XX\' .2
(k
'
m) \ I
I,h ere [ _ depends on E and tends to zero as [ ->0. Tfiis completes the nroof. The theorem has the folloh'ing use. Assume that we have adopted a very small positive value for E. If (29) is satisfied, then when we calculate the fourth cascade, we change the value of m in defining h 4 (n) from the gene ral formula (28). Again, if - 2
=~ +
z5
< E
we would change the va lue of m in defining
(34)
o(n-m?)
(r.,m l ,m 2 )+ 0 (n-m l ) xxxYK+I - 0 (n-m?) (35)
hK+2(n)=cPxXXYK+2 (n, ml,m2)-0(n-ml ) - 0( n-m 2 ;
(36)
hK 3(n)= (n,m l ,m 2 )- 0 (n-m l ) + xxxYK+3 +0( n-rn ) (37) 2
The signs of the 0 functions in the formula for ~K+4 are as in (34), and then the sequence of SIgn changes is repeated as just set out. It can be shown that if +Z2
K+2
+ --;:'K+3 2
(38)
where the delta responses are defined in (34) to (37), then for all k (39) (k,m 1 ,m 2 )1 < E" XXXYK ",here E" depends on E ancl. tends to zero as E+(J. 1<1>
-2
(n,r'l,r1 2 )+0(n-m l )
ensuring that hK(n) is not identically zero (by changing the sign of one of the 0 function s in (34) if necessary). In representing the static nonlinearity fK' the zero, first and second order Hermite polynomials are omitted since YK=(J, =0, cP ~O . xYK xXYK The procedure is continued analogously to that set out above for th~ second and third cascades. In defining t~le delta responses for the linear systems in cascades K+l, K+2, and K+3, the sign of one of t~le 0 functions in (34) is reversed each time as follows: hK+l(n)=cP
2
Then set
M. Korenberg
674 FUNCTIONAL EXPANSIONS AND DIFFERENTIAL EQUATION EXPANSIONS
Once the above identification is complete, the system output (actually y-y) can be expressed using the sum of the cascade outputs, and by collecting Hermite polynomials of like degree, we immediately obtain the Wiener functional expansion for y in terms of orthogonal functionals. Alternately; each of the static nonlinearities can be rearranged as an ordinary polynomial or power series, so that the output of the jth
casca:~(:~ :quyva:~~t[l~ J
i=O
J 1 k=O
h.(k)X(n_k)]i J
( 40) Usinr this representation for z., summing over all the cascades and collecting like powers, ~e immediately obtain the Volterra functional expansion for y. A method of developing a differential equation expansion for the nonlinear system will now be outlined. Although a discrete representation was obtained for the nonlinear system, it will be assumed below that the real system has an output yet) defined in continuous time. For simplicity, it will be supposed that y=O; otherwise the mean would be subtracted from the output. Let the input to the nonlinear system be x(t) = A sin wt
(41 )
From the above identified representation, we can immediately write down the corresponding output as the sum of the individual cascade outputs: 00
yet)
L
Bli(A'H l Sin( wt+ 8 1 ))i
L
B2i (A'H 2 sin(wt+8 2 ))i+ ...
\" L
A si(t)
i=O +
i=O
i
(42)
i=O
output derivatives respectively. Substituting (41) and (42) into (44) and equating the coefficients of the linear terms in A yields:
I
cisii)(t) =
i=O
~
i=O
b. sin(i)wt (45) 1
Since sl(t) is known from (43), the parameters c. and b. (and the order of the differenfial eqGation) may be estimated from Bode plots, or from linear regression on (45) using various values of w. Similarly, substituting (41) and (42) i~to (44) and equating the coefficients of A- would yield an equation enabling us to estimate the parameters c .. and b... The re~aining parameters a}J estimited rulalogously so that expansion (44) may be developed to an arbitrary extent. CONCLUSION The general identification procedure proposed here was illustrated using successive approximations of a simple cascade model. In a future paper, the proc~dure will be implemented when the basic model is a simple non linear differential equation rather than the cascade model employed above. REFERENCES Billings, S .A. (1980). Proc. lEE, llZ., Pt. D. 272 -28 5. Billings, S.A., and S.Y. Fakhouri (1978). Proc. lEE, 125, 69l-6~7. Korenberg, H. (1973a). Proc. 16th Hidwest Symp. Circuit Theory, 18.2, 1-9. Korenberg, ~1. (1973b) . Proc. 10th Ann. Rocky Hountain Bioeng. S~P .. 47 -52. Lamme~H.C . and E. De"Boer (1979). Proc. IEEE, 67, 432-434. Harmarelis,P.Z., and V. Z. ~larJTlarelis (1978). Analysis of Ph ys iological Systems - The White-Noise Approach . Plenum Press, New York. Spekreijse, H . (1969). Vision Res. ~, 14611472 . Spekreijse, H., and H. Oosting (1970). Kybernetik, 7 , 22 - 31.
where 00
Si (t) =
L
j=O
S·· (H. Sin( wt+ 8 .))i (43) J Jl J
Here the B. . , H. and 8. are known, the coefficients B~~ de~ining the static nonlin earitv in the j-t~lcascade, and the H. and 8. de- . fining the gain and phase char~cteristics of the linear system in this cascade. The nonlinear system may be approximated to an arbitrary accuracy using a differential equation expansion of the form:
I
i=O
c.y(i) + 1
I I
c .. y(i)/j) + ... i=O j=O IJ
~
b.x(i)+
~
~
b .. x(i)x(j)+ .. . IJ . (44) where xCi) and y(l) are the i-th input and i=O
1
i=O j=O
k(~)l-~ ~-~t~~ ~l ~t~:_Plt2. ~ '--.J
1_ _ _
:.r__
I
g (C1) =p=-CCr::tT)X""'C~tT).,...x'C7t -: _oC"C)
Fig. 1: Cascade of a linear, a static nonlinear, and a linear system. The dotted portion illustrates feeding the output of the cascade through an additional linear filter in order to identify get).