MarhI Compur. Modeliing, Vol 12, No. 8, pp. 991-995, Printed in Great Britain. All rights reserved
0895-7177189 93.00 + 0.00 Copyright 0 1989 Maxwell Pergamon Macmillan plc
1989
IDENTIFICATION AND MODELLING OF A CLASS OF NONLINEAR SYSTEMS V. Z. Departments
of Biomedical
MARMARELIS
and Electrical Engineering, University CA 90089-1451. U.S.A.
of Southern
California,
Los Angeles,
(Received February 1988; accepted for publication November 1988) Communicated
by X. J. R. Avula
Abstract-This paper presents an identification procedure based on Volterra-Wiener functional expansion representations of a class of nonlinear systems often encountered in practice. Approximate analytical expressions are derived for the Volterra and Wiener kernels of these systems in the case of weak nonlinearities, i.e. for small magnitude coefficients of the polynomial nonlinearity. This identification method is based on Wierner kernel measurements obtained from input-output data when the input is Gaussian white noise. This method can be used for modelling and validation studies of this class of systems from experimental input-output data.
1.
INTRODUCTION
Consider the class of nonlinear systems that can be described by the ordinary nonlinear differential equation: L[y] + f cc;y’= x
(1)
r=2
where L[ .] is a linear time-invariant differential operator, x is the input and y the output of the system. We will examine the case where 1clil < c < 1. For the region of stable solutions of this equation, there exists a Volterra functional expansion [l-4]: y(t) =
$ ...k,(r,,
. . . ,z,)x(t
-r,),
. .x(t
-r,)dT,
. . .dr,
(2)
which represents the system output in terms of a series of multiple convolution integrals of the input. The kernel functions {k,) characterize the dynamics of the nonlinear system and are called the Volterra kernels of the system. They are symmetric functions of their arguments, i.e attain the same value for all permutations of given (5,) . . . , 5,) values. A method of “generalized harmonic balance” can be used to derive analytically the Volterra kernels that correspond to the system of equation (1). This method is based on ideas found in [l, 51. According to this method, the mth order Volterra kernel can be evaluated in the m-dimensional Laplace domain by considering generalized harmonic inputs: x,(t) = exp(s, 2) + . . . + exp(s,t)
(3)
where s, are distinct complex (Laplace) variables. Substitution of the input x,(t) into equation (2) results in the output expression: y,,(t)=
f
i
t7=1J,=i
x
...
f j, =
exp( -s,,T, - . . . -
s ~~~7,)
k,(~,,...,5,)
...
exp[(sj,+...+s,n)t]
1
s
ds, . . . dT,,
= n$, j!, . ’ ‘j$, exp[@j, + ’ ’ ’ + sjn)rlKn(sj! 2. *. 3sjn> where K, is the n-dimensional
Laplace transform of the kernel k,. 991
(4)
992
V.Z.
MARMARELIS
Note that the output terms that contain the complex exponential with all distinct complex frequencies of the input (i.e. terms of the form exp[(s, + . . . + sm)r]) are terms associated with the mth order kernel. Substitution of the output expression of equation (4) into equation (1) yields an equation that can be used for the evaluation of the m th order kernel K,,(s,, . . . . , s,) on the basis of harmonic balance, i.e. by selecting only those terms in the equation that contain the complex exponential exp[(s, + . . . + s,)t]. The contributions of these terms in the equation must balance out, according to the harmonic balance approach. Let us follow this approach for increasing values of m. For m = 1 we have: (5) and A(t) = .zj . . .n$, expb, Therefore,
+ . . + nib, flK, b,
. . . , sI 1.. .K, Gl,
. , %I.
(6)
balancing terms with exp(s, t) we obtain: Uexp(s, t)K(s, >I= eW, t).
(7)
If L is a Ith order differential operator:
4.1
dC.1
L[.]=c,+c,~+-+C,~
(8)
then: K(s) =
1 c, + c,s + . . . + c,s’.
(9)
Note that the power terms in equation (1) do not contribute to K,, i.e. the first order Volterra kernel represents strictly the linear portion of the nonlinear differential equation (1). The second order kernel can be evaluated for m = 2. Then: y*(t)=
f ?f=lJ,=l
i
.
..$ eXp[(Sj,+*.
’+
sjn
In= 1
)flKn(sjg 3. . . 9sj, 1.
(10)
Since the resulting expressions for y:(t) are rather unwieldy, we focus only on those terms that will contain the exponential exp[(s, + sz)t]. No such terms will be present in y;(t) for i > 2. The only such terms will come from y:(t) and, of course, from yz( t) itself. Observing that the kernels are symmetric functions of their arguments, we have: L (2 exp[(s, + ~d~l&(sI1~~1)+ 2Gh
(11)
VG(~devK~,+ ~~)tl= 0.
Solving equation (11) for K2, we obtain: K,(s, 3~2) =
-~~KI(~,K(J~W,(J,
(12)
+ ~2).
Continuing on with m = 3, we see that terms containing the exponential exp[(s, + s2 + s3)t] are contributed by y3(t),y:(t) and y:(t). The resulting equation is: L{3! exp[(s, + s2 + ~&MS,, s2,s3)) + 2~~w-G, +x2 + s,)~llK,h)K2(~2,s3) + 4 ts2K2C%,sl>
+K,(s,)K2(~,r~2)l+3!~3exp[(s,+s,+~~)~lK,(~,)K,(~z)K,(s3)=0. (13) Solving equation (13) for K3 and substituting K2 from equation (12) we obtain: 4 K~(S,~S~,~~)=~K,(S,)K,(S~)K~(S))~K,(S~+S))+K,(S+S,)+K,(S,+S~)IK,(S,+S~+S))
- a~K,(s,)K,(q)K,(s~)K,(s,
+ s2 + 4
(14)
Since we have assumed that (ai/ < L < 1 for all i, the first term of the right-hand side of equation (14) may be considered negligible. We observe that this term constitutes the contribution of the
Nonlinear systems identification and modelling
993
square term of equation (1) to the third order kernel. This observation can be generalized by considering the mechanics of this harmonic balance approach to state that: the contributions to the mth order kernel of power terms of equation (1) with degree less than m, are of order (6’) or higher and, consequently, negligible under the assumption c 4 1. Observing also that power terms of equation (1) with degree greater than m do not contribute at all to the mth order kernel, we conclude that, under the stated assumption for the coefficients OL~, we have in first approximation (for m 3 2):
K,(s, , . . .,s,)=
-cc,K,(s,)...K,(s,)K,(K,(s,+...+s,)).
(15)
The form of the Volterra kernels given by equation (15) indicates that this class of systems can be approximated with “S, models” studied in [6,7].
2. THE
IDENTIFICATION
PROCEDURE
Since the Volterra kernels of a “black-box” nonlinear system cannot be directly evaluated from input-output data, Wiener proposed the orthogonalization of the Volterra functional series for a Gaussian white noise (GWN) input. The resulting orthogonal functional expansion (termed the Wiener series) allows the practical identification of the associated Wiener kernels by isolating one kernel at a time through the use of appropriate instrumental functionals. This approach has been studied extensively in the last 30 yr and the reader is referred to the original Wiener monograph [S] and recent books on the subject [2,3,93 for details. A summary, necessary for the developments of this paper, is given below. The Wiener orthogonal series is: y(t) = :
G,,[h,, x0’), t’ < t]
(16)
ll=O
where the n th order Wiener functional
x(t -r,)*
..x(t
-Tn_Zm)d5,...dz,_2,da,...do,
(17)
contains the n th order Wiener kernel h, and the power level P of the GWN input. The most widely used method for the estimation of the nth order Wiener kernel is the Lee-Schetzen cross-correlation technique [lo], according to which:
MT, . . . L) =
-&
E[y’“‘(t)x(t
- 5,) . - . x0 -
L)l
(18)
where It-1
y’“‘(t) = y(t) -
c
G,[h,;
x(t’,
I’ < t)].
(19)
m=O
As indicated by equations (18) and (19), the kernels are estimated successively, in ascending order. In practice the ensemble-average of equation (18) is evaluated as a time-average over the experimental input-output record under the assumption of system stationarity. Issues of actual implementation and estimation accuracy in this approach have been addressed extensively in [3,9]. The Wiener kernels are, in general, different from the Volterra kernels of a given system (for which both expansions exist.) The n th order Wiener kernel can be expressed in terms of the Volterra kernels of the same and higher order as [2, 91:
UT,,
’
.
.
3
k (z z (n+2m)!P" 00 . *. T,) = 1 n+2m I)...1 n!m!2" s0 ItI= s
t,,
Q,,u,,
. . . . t~,,,,o,,,)dtr,
* 9 * da,,, (20)
994
or in the frequency
V.Z. MARMARELIS
domain:
-u,,, For the specific than first have the with equation (21) Wiener kernels of
. . ,U,,-z4U,)du,.‘.du,
(21)
class of systems considered in Section 1, the Volterra kernels of order higher approximate form given in equation (15). Therefore, combining equation (15) we obtain, in first approximation, the following expression for the high-order this class of systems (n 2 2): H,(,f,,
. . .,fn)=
-~,K,(,f,)..,K,(f,)K,(f,$-...+f,)
(22)
where E = f ‘I?
(n + 2m)!P”Q”
??l=O
Q= For the first order
Wiener
kernel,
s=
n!m
!2”
CIn+ Zm
IK,(u)l’du.
(23)
(24)
- 7.
we have: H, (f) = K, 0
- A, K:(Y).
(25)
Equations (22) and (25) are pivotal in identifying the considered class of systems from input-output data. From the identification point of view, the Wiener kernels are estimated from It would be desirable to have input-output data, while K,(f) and (A,, &, . . , A,,,) are unknown. a way of estimating these unknown quantities by use of equations (22) and (25). The unknown coefficients (a,, . , ‘xM) could then be evaluated from the estimates of (&, , AM) using equation (23), leading to complete identification of the system in equation (1). There are many ways to utilize equations (22) and (25) in order to estimate K,(,f) from Wiener kernel measurements. Consider, for instance, the following slice of H,: H, (L 0) = - GK,
(O)K:(f).
Since H, is known (i.e. computed from input-output scalar y’ = l/(-&K, (0)), and substituted into equation H,(f)
= a~[H,(f,
data), K:(f) can be evaluated (35) to yield:
0)1”2- h~~ffz(f, 0).
(26)
within
a
(27)
The unknown quantities in equation (27) y and i, , can be estimated from measurements H, (f) and H,(f, 0) through least-squares fitting. Note that linear regression can yield estimates of 7 and (A,y2), from which II, can be obtained; and that two opposite values of y estimates result, corresponding to the two branches of [H,(f, O)]“‘. Having estimated y, we now obtain two possible solutions for K,(f) (fixed phase difference of 71 rad): K, Cr) = + “r’[%(f,
0)1”2.
(28)
Using the estimates of K,(f), we can obtain estimates of unknown scalars (A*, . . , A,,,) from equation (22). This, of course, is a burdensome task in practice, since it requires computation of high-order Wiener kernels. However, since the problem is highly overdetermined (i.e. only a single parameter need be estimated from each kernel equation) we can limit kernel computation to a small portion of its domain (in principle, only at one point.) Having estimated the scalars (AZ, . . . , A.,), we can use equation (23) to evaluate the unknown coeflicients (c(*, . . . , mH) through linear inversion. Note that the matrix that need be inverted is upper triangular. This completes the identification of the system in equation (1). The outlined procedure for the estimation of K, (f) from measurements of H,(f, 0) is not unique. Many variations of this technique can be developed. For instance, if H,(f, 0) = 0, then we can use measurements of H,(f, 0,O) (assuming it is nonzero) or any other nonzero kernel “slice” to accomplish the same task.
Nonlinear
systems
identification
and modelling
995
It is also worth noting that the first order Wiener kernel H,(f) is equivalent to the “apparent transfer function” often evaluated in practice from input-output data and viewed as the best linear approximation of a given nonlinear system. It is evident from equation (25) that, for this class of systems, the apparent transfer function may depart significantly from the linear portion of the system [represented by K, (f)] depending on the value of the scalar II, and the function K:(f). From equation (23) we see that this scalar depends on the GWN input power level and all the odd coefficients of the polynomial terms. 3.
CONCLUSIONS
The class of nonlinear systems, described by equation (1) under the assumption of small magnitude coefficients IX,,has a Volterra functional expansion with kernels of the approximate form given by equation (15). This relatively simple form of the Volterra kernels leads to simple expressions for the Wiener kernels of these systems, given by equations (22) and (25). Since the Wiener kernels can be estimated in practice from input-output data, when the input is Gaussian white noise, an identification procedure is proposed that allows estimation of the unknown transfer function of the linear portion of the system and the coefficients of the polynomial nonlinearity. This method may find application to modelling and validation studies of systems of this class that can be tested with Gaussian white noise inputs. AcknoM,ledRemenrr-This
work was supported
by Grant
No. RR01861
from the National
Institutes
of Health
REFERENCES 1. E. Bedrosian and S. 0. Rice, The output properties of Volterra systems (nonlinear systems with memory) driven by harmonic and Gaussian inputs. Proc. IEEE 59, 1688-1708 (1971). 2. W. J. Rugh, Nonlinear System Theory: !he Volterra/Wiener Approach. Johns Hopkins Univ. Press, Baltimore, Md (1981). 3. M. Schetzen, The Volterra and Wiener Theories of Nonlinear Systems. Wiley, New York (1980). 4. V. Volterra, Theories of Functionals and of Inregral and Integro-differential Equations. Dover, New York (1930). 5. M. J. Korenberg, Identification of nonlinear differential systems, Proc. Joim Aufomatic Control Conf., pp. 597-603 (1973). 6. S. Baumgartner and W. J. Rugh, Complete identification of a class of nonlinear systems from steady-state frequency response. IEEE Trans. Circuits Syst. CAS-22, 753-759 (1975). 7. E. Wysocki and W. J. Rugh, Further results on the identification of problem for the class of nonlinear system S,,,. IEEE Trans. Circuits Syst. CAS-23, 664-670 (1976). 8. N. Wiener. Nonlinear Problems in Random Theory. Wiley, New York (1958). 9. P. Z. Marmarelis and V. Z. Marmarelis, Analysis of Physiological Systems: the Whire Noise Approach. Plenum Press, New York (1978). 10. Y. W. Lee and M. Schetzen, Measurements of the Wiener kernels of a nonlinear system by cross-correlation. Inr. J. Control 2, 237-254 (1965).