Copyright © 19% IFA C 13th Triennial World COlIgrc~s . San Francisco. USA
3a-063
Identifying MIMO Hammerstein systems in the context of Subspace Model Identification Methods Michel Verhaegen Dclft University of Technology, Department of Electrical Engineering, 2628 CD Delft, The Netherlands
[email protected] and David Westwick McGill University, Dept. of Biomedical Engineering, Montreal, Canada
Abstract
such as in [5 J, as well ment, such as in [7].
as in the repetition of che experi-
In this paper, we outline the extension of the MOESP (standing for Multivariable Output Error State sPace model identifi- (b) one-step non-iterative solutions: The key part here is to expand the static non -linearity into a series, which is gencation and introduced in [I].) family of subspace model idenerally a polynomial one. In this way the non-linear systification schemes to Hammerstein type of non-linear systems. is transfonned to a linear time-invariant system , tems One type of identification problem is considened. This type adenabling the application of standard linear identification dresses the identification of both the linear dynamic part and techniques. An example of this approach is [R] . the static nonlinearity, where only limited a priori information regarding the struct ure o f the nonlinearity is available. Another (c) rwo·step lion-iterative solutions: The first step is the iden· type of Hammerstein identification problem. considered in l2J. ti fication of the (linear) dynamic part and this infonnation assumes the (polynomial) structure of the static non-linearity to is then used in the second step to model che non -linearity. be given and the task here is to identify similarly the linear sysAn example of this approach is [9]. tem dynamics and the unknown proportio nal constants in the parametrization of the static non-linearity. The improved ro- In this paper, we focus on the class of non-iterative solutions bustness properlies of the algorithms developed in this paper (c). For the other solution helonging to class (b) we refer tu [2], The contributions in these classes respectively are: over existing corre lation based schemes is illustrated in Pl Keywords. Subspace Model Identification , Hammerstein (for class b) Generalizacion of the linear time-invari ant subNonlinear systems. space model identificatio n algorithms devel oped in the MOESP family [I. 10, 11. J 2, 13] for Hammerstein sys· terns. These algorithms in a natural way address the 1 Introduction. MIMO case and consider more general types of instrumental variables then the solution proposed in [8]. A problem of praccical relevance in Ihe area of identifying no n-linear systems ili the identification of so-called Hammer- (for dass c) Within this context, the first step is solved by a particular representative of the MOESP family, namely the stein systems. Such systems consist of a series connection of PI I scheme, instead of the widely used correlation apa memory-less (static) no n-linearity followed by a linear timeproach. furthermore, attention is given to incorporate invariant dynamical system. Examples of physical phenomena partial information on the structure of the non-linearity, that can accurately be described by this class of systems occur such as the presence of a dominant quadratic term. in solvin the area of non-lin ear wave force model s [3], l4J . ing this first step. Hammerstein models seems 10 be first propo.rred by Narendra and Gallman [5] . Since then a 10[ of contributions have been In this paper we assume (he stochastic sequences to be ergodic made. see (6 J for an overv iew. A crude way to subdi vide these and as such their analys is is performed within the framework contributions is in the following three classes : outlined in [I J. The outline of the paper is as follows. The Hammerstein identification problem is stated in section 2. Its (a) iterative solutions: These schemes are iterative in the way of computing an estimate of the ~tatic and dynamic part. 1 The acronym PI stallds for Pasllnpul imtrumental vari:lbl e~ Ill ].
4092
solution is treated in section 3. Some of the improvements of che schemes proposed in section 3 over (he existing correlation based approaches are highlighted in [2] . Finally section 4 contains some concluding remarks_ For the proofs of the Lemmas and Theorems stated in Ihis paper, we refer to [2]_
2
Hammerstein type of non-linear systems in state space form.
Example J Lelthefollowing Hanrmerslein system be given : Xk+ 1
AT-I:
+ 8uJ.- + uBui
zk
C:J:k
+ OUk + (l'Du~ + Vk
where Uk and Vk are independent scalar zero-mean while noise sequences with a Gaussian distribution. Now we attempt tu identify the quadruple {A. , B , C. DJ usillg the ilo sequences {u., z. /n Ihis case, let Ihe RQ faclOriZlllion step (cfr. Eq. (25) of {/ JI of the ordinary MQESP scheme be denoted as,
J.
In this paper, we assume the MIMO Hammerstein system to be represented ao;; follows:
A.F-, ("x,
+ BI" + DI" + v,
(1)
(ukl
(2) with Uj,~ , r.." and Zj .&,N block-Hankel matrices constructed from {u.} and {zd respectively. On the other hand, the data equation relating the Hankel mafriC.-·es Zj , ~ , N and Ui ,$,N now become~':
l
wherex-" E fR l) ./lk , ·U·k E IRm ,Yt,Vk E fR and<1>( _): fR m _ fA. III is a non -linear vector function . Unless stated more precisely, we assume Vir to be a zero-mean stationary stochastic process of arbitrary calor. Ba...ed on this system description, a formulation of the two identification probl ems analyzed in this paper is given next. Formulation: Let the following data sequences [Uj tI-j + I ... Uj +N-l ] and [ :j 2'j+l Zi+N - l ] of input and output (i/o) data of the system (I) be given. assume that the input sequence {Uk} is (strongly) persistently exciting, as defined in Definition 1(2) and statistically independent from the perturbation Vk, the task is to consistently estimate the quadrople (A, B , C, D) (up to a similarity transformatio n and an input transformation) without ao;suming detailed information on the nature of the non-linear map ~(.) _ The 'approximation' ofthis map is addressed in a succe.
3
The identification problem when restrictive information is available on the nature of the non-linear map <1>( ,).
Although the assumplion on the fixed polynomial structure of the static nOIl-linear map cIl( .) is quite standard in the literature. it is generall y not very realistic. Furthermore the direct approach discussed in r21 requires all the <.;omponents of Wk to be included. This may certainly for the MIMO case lead to a si gnificant increase in computational complexity, [n reality, instead of the detailed information a ... sumed in the type I of Hammerstein identificatio n problem considered in (21. one only has partial information on the nature of the non-linear map
7,j, .• ,N = r ,Xj ,N
+ H ,Uj,•. N + "H.U/"N + V; ",N(3)
with U/~ , N the block-Hallkel nultrix constructed/rom the sequellce {un and the matrices I'.. and H, defined as:
[
~~ 1
TJ CB
H.
o D
=[
CA", - l
CB
Using thejacturizarion (2), we can write Eq_ (3) as,
R21Ql
+ R'2'2Q'2 = r $X j, j\ ' + H8R.l1Ql + CkH~U/8,N
+V; ,8,i\' Mulliply this equation olllhe right by
R':l1
Qr yields,
= r"xj ,i\'Qf + H,H. l1 +aH~ Uj~ , N QT + V} .3 ,NQi
alld therefore,
RnQ, = r,XjN (r - QfQl)+aH.• U/, .N (r - QrQd(4) Hj ,. N (! - Q;Ql ) The assumption orl the Gaussian distribution of the input UJ; leads 10 Ihe following result /14/:
. I?,') -T J~oo N{jj~;.<,N Uj,$, N
=0
and when assuming that the matrix HIt is non-singular(guaranteed by the white lIoise property). this can be d enoted as,
.)wu/. ,NQf = O"' « ) where ON (() is a matrix of norm f. In the same way, the independenl'Y between Uk and vi /eads to,
_I_V .IN },."N QT1 = 0 N ('l·
4093
Sllb:,titulinR Ihese two expressions into Eq. (4) yields, . -r . ,2
tenns. In the next subsection, we show that the PI scheme, [111 , generates the same result under R2'JQ '2:::: r " X j', NQz Q '1 + UHtlr'j, ~ . N + Vj , ~ , N + O N(f ) these constraints. Fur an illustration of the improvements of The presence of th e n(m-varzishing term (for N 00) the PI scheme over the classical finite impulse response (FIR) o'H 6 U/$,N prevelUs the successful operation of lhe ordinary estimation strategy we refer to [21. MOESP scheme. Even when the perturbation Vk is zero-mean while lIoise, the column space of R'12 never becomes a con- 3.2 The non-linear map ~ (.) dominantly consislenl estimate of that of r,. As If consequence the estimated tains odd terms. quadruple n/system malri!-es will be biased. The first step of the PI scheme [11 J is the following RQ factorRemark 2 In the same way it can be shown that the PO scheme ization : fails to operate succes,lfully when use is made o/partial inforQ'~ ] (7) malion on the structure of the non-linear map ~( . ).
as derived by Verhaegen
Rt:, ] [
3.1
The data equation of the Hammerstein identification problem.
[n this and subsequent sections , we assume without loss of generality that the non-linear map ~(.) is parametrized as, N.,
(u.)
=L
(5)
Wq""
q:::; l
' JP
.T. q E IR mxm . ·, [ 'Ul,t ' h were Us: Urn ." and "' This pararnetrization. which does not allow for any cross-terms in the nonlineari ty. will simplify the upcoming theorems and derivations. In sect io n 3.4 we will consider the identification of systems with a more general nonlinearity. For this particular case, the data equation can now be denoted as.
The next theorem shows that fo r under the present constraint on Ihe non-linear map <1>( .). the twosubmatrices I t:;' and Rf:, play the same key role in identifying the LT! part of Hammerstein systems as in the original scheme to identify LTI syscems only [11]. As a consequence the calculation of the system matrices would proceed in exaccly the same way as outlined in [11]. Therefo re in the following we focus the discussion on this uniform role of the submatrices appearing in the third block row and first and second block column of the proposed RQ factorizations. Before presenting Ihis theorem. we need the following lemma.
Lemma 3 Let the input 'Uk be a stochastic random sequence with gaussiarz distribution a,ui id the RQ factorization in (7) be give,!. If -j;;R'0. and -;};:;n{,; are non-singular for large N , then,
~~~ ~U'+I .•. N(Q~ f =
N~
Zj , ~ , N
:::: 1',X j .•'Io·
+ H s L \[I q .tl U/~ , N + Vj ,I ,N
(6)
Dq limN~ ~
q==l
o
{ where U/~.N is the bloc k-Hankel matrix constructed from the q ~equence { tlir } and the matrix \If >J.~ equals:
IJ {I
1
Wq In the next two subsectio ns, we investigate the case where for some q the matrices "" are such that the pair (A liLW,) rt::mains controllable. Two ca..;;es will be considered, one where the this assumption hol ds for udd q (see section 3.2) and one where the assumption holds for even q (see section 3.3). When the input sequence is Gaussian white noise, it is well known. see e.g . l15J , that correlation (or spet..'tral analysi s) methods a1low to calculate an asymptotic unbiased estimate of the impulse response of the LTI (stable) system up to some proportional constant. Here one generally imposes the constraint that the non-linear map $(.) contains odd non-linear
g~
f,; R~l
q odd If e'ue n
~!~~ .Jr;,U'+l., .N (Q~l = 0 Theorem 4 Lel the input Uk la th e Hammerste;n system ( I) with the nOli-linear map 4l-(.) given as in Eq. (5) be zero·meal1 Gaussian lIoise, which is persislently exciting. Furthermore.
let u. be independent of Vi ('if k. () and let the RQfacto rization of Eq. (7) be given. then: .
lun
· Illn
N - oo
L N
N~' ~ hm
",R3 l
N - oo V
N __ =
1
NT
",f., X.,+l,N(Ql)
v :'V
+
(8)
I H - uN JR: ~ "=' 6 HII v .N
where ::: .. i.l· block diagonal with equal blocks on the diagonal, and .
hm
N_ =
4094
IN
",fl" v 11{
I N T = N-' co v""r,X Hl .N (Q2) .V 11111
(9)
When p (X,+I,N(Q:i?) = n, Eq, (9) shows that the key numerical step of the MOBSP family of algorithms, yielding a consistent estimate of the column space of r $. is now performed by a SVD of the matrix R~;. This is all in line with the conventional case of identifying the LTI system with algorithms belonging to the MOESP family ([1],[10],[ 11]), Theorem I does not guarantee identifiability, as it is entirely possible for 2$ to be singular, even though one or more of the qr q, (q odd), are themselves non-singular. To illustrate this point, consider the following example.
where the matrices [f~·~~,N correspond to the matrices 5 ,N for j = 1 and s + 1 but with zero-mean entries. For this factori7.ation, we have the following theorem. Since the proof of this theorem is completely analogue to that of Theorem 2, we state it without proof.
U/
Theorem 6 Let the input tl.k to the Hammerstein system (J) with the non-linear map >£1>(.) given as in Eq. (5) be persistently exciting zero-mean Gaussian noise. Let Uk be independent of Vi (Vk, t) and let the RQ factorization in Eq, (/0) be given, then:
Example 5 Let a SISO Hammerstein system (J) with a nonlinear map be described by Jlk
, hm
N ...... oo
= (tld == CIUk + CilHr
and let the variance afthe input (scaiar) be O'~. Then, as shown in [2 J we can express 3~ as, 2!
2('1 + 2'2 Thus. if O"~
, Ilm
will be zero, and the estimation
of the matrices Band D will fail, Clearly, this particular example will only be probLematical have opposite signs.
if the coefficients,
=
,1 -N T hm rvr,X'+1,N(Q,) + N-+(Xl
v.V
(11)
I -N ""H,'5."R l1 yN
where 2" is hlock diagonal with equal blocks on the diagonal, and,
4! 'I x 2! e3 u ;
= -:fZ;. then 25
lim
N ...... =
1 -N y N
""R 31
Cl
and
C3
N-+CQ
V
1 -RN = 32 = I'1111 N
N .... oo
1 r X "", ,+I,N (-QN)T 2 y.V
(12)
Remal'k 7 For the case, the non-linear map q,(.) contains both even and odd terms (in a dominant way). the above two variants of the PI scheme may be combined to yield a single estimate of the LT! system. One direct way to merge the two schemes is to use the original PI scheme with the extended in-
This particular difficulty also arises when correlation-based methods are used. The same input variances which result in 0 also cause the input-output cross-correlation to vanish. put u Tk ( Uk~")T JT . 25 From the preceeding example, we have seen how the identification of a SlSO Hammerstein system can fail. even though Remark 8 In the 'classical' case of estimating an FIR by the conditions of Theorem I are satisfied, and the nonlinearity means of least ...quares methods, Lemma 1 can be used to show contains only odd terms. Furthermore, repition of the experi- that the FIR can be identified asymptotically unbiasedly when ment, using a ditlerent input variance, can be used to restore the underlying system has a FIR (or is asymptotically stable), As such there i no diffirence fmm an asymptotica accuracy identifiability. This example suggests that a MIMO Hammerstein system point of view benveen this classical approach and variants of may not be identifiable, even though the input is persistently the PI scheme presented in this section. However; we refer the exciting. and the nonlinearity contains odd tenns. Identifica- interested reader to 12J for an illustration of the difference betion will fail if the diagonal blocks of the matrix 25 become tween the two procedures whenfillite data length sequences are singular. This, in turn, will depend on both the nonlinear map used. of the unknown system and the t:ovariance matrix of the test signal. A more detailed study of this problem, as well as the development of procedures to reslore identifiability, is a topic 3.4 General MIMO static nonlinearities. of ongoing research. So far we have restricted ourselves to nonlinear maps of the form given in (5). In this subsection, we will show that the 3.3 The non-linear map <1'>(,) dominantly con- PI method ([ 11]) can be used when the nonlinear map includes cross-terms. For simplicity of notation. we will deal with nontains even terms. linear maps that include a single cross-term, under the addiIn this case, we propose to use the following RQ factorization tional assumption that the input are mutually independent. but as the first step in PI scheme: not necessarily white. By superposition. this result can then be applied to maps with an arbitrary number of cross-terms. finally, for the sake of notational simplicity, we use a two-input example to outline the extension of these last results to the case where the input components are correlated.
=
l
4095
Thus. let the nonlinearity contain a single cross-tenn of the form
Then a (non-minimal) descriptor representat ion of the inverse of this system is given straightforwardly as:
[~ ~ ][ ~::,' 1 We will now consider the identification of a Hammerstein system, whose nonlinear map consists of a single cross-term, using the PI method ([11]). By superposition, thi s will describe the effect that the presence of such a cross-term will have on the identification of a Hammerstein system whose nonlinearity includes both cross-term and the additive parallel structure describes in (5). Considerthe expected value, El ","('11 , Q2, T)UJl. If q1 + q, is even, then this will be zero, and the tenn will have no effect on Ihe ideolification. Thus, without loss of generality. assume 11 to be evcn. and q2 to be
ur
Lemma 9 Let the input a non/in ear map be
tt k 10
I .un
I / ."
r;;; ) 5+1 s N
· D er I101
N- oo
"
(causal part) ( 14)
( ql,
q2) T
)(QN),!, _ 1 -
1 RN r;;r i 1 v ." r
where Dc I' i.f a block-diagonal matrix with equal blocks on the diagollal. Then given this lemma. we may prove Theorem I provided the input has independent components.
3.5
=
This representation is valid for the case D is si ngular. A more general approach which is valid for the case D is singular as well as when the LTl system is non-minimum phase is based on the variant of the MO ESP fami Iy of algorithms described in [17]. For that purpose, we first generate a data record {'Id fora selected record {~d using the model ( 13) with zero initial condition s. Then, second we identify a mixed causal, anti-causal system of minimal order using Ihe record {1J.t} as input and {{.} lIS output. The identified system has the following structured tonn:
a Hammerstein .system ( } ) with
where I} I is even., q'l is odd, and rtO, have independent, possibly nOIl-white, components, and let the RQjactorization (7) be given. Then, N - ~ vtv
~k
=
Approximating the non-linear map q>( uk l,
where .I(A) :S 1 and .I(E) < 1 for .I(.) denoting the spectral radius ofthe matrix (.). Having this mixed causal, anti-causal system description, we can use the output Zk to calculate an approximation of the input to the LT! system. This quantity, denoted by ijk. is affected by the (zero-mean) perturbati on Vi:. However, caused by the independency between lit and Ul' and the zero-mean property of 1Ir.. the effect of Vi; a'iymptotically cancels out in the process of modeling the non -linear map
4
Concluding remarks.
Once the LTl part of the Hammerstein system has been determined (up to some input scaling), the approximation of the In this paper we have revised the identification of Hammerstein non-linear map ~(.) may proceed via the calculation of the in- non-linear systems in the context of subspace model identifiverse of this LTJ system. see [161 for a more elaborate discus- cation. Apart from the straightforward application to MIMO sion on this topic. In this section , we embed this procedure in identification problems, the improvements of the presented a state space framework and consider MIMO systems. schemes over existing correlation based schemes have been Let the obtained system matrices {.4T,HT,CT,Jh} , hi ghlighted by a numerical experiment. The improvements eswhere BT = T B1'~- l , DT :::::: H'I~- J for non·singular square pecially hold for the practically interesting case of fini te data matrit:es T and 1;" represent (he foltowing dynamical system; length sequences. In the identification of the type of Hammerstein systems conr.k+ I :1-,.;):.1,: + R'I'f.k sidered in this paper, that is when only partial infonnation is ( 13) available of the non-linear map (~ ( . ) we have treated the ca')e W + UT~"
= = er"'.
4096
of odd and even terms in this map independently. In practice it might also be important to detect the dominance of either or both terms. This topic has to be considered for the schemes developed in this paper as well as for the classical correlation based methods. This broader analysis is a topic of further research. For an analysis of the case when precise information is available on the structure of the static non linearity ~(.) as well as for an illustration of the obtained results by simulation examples we refer the interested reader to [2].
Acknowledgement The research of Dr. Michel Verhaegen has been made possible by the funding of a senior research fellowship from the Royal Dutch Academy of Arts and Sciences.
References [I] M. Vcrhacgen and P. Dewilde, "Subspace Model Identification. Part!: The Output-Error State Space model identification class of algorithms," Int. J. Control, vol. 56, no. 5, pp. 1187-1210,1992. [2J M. Verh.egen and D. Westwick. 'identifying MIMO Hammerstein systems in the context of a subspace Mode1 Identification Strategy." to appear in Int. J. Control, 1996. l3J J.S. Bendat, Nonlillear System Analysis and Identificationfrom Random Data. Wiley-Interscience Publication, 1990. [4] D. Choi, R.WMiksad, and E.J. Powers, "Application of digital cross-bispectral analysis techniques to model the nonlinear response of moored vessel system in random seas," 1. (~f Sounds and Vibration. vol. 99, no. 3, p. 309, 1985. l5 J K.S. Narendra and P.G. Gallman, "An iterative Method for the Identification of Non linear Systems Using a Hammerstein Model," I.E.E.£' Trans. on Autom. Control, vol. AC-II, pp. 546-550, 1966.
[9] M. Pawlak, "On the Series Expansion Approach to the Identification of Hammerstein systems," l.E.E.E. Trans. on A"tom. Control, vol. AC-36, pp. 763-767, 1991. llO] M. Verhaegen and P. Dewilde, "Subspace Model Identification. Part IT: Analysis of the elementary Output-Error State Space model identification algorithm," Int. J. Control, vol. 56, pp. 1211-1241,1992. [11] M. Verhaegen, "Subspace Model Identification. Part Ill: Analysis ofthe ordinary Output-Error State Space model identification algorithm," Int. J. Control, vo!. 58, no. 3, pp. 555-586, 1993. l12] M. Verhaegen, "Application of a Subspace Model Identification Technique to identify LTI systems operating in closed loop," Automatica, vol. 29, no. 4, pp. 1027-1040, 1993. [13] M. Verhaegen, "Identification of the deterministic part
of MIMO state space models given in innovation form from input-output data," Special Issue on Statistical Signal Processing and Control ofAutomatica, vol. 30, no. I. pp. 61-74, 1993. ll4J J.S. Bendat and A.G. Piersol, Random Data: Analysis and Measurement Procedures. Wilcy-Interscience, 2nd ed., 1986. [15] W Greblicki and M. Pawl,lk, "Identificationof Discrete Hammerstein Systems using Kernel Regression Estimates," I.E.E.E. TraIlS. on Autom. Control, vol. AC-31, no. 1, pp. 74-77, 198fi. [16] M.J. Korenberg and I.W Hunter, "The identification of Nonlinear Biological Systems: LNL Cascade Models," Biological Cybernetics, vol. 55, pp. 125-134,1986. [17] M. Verhaegen, "A Subspace Model Identification solution to the identification of mixed causal, anti-causal LTl systems," to appear in SIAM J. Matrix AnalysiS and Applications, 1996.
[6) R Haber and H. Unbehauen. "Structure Identification of Nonlinear Dynamic Systems - A Survey on Input/Output Approaches," Autonwtica. voL 26, no. 4, pp. 651-677, 1990. [7] L. Zi-Qiang, "Controller Design Oriented Model Identification Method for Hammerstein Systems," Automatica, vol. 29, no. 3, pp. 767-771,1993. [8] P. Stoica and T. Sbderstrt)m, "Instrumental-variable methods for identification of Hammerstein systems," Int. J. Conllvl, vol. 35, no. 3, pp. 459-476, 1982.
4097