- Email: [email protected]
Identification for Wiener Systems with RTF Subsystems*
European Journal of Control (2006)6:581–594 # 2006 EUCA
Identification for Wiener Systems with RTF Subsystems Xiao-Li Hu and Han-Fu Chen** Institute of Systems Science, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, P. R. China
Recursive algorithms are proposed for identification of the Wiener system composed of a subsystem with rational transfer function (RTF) cascaded with a static nonlinearity. The nonlinear block is not necessary to be invertible. The strong consistency is proved for estimates for the coefficients of RTF and for the value of the nonlinear function at any point. Some numerical examples are given showing that the simulation results are consistent with theoretical analysis. Keywords: Recursive identification; Wiener systems; rational transfer function; strong consistency; stochastic approximation
1.
Introduction
The Wiener system consisting of a linear subsystem followed by a static nonlinearity is an important model in practice. For recent years, its identification issue has been attracting a great attention from researchers, e.g., [2,4,7,8,10–20] among others. The SISO Wiener system considered in the paper is presented by Fig. 1. To present the nonlinearity there are the parametric [1,5,9,12,14–20] and the nonparametric approaches [2,7,8,10,11,13]. The present paper belongs to the latter and is the continuation of [10,11]. If vk were observed, then identification of the linear subsystem would be an easy task. Therefore, it is important to extract information concerning vk from observations. In the case where fðÞ is invertible, vk can Supported by the NSFC under Grants No. G0221301, 60334040, 60474004. Correspondence to: H.-F. Chen, e-mail: [email protected]
be expressed by system output: vk ¼ f1 ðyk Þ, and such an expression may simplify analysis in many cases, in particular, for the case where the linear subsystem is also nonparametric [7,8,13]. However, the invertibility condition on fðÞ excludes many important nonlinearities including saturations from consideration. This motivated us to consider identification of Wiener systems without requiring invertibility of fðÞ in [10,11], where the recursive identification algorithms and their strong consistency were established for Wiener systems with general nonlinearities. However, the linear subsystem considered there is limited to an finite impulse response (FIR) system. In this paper the FIR system is extended to the RTF system, which clearly is with infinite impulse response. Similar to [10,11], the parameters of the linear part are first estimated, then the intermediate signals vk are approximated by using the estimated linear subsystem. Finally, the nonlinear block is identified on the basis of the approximated intermediate signals and the observations fzk g. Such an extension is technically nontrivial by the following observation. In the convergence analysis itP is of importance to have convergence of series like 1 k¼1 ak vk where fak g are the stepsizes used in the algorithms to be defined later on. P For an FIR linear subsystem the series 1 k¼1 ak vk can be split into a finite number of convergent subseries (due to independence), while for an RTF subsystem such a direct separation is impossible and the treatment is more delicate. The estimates for coefficients of the linear subsystem and for the value fðvÞ of the static nonlinearity Received 10 July 2006; Accepted 15 January 2007 Recommended by A. Karimi and L. Ljung
582
X.-L. Hu and H.-F. Chen
Fig. 1. Wiener system.
at any v are given by the stochastic approximation (SA) algorithms with expanding truncations [3]. These algorithms together with conditions for their convergence are presented in Section 2. The limiting behaviors of some correlation functions and series are demonstrated in Section 3. In order to deal with the difficulty arising from the infinite impulse response of the linear subsystem mentioned above, the totality of integers is separated into a countable union of disjoint subsets, and convergence of several series is established with the help of such a separation. This is done in Section 4. The strong consistency of the estimates given in Section 2 is proved in Section 5, and a few numerical examples are given in Section 6. In the last section some concluding remarks are given.
2. 2.1.
Problem Setting and Estimation Algorithms Problem Setting
Let the linear subsystem be with RTF CðzÞvkþ1 ¼ DðzÞuk ,
ð1Þ
where CðzÞ ¼ 1 þ c1 z þ . . . þ cp zp ,
ð2Þ
DðzÞ ¼ 1 þ d1 z þ . . . þ dq zq
ð3Þ
are polynomials in backward-shift operator zuk ¼ uk1 . It is noticed that in the special case CðzÞ ¼ 1, the linear subsystem is with FIR and is considered in [10,11]. The purpose of this paper is to estimate the coefficients ci , i ¼ 1, . . . , p and dj , j ¼ 1, . . . , q and the value fðvÞ of the nonlinearity at any v on the basis of the designed inputs fuk g and observations fzk g: zk ¼ yk þ " k ,
yk ¼ fðvk Þ,
where f"k g is the observation noise. We now list conditions to be used in the sequel.
ð4Þ
H0 CðzÞ and DðzÞ are coprime and CðzÞ is stable. H1 The function fðÞ satisfies the following condition jfðxÞ fðyÞj < Lð1 þ jxj þ jyj Þjx yj 8x,y 2 R; ð5Þ where L > 0 and 0 are some constants. H2 fk g is a sequence of mutually independent random variables with Ek ¼ 0 and supk E2k < 1, and is independent of fuk g. Remark 1. Condition H1 proposed here is probably used for the first time. H1 is more general than the global Lipschitz condition, which requires ¼ 0 in (5), but it is stronger than the local Lipschitz condition, by which in lieu of Lð1 þ jxj þ jyj Þ in (5) there may be any continuous function of x, y. Under H1 the growth rate of fðzÞ is not faster than a polynomial as z tends to infinity. This property is used in (40) to guarantee an acceptable approximation to yk whenever vk is well estimated, which is needed for deriving convergence of some stochastic series in Lemma 8. Remark 2. We do not claim any necessity of H1. As shown later (Section 6, Example 3), even if H1 does not hold, the identification algorithm given in the paper still provides satisfactory estimates for both linear and nonlinear subsystems. However, it is not clear if it works in a feedback case. Since fuk g is at designer’s disposal, let uk ¼ 0, 8k 0. Then under Condition H0, (1) can be expressed as k1 X hj ukj , h0 ¼ 1, ð6Þ vkþ1 ¼ C1 ðzÞDðzÞuk ¼ j¼0
where the coefficients hj exponentially decay: jhj j ¼ O ej , > 0: 2.2.
ð7Þ
Stochastic Approximation Method
Since our identification algorithms are based on an improved version of classical stochastic approximation (SA) algorithm, we first briefly introduce SA and its improved version.
583
Identification for Wiener Systems
For an unknown function gðÞ with single root denoted by x0 , the classical Robbins-Monro (RM) algorithm, being the first version of SA, suggests to seek the root x0 by the simple recursion: xkþ1 ¼ xk þ ak gkþ1 ,
k ¼ 0, 1, . . . ,
ð8Þ
where
gkþ1 ¼ gðxk Þ þ kþ1
ð9Þ
is the noisy observation of gðxk Þ, xk is the k-th estimate for x0 , kþ1 is the observation noise, and fak g is a sequence of positive real numbers tending to zero as k ! 1. For convergence of xk to x0 , a certain growth rate restriction on gðÞ or an a priori boundedness assumption on fxk g are required in the classical theory in addition to conditions on the noise fk g. These requirements quite restrict the application of SA. Moreover, the noise is usually required to analyze for the whole sequence of estimates fxk g. In order to extend the range of applicability of SA method, the RM algorithm (8) is modified to be truncated at expanding bounds in [3]. By this the increasing rate of fxk g is controlled. Let fMk g be an increasing sequence of positive real numbers with Mk ! 1, and let xk be generated by k!1
the following algorithm: xk þ ak gkþ1 , if jjxk þ ak gkþ1 jj Mk , xkþ1 ¼ otherwise, x ; ð10Þ k ¼ k1 þ I½jxk1 þak1 gk j>Mk1 ,
0 ¼ 0,
ð11Þ
where function IA denotes the indicator function of a set A, and k is the number of truncations occurred until time k 1. From (11) it is seen that k k 1, and from (10) jjxkþ1 jj jjx jj ^ Mk1 . This means that the growth of jjxkþ1 jj should not be faster than that of Mk . At any time if jjxk þ ak gkþ1 jj exceeds the truncation bound Mk , then we pull xkþ1 back to the fixed point x and simultaneously extend the truncation bound from Mk to Mk þ1 . Otherwise, it develops as (8). Since the region where x0 is located is unknown, it is important to allow xk to grow up in order to have possibility to reach x0 . For convergence of fxk g we must have ak gðxk Þ ! 0. Thus, it is reasonable to require k!1
ak ! 0. However, if ak tends to zero too fast so that Pk!1 1 may be unable to k¼1 ak < 1, then the estimate xkP reach x0 . Thus, we have to require 1 k¼1 ak ¼ 1 (see 1 [3]). A common choice is ak ¼ k. The selection of
fMk g is defined by that how fast we allow fxk g to grow up. The following fact to be used in what follows is a special case of the general convergence theorem (GCT) Theorem 2.2.1 of [3]. Let gðxÞ ¼ ðx x0 Þ and ak ¼ 1=k. Then xk defined by (10)(11) converges to x0 for any sample paths satisfying mðnX k , Tk Þ ai iþ1 ¼ 0 8Tk 2 ½0, T lim lim sup T1 i¼n T!0 k!1 k
ð12Þ for any nk such that fxnk g converges, where m X aj T : mðk, TÞ ¼ max m : j¼k
Further, if fk g in (9) for the sample path ! under consideration can be decomposed into two parts k ¼ k0 þ k00 such that 1 X
0 a1 k k < 1,
k00 ¼ Oðak Þ
k¼1
for some 2 ð0, 1, then by the convergence rate theorem (CRT) Theorem 3.1.1 of [3] xk given by (9)(10)(11) converges to x0 with the following rate: k xk x0 k ¼ oðak Þ: It is worth noting that Condition (12) is required to verify only along convergent subsequences fxnk g rather than along the whole sequence fxn g. As pointed in [3], in many cases (12) cannot be directly verified along the whole sequence fxn g, but it can be done for convergent subsequences. This is the case in Theorem 2 of Section 5. The method verifying Condition (12) along convergent subsequences has shown a great advantage over verification along the whole sequence, and it is called as trajectory-subsequence (TS) method in [3]. 2.3.
Identification Algorithms
As stated in introduction, we first estimate the parameters ci , i ¼ 1, . . . , p, dj , j ¼ 1, . . . , q of the linear part, then construct v^k to approximate the output vk of the linear subsystem; finally, identify the nonlinear block by using the data f^ vk , zk g. We take the input sequence fuk g1 1 to be iid and uk 2 N ð0, 1Þ. Under condition H1, let us set Z 1 X t2 1 2 2 hi , and ¼ pffiffiffiffiffiffi tfðtÞeð22 Þ dt: ¼ 2 3 R i¼0 ð13Þ
584
X.-L. Hu and H.-F. Chen
P 2 Clearly, vk obeying N ð0, k2 i¼0 hi Þ is asymptoti2 cally distributed as N ð0, Þ, and is the asymptotic correlation coefficient of the input and output of the nonlinearity. Define the real sequences to be used in algorithms: ak ¼ k1,
bk ¼ k1 ,
Mk ¼ M0 þ k
and
ð14Þ
with > 0, > 0, M0 > 0, and 3 þ < 12. Similar to [4,10,11], hi , i ¼ 0, 1, . . . , p þ q are estimated by the SA algorithms with expanding truncations: kþ1 ðiÞ ¼ ½ k ðiÞ ak ð k ðiÞ uk zkþiþ1 Þ I½j k ðiÞak ð k ðiÞuk zkþiþ1 ÞjM k ðiÞ ,
ð15Þ
Then (19) can be written in the vector form ½hqþ1 hqþ2 . . . hqþp T ¼ ½c1 c2 . . . cp T ,
which is called the Yule–Walker equation. Replacing hi , i ¼ q, . . . , p þ q, in (21) with their estimates ðiÞ k if k ð0Þ 6¼ 0 hik ¼ k ð0Þ 0 otherwise and denoting by k the resulting matrix corresponding to , we derive a linear equation with respect to ½c1 . . . cp T at time k. The solution ½c1k . . . cpk T whenever exists to the obtained linear equation serves as an estimate of ½c1 . . . cp T , i.e., T ½c1k . . . cpk T ¼ 1 k ½hqþ1, k . . . hpþq, k :
k ðiÞ ¼ k1 ðiÞ þ I½j k1 ðiÞak1 ð k1 ðiÞuk1 zkþi Þj>M k1 ðiÞ , ð16Þ
0 ðiÞ ¼ 0
E½vkþlþ1 uk ¼ hl ,
l ¼ 0, 1, . . . :
ð17Þ
Multiplying both sides of (1) by uki and taking the mathematical expectation lead to hi ¼ E½vkþ1 uki ¼
i^p X
cl hil þ di ,
i ¼ 1, . . . , q,
l¼1
hj ¼ E½vkþ1 ukj ¼
i^p X
clk hil, k ,
i ¼ 1, . . . , q:
Thus, we have estimated all coefficients of the linear subsystem. It remains to estimate fðvÞ at any fixed v. Let us use these estimates to approximate the intermediate signal vk first. Define 0 1 1 0 0 c1 B c 0 1 0 C B C 2 B C B C C ¼ B C , B C 0 1 A @ cp1 0 cp 0 0 0 pp 0 1 1 d 1 dq B 0 0 0 C B C D¼B , C @ A 0
0
Uk ¼ ½uk . . . ukq T ,
0
pðqþ1Þ
H ¼ ½1 0 . . . 01p
xkþ1 ¼ Cxk þ DUkþ1 :
cl hjl ,
ð23Þ
l¼1
ð18Þ j^p X
ð22Þ
Motivated by (18), the estimates for di are given by dik ¼ hik þ
with arbitrary initial values 0 ðiÞ, i ¼ 0, 1, . . . , p þ q. We note that k ð0Þ serves as the estimate for , and ð k ðiÞ= k ð0ÞÞ if exists serves as the estimate for hi , i ¼ 1, . . . , p þ q. In order to estimate the coefficients of the linear subsystem let us first derive the corresponding Yule– Walker equation. Since fuk g is iid with Euk ¼ 0, from (6) it follows that
ð21Þ
ð24Þ
l¼1
j ¼ q þ 1, . . . , p þ q,
ð19Þ
where and hereafter a ^ b denotes minða, bÞ. Define 0
hq h B qþ1 ¼B @ hpþq1
hq1 hq hpþq2
1 hqpþ1 hqpþ2 C C, A hq
where hi ¼ 0 for i > 0 by definition.
It is clear that with x0 appropriately chosen the output vk of the linear subsystem coincides with Hxk : vk ¼ Hxk :
ð20Þ
ð25Þ
Replacing the elements of matrices C and D with their estimates derived above gives Ck and Dk which naturally serve as estimates for C and D, respectively. Define x^kþ1 ¼ Ckþ1 x^k þDkþ1 Ukþ1 ,
v^k ¼ H x^k
ð26Þ
with an arbitrary x^0 . We use v^k as an estimate for vk .
585
Identification for Wiener Systems
For estimating fðvÞ we need the kernel function
wk ¼
pffiffiffiffi 2 bk exp
2 vkbv þ k
v2k 220
ð27Þ
,
where 0 > 0 is a fixed constant. The kernel function wk can be estimated by
^k ¼ w
pffiffiffiffi 2 bk exp
(
v^k v 2 þ bk
) v^2k 220
:
ð28Þ
Finally, fðvÞ is estimated by k ðvÞ given by the following recursive algorithm:
Consequently, we have lim Ewk yk
k!1
Z pffiffiffiffi
x2 22 0
x2
2 ðxv bk Þ þ
1 22 k fðxÞpffiffiffiffiffiffi e dx 2 k (
) Z ðvþbk tÞ2 1 1 t2 1 dt ¼ lim e fðvþbk tÞk exp k!1 R 2 20 2k pffiffiffi ¼ fðvÞðvÞ: &
¼ lim b2 e k!1 k R
Lemma 2. Assume H1 holds. Then lim E½ukj ykþ1 ¼ hj ,
j ¼ 0, 1, . . . , p þ q,
k!1
^k ðk ðvÞ zk Þ kþ1 ðvÞ ¼ ½k ðvÞ ak w I½jk ðvÞak w^k ðk ðvÞzk ÞjMk ðvÞ ,
ð29Þ
k ðvÞ ¼ k1 ðvÞ þ I½jk1 ðvÞak1 w^k1 ðk1 ðvÞ yk1 Þj > Mk1 ðvÞ ,
0 ðvÞ ¼ 0 ð30Þ
with any given initial value 0 ðvÞ.
3.
Limiting Behaviors of Some Random Variables
We now establish the limiting behaviors of some random variables connected with the kernel and the a.s. convergence of some series. Lemma 1. Assume H1 holds. For wk defined by (27) the following limits take place: pffiffiffi pffiffiffi ðvÞ, lim E½wk yk ¼ fðvÞðvÞ, k!1 k!1 pffiffiffi ð31Þ lim E½wk jyk j ¼ jfðvÞjðvÞ, lim Ewk ¼
k!1
pffiffiffiffiffi pffiffi lim E½ bk wk 2 ¼ 2ðvÞ, k!1 pffiffiffiffiffi pffiffi lim E½ bk wk yk 2 ¼ 2fðvÞðvÞ,
where and are defined by (13). Proof. Noticing that ukj and vkþ1 are jointly Gaussian and by (6) h i h E ukj 2j vkþ1 vkþ1 ¼ 0, k
h
we find that ukj 2j vkþ1 is independent of vkþ1 . k Taking the conditional expectation with respect to vkþ1 leads to
hj E ukj jvkþ1 ¼ 2 vkþ1 : k By this we have
E ukj fðvkþ1 Þ ¼ E E ukj fðvkþ1 Þjvkþ1
¼ E fðvkþ1 ÞE ukj jvkþ1 ¼ hj
Lemma 3. Assume H2 holds. Then the following series converge a.s.:
n 2 v 2
1 0
12
ak w k k ,
k¼1 1 X
k!1
where ðvÞ ¼ 1 exp
E½vkþ1 fðvkþ1 Þ , 2k
which by (13) yields the desired result by letting k ! 1 because 2k ! 2 and vk is asymptotically & distributed as N ð0, 2 Þ.
1 X
ð32Þ
ð33Þ
1 X k¼1
ak wk ðjk j Ejk jÞ, ð34Þ
ak uk kþj :
k¼1
o .
Proof. It suffices to show limk!1 Ewk yk , and the rest can be provedP in a similar manner. k2 hj ukj1 , and hence vk 2 N ð0, 2k Þ By (6) vk ¼ j¼0 P k2 2 with 2k ¼ j¼0 hj .
Proof. The proof is essentially based on the following fact (see Chapter 11.1, Theorem 3 of [6]): Let fSn , F n , n 1g be a martingale with supn EjSn j < 1, and let fYn , F n1 , n 1g be with supn jYn j < 1 a.s. 1 P Xn Yn < 1 a.s. where Xn ¼ Sn Sn1 . Then n¼1
586
X.-L. Hu and H.-F. Chen
For example, ak wk k can be written as ak wk k ¼ Xk Yk 2 . Since and Yk ¼ exp vkbv k
pffiffiffiffi 2 with Xk ¼ k1
k
0 < < 1/6, applying the above cited result leads to 1 P Xn Yn < 1 a.s. n¼1
The proof for convergence of the second series is similar. The convergence of the last series follows from the fact that fuk kþj g is a martingale difference sequence. In what follows we need the divergence rates of uk and yk as k ! 1. Lemma 4. For any > 0 the following rate takes place uk a:s: ! 0: k k ! 1
ð35Þ
If C(z) is stable and H1 holds, then k a:s: ! 0, k k ! 1
yk a:s: ! 0: k k ! 1
vk a:s ! 0, k k!1
1 X juk j P > " < 1: k k¼1
4.
&
Convergence Analysis via Separation of Integer Set
Since the linear subsystem is with RTF, its output is no longer a combination of inputs of the system. This makes convergence analysis for some series difficult. To overcome the difficulty we truncate the expression (6) for vkþ1 as follows vkþ1 ¼
k X i¼0
hi uki ,
k ¼ ½k , 0 < < 1,
ð39Þ By (5) of Condition H1 and Lemma 4, we further have vkþ1 Þj jykþ1 ykþ1 j ¼ jfðvkþ1 Þ fð < Lð1 þ jvkþ1 j þ j vkþ1 j Þjvkþ1 vkþ1 j ¼ Lð1 þ jvkþ1 j þ j vkþ1 j ÞO e2k a:s:, ð40Þ ¼ O e 4 k and similarly,
which incorporating with H1 implies
yk a:s ! 0. k k!1
Let us explain the idea of truncations applied in (37). Assume ¼ 1=2. Then vk2 þ1 is a linear combination of fuk2 k , uk2 kþ1 , . . . , uk2 g for any positive integer k. By the independence of the sequence fuk g and the fact k2 k > ðk 1Þ2 , it is clear that the subsequence f vk2 þ1 , k ¼ 0, 1, . . .g now turns to be a sequence of mutually independent random variables. On the other hand, when the linear subsystem is stable, the slowly truncated vkþ1 approximates vkþ1 as k ! 1 based on the following analysis. By (7) and (35) it is clear that X k1 hi uki ¼ O e2k jvkþ1 vkþ1 j ¼ a:s:
ð36Þ
Then by the Borel–Cantelli lemma (35) follows. By stability of C(z) from (24)(25) and (35) we have
ð38Þ
i¼kþ1
Proof. For any > 0 by the Tchebychev inequality " # juk j juk j2= 1 P >" ¼P > "2= < 2= 2 Ejuk j2= , k2 k " k which implies
where and hereafter ½x denotes the integer part of a real number x. Correspondingly, we set (
) pffiffiffiffiffiffi vk v 2 v2k 2 k ¼ yk ¼ fð vk Þ, w exp þ 2 : bk bk 20
kþ1 j ¼ O e4k jwkþ1 w
a:s:,
ð41Þ
2kþ1 j ¼ O e4k jw2kþ1 w
a:s:
ð42Þ
PTo analyze convergence of the series like ak ðuk ykþ1 Euk ykþ1 Þ, we now separate the totality of integers into a countable union of disjoint sets. Let f k g be a sequence of strictly increasing integers such that 0 ¼ 0, 1 ¼ 1, f kþ1 k g also nondecreasingly diverges to infinity, and 0 < lim inf k!1
k ka
lim sup k!1
k ka
<1
for some a > 1: ð43Þ
Define ð37Þ
A0 ð0Þ ¼ f k : k ¼ 1, 2, . . .g;
ð44Þ
587
Identification for Wiener Systems
Aj ðiÞ ¼ f k þ i : k ¼ j, j þ 1, . . .g,
Aj ¼
i[ 2 ðjÞ
Aj ðiÞ,
ð45Þ
j ¼ 1, 2, . . . ,
l ðj, iÞ ¼
i¼i1 ðjÞ
where i runs from i1 ðjÞ ¼ j j1 to i2 ðjÞ ¼
jþ1 j 1. If i2 ðjÞ < i1 ðjÞ, then set Aj ¼ . We denote by J the totality of integers f1, 2, . . .g. Lemma 5. The integer sets Aj ðiÞ, j ¼ 0, 1, . . . , compose a disjoint separation of J, i.e., Aj ðiÞ, i ¼ i1 ðjÞ, . . . , i2 ðjÞ, j ¼ 0, 1, . . ., are disjoint and A0 ð0Þ [
1 i[ 2 ðjÞ [
Aj ðiÞ ¼ J:
Set
ð46Þ
jas1 1=2 s j1þl þi :
j1þl þ i
By (43) we have ! jas1 1=2 jas1 1=2 1 s < ¼O
sj1þl
j1þl þ i ðj 1 þ lÞ1 þ1=2
1 ð50Þ ¼ O þ1=2 , l1 which combining with (49) yields El ðj, iÞ ¼ 0 and
j¼1 i¼i1 ðjÞ
El2 ðj, iÞ < Proof. The assertions of the lemma follow from the following set operations by noticing that f k þ 1, . . . , kþ1 1g, k ¼ 1, 2, . . ., are disjoint: 1 i[ 2 ðjÞ [
Aj ðiÞ ¼
j 1 [ 1 jþ1[ 1 [
f k þig
j¼1 i¼ j j1 k¼j
j¼1 i¼i1 ðjÞ
¼
j 1 1 [ 1 jþ1[ [
f k þig
¼
C l1þ21
ð51Þ
for some constant C > 0. By the condition of the lemma for fixed j and ifl ðj, iÞg is a sequence of mutually independent random variables with possible exception of a finite number of terms. Then, by the convergence theorem for martingale difference sequences (mds)
j¼1 k¼j i¼ j j1 1 [ k [
Sji ðmÞ ¼
f k þ j j1 ,
k þ j j1 þ1, ..., k þ jþ1 j 1g 1 [ & ¼ f k þ1, k þ2, ..., kþ1 1g: Lemma 6. Let fk g be a sequence of random variables with zero mean and supk Ek2 < 1. If for any fixed i and j, i ¼ i1 ðjÞ, . . . , i2 ðjÞ whenever i2 ðjÞ i1 ðjÞ, j ¼ 1, 2, . . . , the subsequence fk : k 2 Aj ðiÞg is composed of mutually independent random variables with possible exception of a finite number of k , then
for s > 3 2
1 a,
l ðj, iÞ
a.s. converges to some random variable Sji as m ! 1. Given > 0 by (50)(51) we have
P jSji ðmÞj > j
k¼1
a:s:
m X l¼1
k¼1 j¼1
1 X 1 <1 s k k k¼1
ð49Þ
ð47Þ
where a is given in (43).
3 1 1 a1 Proof. From a1 s > 2 a 1it is clear that as a þ 2 > 2 . Let 2 2 , as a þ 2 and let 1 > 0 such that 1 ð48Þ as a þ 1 > 0: 2 Since s > 12, under the conditions of the lemma it is P 1 clear that ks k < 1 a.s. Consequently, the k2A0 ð0Þ
index set A0 ð0Þ can be excluded from consideration.
1 2 j2
EjSji ðmÞj2
1 C X 1 C0 ¼ 2 2 , 2 2 1þ2 1 j j l¼1 l
P 1 where C0 ¼ C 1 l¼1 l1þ21 . Letting m ! 1 we have
P jSji j > j C2 j02 :
ð52Þ
Define j ðÞ¼ ! : jSji j j , i ¼ i1 ðjÞ, . ..,i2 ðjÞ , ð53Þ 1 \ j ðÞ: ðÞ¼ j¼1
In the case a 2 it is clear that i2 ðjÞ i1 ðjÞ 1 ¼ jþ1 2 j þ j1 ¼ Oððj þ 1Þa 2ja þ ðj 1Þa Þ ¼ Oðja2 Þ:
588
X.-L. Hu and H.-F. Chen
Thus, similar to (55) for ! 2 j ðÞ by Lemma 5 it follows that
By (52)(53) it follows that !
1 [
PððÞÞ ¼ 1 P
cj ðÞ
>1
j¼1
>1
1 X j¼1
1 X P cj ðÞ j¼1
i2 ðjÞ X 1 1 X X 1 jSji j k s as k¼1 k j¼1, A 6¼ j 1 1=2 i¼i ðjÞ j
C0 O ja2 2 2 j
¼
! 1 1 X 1 >1O 2 : j¼1 j2þ2a
j¼1, Aj 6¼
ð54Þ
By the definition of , we have that 2 þ 2 a > 1 and, hence, PððÞÞ ! 1. !1
For ! 2 ðÞ by Lemma 5 it follows that i2 ðjÞ X X 1 1 X 1 X 1 1 ¼ k s jþl1 þi s k¼1 k j¼1 ð jþl1 þ iÞ i ðjÞ l¼1
¼ O ð 1Þ
j¼1 i¼i1 ðjÞ
¼O
!
jas1 ða2Þ1=2
i
< 1:
< 1:
!1
Therefore,
0
>1
jas1 1=2
j¼1
i¼1
Similar to (56), for j ðÞ defined by (53) we have ! : jSji j kj , , if i2 ðjÞ ¼ i1 ðjÞ ¼ i, j ðkÞ ¼ otherwise:
1 [
1 X
j¼1, Aj 6¼ 1 X
1 cj ðÞA
P cj ðÞ
C0 2 j2 j¼1, Aj 6¼ 1 X 1 ¼ 1 O ð 1Þ 2 ! 1: 2 ia1 !1 i¼1
>1
1
1 1 a1 ðas1 1=2Þ
j¼1, Aj 6¼
jSji j
1 X
1 X
PððÞÞ ¼ 1 P@
X i2 ðjÞ X 1 1 X l ðj, iÞ ¼ as 1=2 j¼1 i1 ðjÞ l¼1 j 1
Oð1Þ jas1 1=2
To complete the proof it remains to show PððÞÞ ! 1.
1
i2 ðjÞ 1 X X
1
1 X
&
ð55Þ
Lemma 7. Assume C(z) is stable. Let s 2 ð12, 32Þ, P2 ð0, s 12Þ, and vk be defined by (37). Then 1 ks vk < 1 a.s.
In the case 1 < a < 2 we note that for sufficiently large j, ð jþ1 j Þ ð j j1 Þ equals either 0 or 1. In other words, i2 ðjÞ i1 ðjÞ equals either 1 or 0, and correspondingly,
Proof. Since fuk g is iid and Euk ¼ 0, we have E vk ¼ 0 1 2 , 32s Þ. It is clear that and supk v2k < 1. Take a 2 ð1 3 1 s > 2 a, which is required in Lemma 6. It suffices to show that f vk : k 2 Aj ðiÞg is mutually independent. If this is done, then the assertion of the lemma follows from Lemma 6. By (37) v½ðkþ1Þa þi is a linear combination of u½ðkþ1Þa þi , . . . , u½ðkþ1Þa þi½ðkþ1Þa þi , where by (31)
Aj ¼
i[ 2 ðjÞ
Aj ðiÞ ¼
i¼i1 ðjÞ
Aj ðiÞ; if i2 ðjÞ ¼ i1 ðjÞ ¼ i, , otherwise: ð56Þ 1
By noticing that i1 ðjÞ ¼ i implies j ¼ Oðia1 Þ, for ! 2 j ðÞ we have jSji j as j 1 1=2
! 0:
jas1 1=2 j!1
½ðk þ 1Þa þ i ¼ ½ð½ðk þ 1Þa þ iÞ . The total number of these terms is ½ðk þ 1Þa þ i þ 1 ¼ Oðka Þ. Let us denote by Uk the minimal subset of fuj , j ¼ 1, 2, . . .g such that v k þi is the linear combination of elements from Uk. On the other hand, the neighboring terms in f vk , k 2 Aj ðiÞg are with subscripts kþ1 þ i and k þ i, and ð kþ1 þ iÞ ð k þ iÞ ¼ Oðka1 Þ. This means that
589
Identification for Wiener Systems
starting from some large enough k0 fUk , k ¼ k0 , k0 þ 1, . . .g are disjoint subsets of fuj , j ¼ 1, 2, . . .g whenever a < a 1. Therefore, f vk : k 2 Aj ðiÞg is a sequence of mutually independent random variables with possible exception of finite number of terms. Thus, the conditions of Lemma 6 are satisfied. Lemma 8. Assume C(z) is stable and H1 holds. Then the following series converge a.s. 1 X ak ðukj ykþ1 Eukj ykþ1 Þ 8j 0, k¼1 1 X
ð57Þ ak ðwk yk Ewk yk Þ,
k¼1 1 X
k ð0Þ ! k!1
a:s:,
k ðiÞ ! hi k!1
ð61Þ
a:s:,
i ¼ 1, . . . , p þ q with convergence rates j k ð0Þ j ¼ oðk Þ a:s:,
j k ðiÞ hi j
ð62Þ
a:s
ð58Þ ak ðjwk Ewk j Ejwk Ewk jÞ,
in addition, 6¼ 0 and det 6¼ 0, then cik ! ci , k!1
a:s
i ¼ 1, . . . , p, dik ! dj , j ¼ 1, . . . , q with the same k!1
rate as that given in (62), where cik and djk are given by (22)(23), respectively. ak ðwk jyk j E½wk jyk jÞ:
ð59Þ
k¼1
Proof. Let us prove convergence of the first series in (57). The rest can be done by a similar treatment. By (35)(40) and H1 it is clear that 1 X k¼1 1 X
ak ukj ðykþ1 ykþ1 Þ < 1 a:s:,
and
ak ½Eukj ðykþ1 ykþ1 Þ < 1 a:s:
k¼1
Thus, for the first series in (57) it suffices to show 1 X
Theorem 1. Assume H0-H2 hold. Then
where 0 < < 1=2 and is given in (13). Further, if,
ak ðwk Ewk Þ,
k¼1 1 X
We now in a position to establish the strong consistency of estimation algorithms given in Section 2.
¼ oðk Þ a:s:, i ¼ 1, . . . , p þ q,
k¼1 1 X
5. Convergence of Identification Algorithms
ak ðukj ykþ1 Eukj ykþ1 Þ < 1 a:s:
ð60Þ
k¼1
Notice that Eðukj ykþ1 Eukj ykþ1 Þ ¼ 0 and sup Ejukj ykþ1 Eukj ykþ1 j2 < 1: k
Similar to Lemma 7, separating fuj , j ¼ 1, 2, . . .g to the disjoint subsets Uk with s ¼ 1 and 2 ð0, 1=2Þ, we find that fukj ykþ1 Eukj ykþ1 , k 2 Aj ðiÞg is a sequence of mutually independent random variables with possible exception of a finite number of terms. Then (54) follows from Lemma 6. & Remark 3. As shown in Lemma 6, the stepsize ak in (57)-(59) may be replaced by k1s with s > 12.
Proof. Rewrite (15) as 8 k ðiÞak ð k ðiÞhi Þak kþ1 ðiÞ, > > > < if j ðiÞa ð ðiÞh Þa ðiÞj k k k i k kþ1 kþ1 ðiÞ¼ > M ðiÞ,
k > > : 0, otherwise, where kþ1 ðiÞ ¼ uk zkþiþ1 þ hi ,
ð63Þ
i ¼ 0, 1, . . . , p þ q:
ð64Þ In the sequel we will use the general convergence theorem (GCT) and the convergence rate theorem (CRT) for the SA algorithm with expanding truncations mentioned in Section 2.2 and cited in [11] with reference to [3] for detailed proof. From (63) it is seen that by GCT with regression functions gðiÞ ðxÞ ¼ ðx hi Þ for (61) it suffices to show mðn, tÞ 1 X lim lim sup aj jþ1 ðiÞ ¼ 0 T!0 n!1 T j¼n
ð65Þ
8t 2 ½0, T,
i ¼ 0, 1, . . . , r, ( ) m P where mðn, tÞ ¼ max m : aj t . j¼n Noticing kþ1 ðiÞ ¼ ðuk ykþiþ1 E½uk ykþiþ1 Þ ðE½uk ykþiþ1 di Þ uk kþiþ1 , we conclude (65) by Lemma 2 and the convergence of the first series in (57).
590
X.-L. Hu and H.-F. Chen
For (62) it suffices to set ¼ 1, F ¼ 1, and 00 k ¼ 0 in CRT stated in [11]. The assertions concerning ci and dj follow from (22), (23) and the assumptions 6¼ 0 and det 6¼ 0. &
By stability of C and by the convergence Ck ! C, there exist N1 > 0 and 2 ð0, 1Þ such that
Theorem 2. Assume that H0-H2 hold, 6¼ 0 and det() 6¼ 0 where and are given in (13) and (20), respectively. Then for k(v) given by (29)(30)
where
k ðvÞ ! fðvÞ
Sð , kÞ ¼ N2
ð66Þ
a:s:
k!1
x0 x0 jj þ Sð , kÞ, jj x^kþ1 xkþ1 jj N1 kþ1 j j^
kþ1 X
ð69Þ
kjþ1
j¼1
ðjjCj Cjj jjxj1 jj þ jjDj Djj jjUj jjÞ
Proof. We rewrite (29) as pffiffiffi 8 k ðvÞ ak ðvÞðk ðvÞ fðvÞÞ þ ak ek ðvÞ, > > pffiffiffi > < if jk ðvÞ ak ðvÞðk ðvÞ fðvÞÞ kþ1 ðvÞ ¼ > þak ek ðvÞj Mvk ðvÞ, ð67Þ > > : 0, otherwise, where ek ðvÞ ¼
pffiffiffi ^k ðk ðvÞ zk Þ: ðvÞðk ðvÞ fðvÞÞ w
Byffiffiffi noticing that f(v) is the root of the linear function p ðx fðvÞÞ, by GCT for the theorem it suffices to show mðki , Ti Þ 1 X ðsÞ aj ej ðvÞ ¼ 0, ð68Þ lim lim sup T!0 i!1 T j¼k i
s ¼ 1, 2, 3, 4
8Ti 2 ½0, T
for any convergent subsequence ki ðvÞ ! ðvÞ, where ð1Þ
ð2Þ
ð3Þ
ð4Þ
ek ðvÞ þ ek ðvÞ þ ek ðvÞ þ ek ðvÞ ¼ ek ðvÞ, ð1Þ
^k Þðk ðvÞ zk Þ, ek ðvÞ ¼ ðwk w pffiffiffi ð2Þ ek ðvÞ ¼ ð ðvÞ wk Þk ðuÞ, pffiffiffi ð3Þ ek ðvÞ ¼ wk yk ðvÞfðvÞ, ð4Þ
ð4Þ
ek ðvÞ ¼ ek ¼ wk k : By the second limit in (31) and the convergence of the second series in (57) it follows that (68) holds for s ¼ 3. Further, by Lemma 3 (the first series in (34)) (68) holds for s ¼ 4. We now show (68) for s ¼ 1. For this we estimate ^k wk j. Since 12 3 > 0, there is a $ > 0 jw such that 12 3 3$ > 0. By (24)(26) it follows that x^kþ1 xkþ1 ¼ Ckþ1 ð^ xk xk Þ þ ðCkþ1 CÞxk þ ðDkþ1 DÞUkþ1 :
and
N2 > 0:
By Theorem 1 jjCk Cjj and jjDk Djj are of order oðk1=2$ Þ, and by Lemma 4, a:s:
xk ! 0. k$ k!1
Uk a:s: ! k$ k!1
0,
Then we have
Sð , kÞ ¼
kþ1 X
kjþ1 o
j¼1
1
j1=22$
X kþ1
1=22$ k 1=22$ j k j¼1
X
1=22$ k 1 k j ¼ o 1=22$
kjþ1 k j¼0 1 0 k
X k 2 X C 1 B ¼ o 1=22$ @ þ A j
k j¼0 k
¼o
1
kjþ1
1=22$
j¼ 2 þ1
k kjþ1
k 1 1 2 þ1 þ 2 ¼ o 1=22$ k2 k 1 k
1 ¼ o 1=22$ : k Noticing that the first term on the right-hand side of (69) exponentially decays, we find that
jj x^k xk jj ¼ o
1
k1=22$
a:s:
Since j v^k vk j and jj x^k xk jj are of the same order, we further have
1 1 ^k wk j ¼ o 3 j v^k vk j ¼ o 3 jj x^k xk jj jw b bk
k 1 ð70Þ ¼ o 1=23 2$ : k
591
Identification for Wiener Systems
By the selection of Mk defined in (14) from (29) it follows that jk ðvÞj k , which together with (70) implies ! 1 1 X X 1 ^k wk Þk ðvÞj ¼ o ak jðw : k3=23 2$ k¼1 k¼1 ð71Þ By Lemma 4 and (70) we have 1 1 X X ^k wk Þyk j ¼ ak jðw
1 yk ^ ð w w Þ k k k$ k1$ k¼1 ! 1 X 1 ¼o < 1: k3=23 3$ k¼1
k¼1
ð72Þ Again by (70) it is seen that 1 X
and s, ki ¼ 1 þ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ
ð77Þ
as i ! 1 and T ! 0. By (76) and the convergence of the second series in (34) it follows that X s s, jþ1 aj wj j ¼ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ j¼k i
as i ! 1 and T ! 0. By the last limit in (31) and (59)(77) it is seen that X s s, jþ1 aj wj yj j¼ki
s X
s, jþ1 aj ðwj jyj j E½wj jyj j þ E½wj jyj jÞ
j¼ki
¼ OðTÞ
^k wk Þk j ak jðw
8s 2 ½ki , . . . , mðki , TÞ. From the recursion (74) it follows that
k¼1
¼
1 X k¼1 1 X k¼1
^k wk Þj jk j ak jðw
o
ð78Þ
ðjk j Ejk jÞ þ Ejk j k3=23 2$
< 1:
ð73Þ
Combining (71)–(73) leads to (68) for s ¼ 1. Finally, we prove (68) for s ¼ 2. For this we first show that there exists an 0 with P0 ¼ 1 such that for any ! 2 0 if ki ðvÞ is a convergent subsequence, then ^s ðs ðvÞ zs Þ, sþ1 ðvÞ ¼ s ðvÞ as w
ð74Þ
sþ1 ðvÞ ¼ s ðvÞ as ws ðs ðvÞ zs Þ ^s Þðs ðvÞ zs Þ þ as ðws w s X ¼ s, ki ki ðvÞ þ s, jþ1 aj wj zj j¼ki s X
þ
^j Þðj ðvÞ zj Þ, s, jþ1 aj ðwj w
j¼ki
which incorporating with (71)–(73) yields sþ1 ðvÞ ¼ ki ðvÞ þ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ: ð79Þ
and jjsþ1 ðvÞ ki ðvÞjj cT, s ¼ ki ,ki þ 1, ...,mðki ,TÞ ð75Þ provided i is large enough and T > 0 is sufficiently small, where c > 0 does not depend on ki . Consider the recursion (74) starting from ki ðvÞ. Let
i, j ¼ ð1 ai wi Þ . . . ð1 aj wj Þ,
i j,
j, jþ1 ¼ 1:
By the first limit in (31) and the convergence of the first series in (58) it is clear that s X
aj wj ¼ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ,
j¼ki
which implies log s, ki ¼ O
s X j¼ki
! aj wj
ð76Þ
This means that the algorithm (29)(30) has no truncation for s 2 ½ki , . . . , mðki , TÞ, when i is large enough and T > 0 is sufficiently small. Since Lemmas 1, 3, 4 and 8 are all proved a.s., the possibly exceptional set is with probability zero. This means that 0 may be taken with P0 ¼ 1 such that for all ! 2 0 (74)(75) are verified for large i and small T > 0. Let ki ðvÞ ! ðvÞ. Then we have mðki , Ti Þ 1 X ð2Þ aj ej ðvÞ T j¼k i mðki , Ti Þ 1 X aj ðj ðvÞ ðvÞÞðwj Ewj Þ T j¼k i mðk i , Ti Þ X 1 aj ðwj Ewj Þ þ TðvÞ j¼k i mðk i , Ti Þ X pffiffiffi 1 aj j ðvÞðEwj Þ: ð80Þ þ T j¼k i
592
X.-L. Hu and H.-F. Chen
By the convergencepofffiffiffi the second series in (58), (75), and limk!1 Ewk ¼ ðvÞ (Lemma 1) we see that the first term on the right-hand side of (80) tends to zero as i ! 1 and T ! 0, while its second term tends to zero as i ! 1 by the convergence of the first pffiffiffi series in (58). Finally, by (75) and limk!1 Ewk ¼ ðvÞ the last term in (80) also tends to zero as i ! 1 and T ! 0. Thus, (68) for s ¼ 2 is proved and meanwhile the proof of Theorem 2 is completed. & Remark 4. Although the estimation algorithms consisting of (15),(16),(26)–(30) do not apply to the feedback control system, they are not complicated and suitable for implementation. Fig. 2.
6.
Numerical Experiments
In this section we give some examples with and without complete satisfaction of conditions required in theorems. For a better comparison in all simulations the linear subsystem is taken to be the same: vkþ1 0:3vk þ 0:6vk1 ¼ uk þ 0:27uk1 0:5uk2 þ 0:79uk3 with noise k 2 N ð0, 0:22 Þ and with iid input fuk g, uk 2 N ð0, 1Þ. Take Mk ¼ 0:1 þ k0:1 , 0 ¼ 2, and ak ¼ 1=k. The inputs for all simulations are taken from the same totality consisting of 6000 data. The difference between examples to be given consists in the nonlinear function fðÞ and in the selection of in bk ¼ 1=k . In all figures shown below the solid lines denote the estimates while the dotted lines the true values. Example 1. Let the nonlinear function be 8 < v 13 v3 , jvj 1, fðvÞ ¼ 2=3, v > 1, : 2=3, v < 1: Take bk ¼ k0:16 , which means ¼ 0.16. Obviously, all conditions required in Theorems 1 and 2 are satisfied. Figures 2 and 3 demonstrate how the estimates given in Section 2 approach to the true values. Figure 2 is for the coefficients of the linear subsystem, and Fig. 3 is for the nonlinear function. The estimate for fðÞ is obtained by separating the interval ½-1:5,1:5 into 100 subintervals with equal length and the values of fðÞ at endpoints of the subintervals are estimated. As shown in the figures, the simulation is consistent with the theoretical analysis.
Fig. 3.
An interesting question is that what will happen if some conditions required in theorems are not met. We show this by examples. Example 2. In (14) it is required that 3 þ < 1=2, which is violated if ¼ 0:6 or ¼ 0:9. We take the same system as that discussed in Example 1 but with different . It is clear the change of does not effect the estimates for coefficients of the linear subsystem. Figures 4 and 5 demonstrate the estimates for the nonlinear function for ¼ 0:6 and ¼ 0:9, respectively. From the figures it is seen that the larger the worse the estimates. As a matter of fact, from (27) we observe that bk ¼ 1=k characterizes the concentrating rate of vk on v. Further, from (28) we see that in order v^k be also concentrated on v, bk ¼ 1=k should not decrease faster than vk v^k , which by Theorem 1 tends to zero not faster than 1=k0:5 . This explains why
should be small.
593
Identification for Wiener Systems
Fig. 4.
Fig. 6.
Fig. 5.
Fig. 7.
Example 3. We now consider the situation where fðÞ is discontinuous, i.e., H1 fails to hold. Let 1, v 0, fðvÞ ¼ 1, v < 0, and let the linear subsystem be the same as that in Example 1. Figure 6 corresponds to the case where the parameters of algorithms completely coincide with those in Example 1, while in Fig. 7 the only change is that ¼ 0:6 rather than 0.16. From these figures we find that the algorithms proposed in the paper work well even H1 does not hold, if is appropriately small.
7.
Concluding Remarks
The paper gives the recursive and strongly consistent estimates for the Wiener systems, when its linear
subsystem is an RTF. From the paper we see that in comparison with the FIR type linear subsystem the infinite impulse response makes a great complication in analysis even under a local Lipschitz-type condition H1 on the nonlinearity. Since Gaussian distribution is even, defined by (13) is zero for any even functions. Thus, Theorem 1 does not apply to identifying Wiener systems with even nonlinearity. Therefore, for further research it is important to solve the even function problem. It is also of interest to weaken the conditions on fðÞ and to seek the conditions guaranteeing nondegeneracy of . To accelerate the algorithms used in the paper is of interest too. It is also worth considering the possible extension to the MIMO case.
References 1. Bendat JS. Nonlinear system analysis and identification from random data. New York Wiley 1999
594 2. Billings SA, Fakhouri FY. Identification of nonlinear systems using the Wiener model. Electron Lett 1978; 13: 502–504 3. Chen HF. Stochastic approximation and its applications. Kluwer, Dordrecht 2002 4. Chen HF. Pathwise convergence of recursive identification algorithms for Hammerstein systems. IEEE Trans Autom Control 2004; 49(10): 1641–1649 5. Chen HF. Recursive identification for Wiener model with nonlinearity being discontinuous piece-wise linear function. IEEE Trans Autom Control 2006; 51(3): 390–400 6. Chow YS, Teicher H. Probability theory. Springer, New York 1978 7. Greblicki W. Nonparametric approach to Wiener system identification. IEEE Trans Circuits Syst-I: Fundam Theory Appl 1997; 44(6): 538–545 8. Greblicki W. Recursive identification of Wiener system. Int J Appl Math Comp Sci 2001; 11(4): 977–991 9. Hasiewicz Z. Identification of a linear system observed though zero-memory nonlinearity. Int J Syst Sci 1987; 18: 1595–1607 10. Hu XL, Chen HF. Strong consistency of recursive identification for Wiener systems. Automatica 2005; 41(11): 1905–1916 11. Hu XL, Chen HF. Recursive identification for Wiener systems using Gaussian signals. to appear in Asian J Control 2007.
X.-L. Hu and H.-F. Chen
12. Hunter IW, Korenberg MJ. The identification of nonlinear biological systems: Wiener and Hammerstein cascade models. Biol Cybern 1986; 55: 135–144 13. Krzyzak A, Sasiadek JZ, Kegl B. Non-parametric identification of dynamic non-linear systems by a Hermite Series Approach. Int J Syst Sci 2001; 32(10): 1261–1285 14. Pajunen GA. Adaptive control of Wiener type nonlinear systems. Automatica 1992; 28: 781–785 15. Nordsjo¨ AE, Zetterberg LH. Identification of certain time-varying nonlinear Wiener and Hammerstein systems. IEEE Trans Signal Proc. 2001; 49: 577–592 16. Verhaegen M, Westwick D. Identifying MIMO Wiener systems in the context of subspace model identification methods. Int J Control 1996; 63(2): 331–349 17. Vo¨ro¨s J. Parameter identification of Wiener systems with discontinuous nonlinearities. Syst Control Lett 2001; 44: 363–372 18. Westwick DT, Kearney RE. A new algorithm for the identification of multiple input Wiener systems. Biol Cybern 1992; 68: 75–85 19. Wigren T. Recursive prediction error identification using the nonlinear Wiener mode. Automatica 1993; 29: 1011–1025 20. Wigren T. Convergence analysis of recursive algorithms based on the nonlinear Wiener model. IEEE Trans Automat Control 1994; 39: 2191–2206
Identification for Wiener Systems with RTF Subsystems Xiao-Li Hu and Han-Fu Chen** Institute of Systems Science, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, P. R. China
Recursive algorithms are proposed for identification of the Wiener system composed of a subsystem with rational transfer function (RTF) cascaded with a static nonlinearity. The nonlinear block is not necessary to be invertible. The strong consistency is proved for estimates for the coefficients of RTF and for the value of the nonlinear function at any point. Some numerical examples are given showing that the simulation results are consistent with theoretical analysis. Keywords: Recursive identification; Wiener systems; rational transfer function; strong consistency; stochastic approximation
1.
Introduction
The Wiener system consisting of a linear subsystem followed by a static nonlinearity is an important model in practice. For recent years, its identification issue has been attracting a great attention from researchers, e.g., [2,4,7,8,10–20] among others. The SISO Wiener system considered in the paper is presented by Fig. 1. To present the nonlinearity there are the parametric [1,5,9,12,14–20] and the nonparametric approaches [2,7,8,10,11,13]. The present paper belongs to the latter and is the continuation of [10,11]. If vk were observed, then identification of the linear subsystem would be an easy task. Therefore, it is important to extract information concerning vk from observations. In the case where fðÞ is invertible, vk can Supported by the NSFC under Grants No. G0221301, 60334040, 60474004. Correspondence to: H.-F. Chen, e-mail: [email protected]
be expressed by system output: vk ¼ f1 ðyk Þ, and such an expression may simplify analysis in many cases, in particular, for the case where the linear subsystem is also nonparametric [7,8,13]. However, the invertibility condition on fðÞ excludes many important nonlinearities including saturations from consideration. This motivated us to consider identification of Wiener systems without requiring invertibility of fðÞ in [10,11], where the recursive identification algorithms and their strong consistency were established for Wiener systems with general nonlinearities. However, the linear subsystem considered there is limited to an finite impulse response (FIR) system. In this paper the FIR system is extended to the RTF system, which clearly is with infinite impulse response. Similar to [10,11], the parameters of the linear part are first estimated, then the intermediate signals vk are approximated by using the estimated linear subsystem. Finally, the nonlinear block is identified on the basis of the approximated intermediate signals and the observations fzk g. Such an extension is technically nontrivial by the following observation. In the convergence analysis itP is of importance to have convergence of series like 1 k¼1 ak vk where fak g are the stepsizes used in the algorithms to be defined later on. P For an FIR linear subsystem the series 1 k¼1 ak vk can be split into a finite number of convergent subseries (due to independence), while for an RTF subsystem such a direct separation is impossible and the treatment is more delicate. The estimates for coefficients of the linear subsystem and for the value fðvÞ of the static nonlinearity Received 10 July 2006; Accepted 15 January 2007 Recommended by A. Karimi and L. Ljung
582
X.-L. Hu and H.-F. Chen
Fig. 1. Wiener system.
at any v are given by the stochastic approximation (SA) algorithms with expanding truncations [3]. These algorithms together with conditions for their convergence are presented in Section 2. The limiting behaviors of some correlation functions and series are demonstrated in Section 3. In order to deal with the difficulty arising from the infinite impulse response of the linear subsystem mentioned above, the totality of integers is separated into a countable union of disjoint subsets, and convergence of several series is established with the help of such a separation. This is done in Section 4. The strong consistency of the estimates given in Section 2 is proved in Section 5, and a few numerical examples are given in Section 6. In the last section some concluding remarks are given.
2. 2.1.
Problem Setting and Estimation Algorithms Problem Setting
Let the linear subsystem be with RTF CðzÞvkþ1 ¼ DðzÞuk ,
ð1Þ
where CðzÞ ¼ 1 þ c1 z þ . . . þ cp zp ,
ð2Þ
DðzÞ ¼ 1 þ d1 z þ . . . þ dq zq
ð3Þ
are polynomials in backward-shift operator zuk ¼ uk1 . It is noticed that in the special case CðzÞ ¼ 1, the linear subsystem is with FIR and is considered in [10,11]. The purpose of this paper is to estimate the coefficients ci , i ¼ 1, . . . , p and dj , j ¼ 1, . . . , q and the value fðvÞ of the nonlinearity at any v on the basis of the designed inputs fuk g and observations fzk g: zk ¼ yk þ " k ,
yk ¼ fðvk Þ,
where f"k g is the observation noise. We now list conditions to be used in the sequel.
ð4Þ
H0 CðzÞ and DðzÞ are coprime and CðzÞ is stable. H1 The function fðÞ satisfies the following condition jfðxÞ fðyÞj < Lð1 þ jxj þ jyj Þjx yj 8x,y 2 R; ð5Þ where L > 0 and 0 are some constants. H2 fk g is a sequence of mutually independent random variables with Ek ¼ 0 and supk E2k < 1, and is independent of fuk g. Remark 1. Condition H1 proposed here is probably used for the first time. H1 is more general than the global Lipschitz condition, which requires ¼ 0 in (5), but it is stronger than the local Lipschitz condition, by which in lieu of Lð1 þ jxj þ jyj Þ in (5) there may be any continuous function of x, y. Under H1 the growth rate of fðzÞ is not faster than a polynomial as z tends to infinity. This property is used in (40) to guarantee an acceptable approximation to yk whenever vk is well estimated, which is needed for deriving convergence of some stochastic series in Lemma 8. Remark 2. We do not claim any necessity of H1. As shown later (Section 6, Example 3), even if H1 does not hold, the identification algorithm given in the paper still provides satisfactory estimates for both linear and nonlinear subsystems. However, it is not clear if it works in a feedback case. Since fuk g is at designer’s disposal, let uk ¼ 0, 8k 0. Then under Condition H0, (1) can be expressed as k1 X hj ukj , h0 ¼ 1, ð6Þ vkþ1 ¼ C1 ðzÞDðzÞuk ¼ j¼0
where the coefficients hj exponentially decay: jhj j ¼ O ej , > 0: 2.2.
ð7Þ
Stochastic Approximation Method
Since our identification algorithms are based on an improved version of classical stochastic approximation (SA) algorithm, we first briefly introduce SA and its improved version.
583
Identification for Wiener Systems
For an unknown function gðÞ with single root denoted by x0 , the classical Robbins-Monro (RM) algorithm, being the first version of SA, suggests to seek the root x0 by the simple recursion: xkþ1 ¼ xk þ ak gkþ1 ,
k ¼ 0, 1, . . . ,
ð8Þ
where
gkþ1 ¼ gðxk Þ þ kþ1
ð9Þ
is the noisy observation of gðxk Þ, xk is the k-th estimate for x0 , kþ1 is the observation noise, and fak g is a sequence of positive real numbers tending to zero as k ! 1. For convergence of xk to x0 , a certain growth rate restriction on gðÞ or an a priori boundedness assumption on fxk g are required in the classical theory in addition to conditions on the noise fk g. These requirements quite restrict the application of SA. Moreover, the noise is usually required to analyze for the whole sequence of estimates fxk g. In order to extend the range of applicability of SA method, the RM algorithm (8) is modified to be truncated at expanding bounds in [3]. By this the increasing rate of fxk g is controlled. Let fMk g be an increasing sequence of positive real numbers with Mk ! 1, and let xk be generated by k!1
the following algorithm: xk þ ak gkþ1 , if jjxk þ ak gkþ1 jj Mk , xkþ1 ¼ otherwise, x ; ð10Þ k ¼ k1 þ I½jxk1 þak1 gk j>Mk1 ,
0 ¼ 0,
ð11Þ
where function IA denotes the indicator function of a set A, and k is the number of truncations occurred until time k 1. From (11) it is seen that k k 1, and from (10) jjxkþ1 jj jjx jj ^ Mk1 . This means that the growth of jjxkþ1 jj should not be faster than that of Mk . At any time if jjxk þ ak gkþ1 jj exceeds the truncation bound Mk , then we pull xkþ1 back to the fixed point x and simultaneously extend the truncation bound from Mk to Mk þ1 . Otherwise, it develops as (8). Since the region where x0 is located is unknown, it is important to allow xk to grow up in order to have possibility to reach x0 . For convergence of fxk g we must have ak gðxk Þ ! 0. Thus, it is reasonable to require k!1
ak ! 0. However, if ak tends to zero too fast so that Pk!1 1 may be unable to k¼1 ak < 1, then the estimate xkP reach x0 . Thus, we have to require 1 k¼1 ak ¼ 1 (see 1 [3]). A common choice is ak ¼ k. The selection of
fMk g is defined by that how fast we allow fxk g to grow up. The following fact to be used in what follows is a special case of the general convergence theorem (GCT) Theorem 2.2.1 of [3]. Let gðxÞ ¼ ðx x0 Þ and ak ¼ 1=k. Then xk defined by (10)(11) converges to x0 for any sample paths satisfying mðnX k , Tk Þ ai iþ1 ¼ 0 8Tk 2 ½0, T lim lim sup T1 i¼n T!0 k!1 k
ð12Þ for any nk such that fxnk g converges, where m X aj T : mðk, TÞ ¼ max m : j¼k
Further, if fk g in (9) for the sample path ! under consideration can be decomposed into two parts k ¼ k0 þ k00 such that 1 X
0 a1 k k < 1,
k00 ¼ Oðak Þ
k¼1
for some 2 ð0, 1, then by the convergence rate theorem (CRT) Theorem 3.1.1 of [3] xk given by (9)(10)(11) converges to x0 with the following rate: k xk x0 k ¼ oðak Þ: It is worth noting that Condition (12) is required to verify only along convergent subsequences fxnk g rather than along the whole sequence fxn g. As pointed in [3], in many cases (12) cannot be directly verified along the whole sequence fxn g, but it can be done for convergent subsequences. This is the case in Theorem 2 of Section 5. The method verifying Condition (12) along convergent subsequences has shown a great advantage over verification along the whole sequence, and it is called as trajectory-subsequence (TS) method in [3]. 2.3.
Identification Algorithms
As stated in introduction, we first estimate the parameters ci , i ¼ 1, . . . , p, dj , j ¼ 1, . . . , q of the linear part, then construct v^k to approximate the output vk of the linear subsystem; finally, identify the nonlinear block by using the data f^ vk , zk g. We take the input sequence fuk g1 1 to be iid and uk 2 N ð0, 1Þ. Under condition H1, let us set Z 1 X t2 1 2 2 hi , and ¼ pffiffiffiffiffiffi tfðtÞeð22 Þ dt: ¼ 2 3 R i¼0 ð13Þ
584
X.-L. Hu and H.-F. Chen
P 2 Clearly, vk obeying N ð0, k2 i¼0 hi Þ is asymptoti2 cally distributed as N ð0, Þ, and is the asymptotic correlation coefficient of the input and output of the nonlinearity. Define the real sequences to be used in algorithms: ak ¼ k1,
bk ¼ k1 ,
Mk ¼ M0 þ k
and
ð14Þ
with > 0, > 0, M0 > 0, and 3 þ < 12. Similar to [4,10,11], hi , i ¼ 0, 1, . . . , p þ q are estimated by the SA algorithms with expanding truncations: kþ1 ðiÞ ¼ ½ k ðiÞ ak ð k ðiÞ uk zkþiþ1 Þ I½j k ðiÞak ð k ðiÞuk zkþiþ1 ÞjM k ðiÞ ,
ð15Þ
Then (19) can be written in the vector form ½hqþ1 hqþ2 . . . hqþp T ¼ ½c1 c2 . . . cp T ,
which is called the Yule–Walker equation. Replacing hi , i ¼ q, . . . , p þ q, in (21) with their estimates ðiÞ k if k ð0Þ 6¼ 0 hik ¼ k ð0Þ 0 otherwise and denoting by k the resulting matrix corresponding to , we derive a linear equation with respect to ½c1 . . . cp T at time k. The solution ½c1k . . . cpk T whenever exists to the obtained linear equation serves as an estimate of ½c1 . . . cp T , i.e., T ½c1k . . . cpk T ¼ 1 k ½hqþ1, k . . . hpþq, k :
k ðiÞ ¼ k1 ðiÞ þ I½j k1 ðiÞak1 ð k1 ðiÞuk1 zkþi Þj>M k1 ðiÞ , ð16Þ
0 ðiÞ ¼ 0
E½vkþlþ1 uk ¼ hl ,
l ¼ 0, 1, . . . :
ð17Þ
Multiplying both sides of (1) by uki and taking the mathematical expectation lead to hi ¼ E½vkþ1 uki ¼
i^p X
cl hil þ di ,
i ¼ 1, . . . , q,
l¼1
hj ¼ E½vkþ1 ukj ¼
i^p X
clk hil, k ,
i ¼ 1, . . . , q:
Thus, we have estimated all coefficients of the linear subsystem. It remains to estimate fðvÞ at any fixed v. Let us use these estimates to approximate the intermediate signal vk first. Define 0 1 1 0 0 c1 B c 0 1 0 C B C 2 B C B C C ¼ B C , B C 0 1 A @ cp1 0 cp 0 0 0 pp 0 1 1 d 1 dq B 0 0 0 C B C D¼B , C @ A 0
0
Uk ¼ ½uk . . . ukq T ,
0
pðqþ1Þ
H ¼ ½1 0 . . . 01p
xkþ1 ¼ Cxk þ DUkþ1 :
cl hjl ,
ð23Þ
l¼1
ð18Þ j^p X
ð22Þ
Motivated by (18), the estimates for di are given by dik ¼ hik þ
with arbitrary initial values 0 ðiÞ, i ¼ 0, 1, . . . , p þ q. We note that k ð0Þ serves as the estimate for , and ð k ðiÞ= k ð0ÞÞ if exists serves as the estimate for hi , i ¼ 1, . . . , p þ q. In order to estimate the coefficients of the linear subsystem let us first derive the corresponding Yule– Walker equation. Since fuk g is iid with Euk ¼ 0, from (6) it follows that
ð21Þ
ð24Þ
l¼1
j ¼ q þ 1, . . . , p þ q,
ð19Þ
where and hereafter a ^ b denotes minða, bÞ. Define 0
hq h B qþ1 ¼B @ hpþq1
hq1 hq hpþq2
1 hqpþ1 hqpþ2 C C, A hq
where hi ¼ 0 for i > 0 by definition.
It is clear that with x0 appropriately chosen the output vk of the linear subsystem coincides with Hxk : vk ¼ Hxk :
ð20Þ
ð25Þ
Replacing the elements of matrices C and D with their estimates derived above gives Ck and Dk which naturally serve as estimates for C and D, respectively. Define x^kþ1 ¼ Ckþ1 x^k þDkþ1 Ukþ1 ,
v^k ¼ H x^k
ð26Þ
with an arbitrary x^0 . We use v^k as an estimate for vk .
585
Identification for Wiener Systems
For estimating fðvÞ we need the kernel function
wk ¼
pffiffiffiffi 2 bk exp
2 vkbv þ k
v2k 220
ð27Þ
,
where 0 > 0 is a fixed constant. The kernel function wk can be estimated by
^k ¼ w
pffiffiffiffi 2 bk exp
(
v^k v 2 þ bk
) v^2k 220
:
ð28Þ
Finally, fðvÞ is estimated by k ðvÞ given by the following recursive algorithm:
Consequently, we have lim Ewk yk
k!1
Z pffiffiffiffi
x2 22 0
x2
2 ðxv bk Þ þ
1 22 k fðxÞpffiffiffiffiffiffi e dx 2 k (
) Z ðvþbk tÞ2 1 1 t2 1 dt ¼ lim e fðvþbk tÞk exp k!1 R 2 20 2k pffiffiffi ¼ fðvÞðvÞ: &
¼ lim b2 e k!1 k R
Lemma 2. Assume H1 holds. Then lim E½ukj ykþ1 ¼ hj ,
j ¼ 0, 1, . . . , p þ q,
k!1
^k ðk ðvÞ zk Þ kþ1 ðvÞ ¼ ½k ðvÞ ak w I½jk ðvÞak w^k ðk ðvÞzk ÞjMk ðvÞ ,
ð29Þ
k ðvÞ ¼ k1 ðvÞ þ I½jk1 ðvÞak1 w^k1 ðk1 ðvÞ yk1 Þj > Mk1 ðvÞ ,
0 ðvÞ ¼ 0 ð30Þ
with any given initial value 0 ðvÞ.
3.
Limiting Behaviors of Some Random Variables
We now establish the limiting behaviors of some random variables connected with the kernel and the a.s. convergence of some series. Lemma 1. Assume H1 holds. For wk defined by (27) the following limits take place: pffiffiffi pffiffiffi ðvÞ, lim E½wk yk ¼ fðvÞðvÞ, k!1 k!1 pffiffiffi ð31Þ lim E½wk jyk j ¼ jfðvÞjðvÞ, lim Ewk ¼
k!1
pffiffiffiffiffi pffiffi lim E½ bk wk 2 ¼ 2ðvÞ, k!1 pffiffiffiffiffi pffiffi lim E½ bk wk yk 2 ¼ 2fðvÞðvÞ,
where and are defined by (13). Proof. Noticing that ukj and vkþ1 are jointly Gaussian and by (6) h i h E ukj 2j vkþ1 vkþ1 ¼ 0, k
h
we find that ukj 2j vkþ1 is independent of vkþ1 . k Taking the conditional expectation with respect to vkþ1 leads to
hj E ukj jvkþ1 ¼ 2 vkþ1 : k By this we have
E ukj fðvkþ1 Þ ¼ E E ukj fðvkþ1 Þjvkþ1
¼ E fðvkþ1 ÞE ukj jvkþ1 ¼ hj
Lemma 3. Assume H2 holds. Then the following series converge a.s.:
n 2 v 2
1 0
12
ak w k k ,
k¼1 1 X
k!1
where ðvÞ ¼ 1 exp
E½vkþ1 fðvkþ1 Þ , 2k
which by (13) yields the desired result by letting k ! 1 because 2k ! 2 and vk is asymptotically & distributed as N ð0, 2 Þ.
1 X
ð32Þ
ð33Þ
1 X k¼1
ak wk ðjk j Ejk jÞ, ð34Þ
ak uk kþj :
k¼1
o .
Proof. It suffices to show limk!1 Ewk yk , and the rest can be provedP in a similar manner. k2 hj ukj1 , and hence vk 2 N ð0, 2k Þ By (6) vk ¼ j¼0 P k2 2 with 2k ¼ j¼0 hj .
Proof. The proof is essentially based on the following fact (see Chapter 11.1, Theorem 3 of [6]): Let fSn , F n , n 1g be a martingale with supn EjSn j < 1, and let fYn , F n1 , n 1g be with supn jYn j < 1 a.s. 1 P Xn Yn < 1 a.s. where Xn ¼ Sn Sn1 . Then n¼1
586
X.-L. Hu and H.-F. Chen
For example, ak wk k can be written as ak wk k ¼ Xk Yk 2 . Since and Yk ¼ exp vkbv k
pffiffiffiffi 2 with Xk ¼ k1
k
0 < < 1/6, applying the above cited result leads to 1 P Xn Yn < 1 a.s. n¼1
The proof for convergence of the second series is similar. The convergence of the last series follows from the fact that fuk kþj g is a martingale difference sequence. In what follows we need the divergence rates of uk and yk as k ! 1. Lemma 4. For any > 0 the following rate takes place uk a:s: ! 0: k k ! 1
ð35Þ
If C(z) is stable and H1 holds, then k a:s: ! 0, k k ! 1
yk a:s: ! 0: k k ! 1
vk a:s ! 0, k k!1
1 X juk j P > " < 1: k k¼1
4.
&
Convergence Analysis via Separation of Integer Set
Since the linear subsystem is with RTF, its output is no longer a combination of inputs of the system. This makes convergence analysis for some series difficult. To overcome the difficulty we truncate the expression (6) for vkþ1 as follows vkþ1 ¼
k X i¼0
hi uki ,
k ¼ ½k , 0 < < 1,
ð39Þ By (5) of Condition H1 and Lemma 4, we further have vkþ1 Þj jykþ1 ykþ1 j ¼ jfðvkþ1 Þ fð < Lð1 þ jvkþ1 j þ j vkþ1 j Þjvkþ1 vkþ1 j ¼ Lð1 þ jvkþ1 j þ j vkþ1 j ÞO e2k a:s:, ð40Þ ¼ O e 4 k and similarly,
which incorporating with H1 implies
yk a:s ! 0. k k!1
Let us explain the idea of truncations applied in (37). Assume ¼ 1=2. Then vk2 þ1 is a linear combination of fuk2 k , uk2 kþ1 , . . . , uk2 g for any positive integer k. By the independence of the sequence fuk g and the fact k2 k > ðk 1Þ2 , it is clear that the subsequence f vk2 þ1 , k ¼ 0, 1, . . .g now turns to be a sequence of mutually independent random variables. On the other hand, when the linear subsystem is stable, the slowly truncated vkþ1 approximates vkþ1 as k ! 1 based on the following analysis. By (7) and (35) it is clear that X k1 hi uki ¼ O e2k jvkþ1 vkþ1 j ¼ a:s:
ð36Þ
Then by the Borel–Cantelli lemma (35) follows. By stability of C(z) from (24)(25) and (35) we have
ð38Þ
i¼kþ1
Proof. For any > 0 by the Tchebychev inequality " # juk j juk j2= 1 P >" ¼P > "2= < 2= 2 Ejuk j2= , k2 k " k which implies
where and hereafter ½x denotes the integer part of a real number x. Correspondingly, we set (
) pffiffiffiffiffiffi vk v 2 v2k 2 k ¼ yk ¼ fð vk Þ, w exp þ 2 : bk bk 20
kþ1 j ¼ O e4k jwkþ1 w
a:s:,
ð41Þ
2kþ1 j ¼ O e4k jw2kþ1 w
a:s:
ð42Þ
PTo analyze convergence of the series like ak ðuk ykþ1 Euk ykþ1 Þ, we now separate the totality of integers into a countable union of disjoint sets. Let f k g be a sequence of strictly increasing integers such that 0 ¼ 0, 1 ¼ 1, f kþ1 k g also nondecreasingly diverges to infinity, and 0 < lim inf k!1
k ka
lim sup k!1
k ka
<1
for some a > 1: ð43Þ
Define ð37Þ
A0 ð0Þ ¼ f k : k ¼ 1, 2, . . .g;
ð44Þ
587
Identification for Wiener Systems
Aj ðiÞ ¼ f k þ i : k ¼ j, j þ 1, . . .g,
Aj ¼
i[ 2 ðjÞ
Aj ðiÞ,
ð45Þ
j ¼ 1, 2, . . . ,
l ðj, iÞ ¼
i¼i1 ðjÞ
where i runs from i1 ðjÞ ¼ j j1 to i2 ðjÞ ¼
jþ1 j 1. If i2 ðjÞ < i1 ðjÞ, then set Aj ¼ . We denote by J the totality of integers f1, 2, . . .g. Lemma 5. The integer sets Aj ðiÞ, j ¼ 0, 1, . . . , compose a disjoint separation of J, i.e., Aj ðiÞ, i ¼ i1 ðjÞ, . . . , i2 ðjÞ, j ¼ 0, 1, . . ., are disjoint and A0 ð0Þ [
1 i[ 2 ðjÞ [
Aj ðiÞ ¼ J:
Set
ð46Þ
jas1 1=2 s j1þl þi :
j1þl þ i
By (43) we have ! jas1 1=2 jas1 1=2 1 s < ¼O
sj1þl
j1þl þ i ðj 1 þ lÞ1 þ1=2
1 ð50Þ ¼ O þ1=2 , l1 which combining with (49) yields El ðj, iÞ ¼ 0 and
j¼1 i¼i1 ðjÞ
El2 ðj, iÞ < Proof. The assertions of the lemma follow from the following set operations by noticing that f k þ 1, . . . , kþ1 1g, k ¼ 1, 2, . . ., are disjoint: 1 i[ 2 ðjÞ [
Aj ðiÞ ¼
j 1 [ 1 jþ1[ 1 [
f k þig
j¼1 i¼ j j1 k¼j
j¼1 i¼i1 ðjÞ
¼
j 1 1 [ 1 jþ1[ [
f k þig
¼
C l1þ21
ð51Þ
for some constant C > 0. By the condition of the lemma for fixed j and ifl ðj, iÞg is a sequence of mutually independent random variables with possible exception of a finite number of terms. Then, by the convergence theorem for martingale difference sequences (mds)
j¼1 k¼j i¼ j j1 1 [ k [
Sji ðmÞ ¼
f k þ j j1 ,
k þ j j1 þ1, ..., k þ jþ1 j 1g 1 [ & ¼ f k þ1, k þ2, ..., kþ1 1g: Lemma 6. Let fk g be a sequence of random variables with zero mean and supk Ek2 < 1. If for any fixed i and j, i ¼ i1 ðjÞ, . . . , i2 ðjÞ whenever i2 ðjÞ i1 ðjÞ, j ¼ 1, 2, . . . , the subsequence fk : k 2 Aj ðiÞg is composed of mutually independent random variables with possible exception of a finite number of k , then
for s > 3 2
1 a,
l ðj, iÞ
a.s. converges to some random variable Sji as m ! 1. Given > 0 by (50)(51) we have
P jSji ðmÞj > j
k¼1
a:s:
m X l¼1
k¼1 j¼1
1 X 1 <1 s k k k¼1
ð49Þ
ð47Þ
where a is given in (43).
3 1 1 a1 Proof. From a1 s > 2 a 1it is clear that as a þ 2 > 2 . Let 2 2 , as a þ 2 and let 1 > 0 such that 1 ð48Þ as a þ 1 > 0: 2 Since s > 12, under the conditions of the lemma it is P 1 clear that ks k < 1 a.s. Consequently, the k2A0 ð0Þ
index set A0 ð0Þ can be excluded from consideration.
1 2 j2
EjSji ðmÞj2
1 C X 1 C0 ¼ 2 2 , 2 2 1þ2 1 j j l¼1 l
P 1 where C0 ¼ C 1 l¼1 l1þ21 . Letting m ! 1 we have
P jSji j > j C2 j02 :
ð52Þ
Define j ðÞ¼ ! : jSji j j , i ¼ i1 ðjÞ, . ..,i2 ðjÞ , ð53Þ 1 \ j ðÞ: ðÞ¼ j¼1
In the case a 2 it is clear that i2 ðjÞ i1 ðjÞ 1 ¼ jþ1 2 j þ j1 ¼ Oððj þ 1Þa 2ja þ ðj 1Þa Þ ¼ Oðja2 Þ:
588
X.-L. Hu and H.-F. Chen
Thus, similar to (55) for ! 2 j ðÞ by Lemma 5 it follows that
By (52)(53) it follows that !
1 [
PððÞÞ ¼ 1 P
cj ðÞ
>1
j¼1
>1
1 X j¼1
1 X P cj ðÞ j¼1
i2 ðjÞ X 1 1 X X 1 jSji j k s as k¼1 k j¼1, A 6¼ j 1 1=2 i¼i ðjÞ j
C0 O ja2 2 2 j
¼
! 1 1 X 1 >1O 2 : j¼1 j2þ2a
j¼1, Aj 6¼
ð54Þ
By the definition of , we have that 2 þ 2 a > 1 and, hence, PððÞÞ ! 1. !1
For ! 2 ðÞ by Lemma 5 it follows that i2 ðjÞ X X 1 1 X 1 X 1 1 ¼ k s jþl1 þi s k¼1 k j¼1 ð jþl1 þ iÞ i ðjÞ l¼1
¼ O ð 1Þ
j¼1 i¼i1 ðjÞ
¼O
!
jas1 ða2Þ1=2
i
< 1:
< 1:
!1
Therefore,
0
>1
jas1 1=2
j¼1
i¼1
Similar to (56), for j ðÞ defined by (53) we have ! : jSji j kj , , if i2 ðjÞ ¼ i1 ðjÞ ¼ i, j ðkÞ ¼ otherwise:
1 [
1 X
j¼1, Aj 6¼ 1 X
1 cj ðÞA
P cj ðÞ
C0 2 j2 j¼1, Aj 6¼ 1 X 1 ¼ 1 O ð 1Þ 2 ! 1: 2 ia1 !1 i¼1
>1
1
1 1 a1 ðas1 1=2Þ
j¼1, Aj 6¼
jSji j
1 X
1 X
PððÞÞ ¼ 1 P@
X i2 ðjÞ X 1 1 X l ðj, iÞ ¼ as 1=2 j¼1 i1 ðjÞ l¼1 j 1
Oð1Þ jas1 1=2
To complete the proof it remains to show PððÞÞ ! 1.
1
i2 ðjÞ 1 X X
1
1 X
&
ð55Þ
Lemma 7. Assume C(z) is stable. Let s 2 ð12, 32Þ, P2 ð0, s 12Þ, and vk be defined by (37). Then 1 ks vk < 1 a.s.
In the case 1 < a < 2 we note that for sufficiently large j, ð jþ1 j Þ ð j j1 Þ equals either 0 or 1. In other words, i2 ðjÞ i1 ðjÞ equals either 1 or 0, and correspondingly,
Proof. Since fuk g is iid and Euk ¼ 0, we have E vk ¼ 0 1 2 , 32s Þ. It is clear that and supk v2k < 1. Take a 2 ð1 3 1 s > 2 a, which is required in Lemma 6. It suffices to show that f vk : k 2 Aj ðiÞg is mutually independent. If this is done, then the assertion of the lemma follows from Lemma 6. By (37) v½ðkþ1Þa þi is a linear combination of u½ðkþ1Þa þi , . . . , u½ðkþ1Þa þi½ðkþ1Þa þi , where by (31)
Aj ¼
i[ 2 ðjÞ
Aj ðiÞ ¼
i¼i1 ðjÞ
Aj ðiÞ; if i2 ðjÞ ¼ i1 ðjÞ ¼ i, , otherwise: ð56Þ 1
By noticing that i1 ðjÞ ¼ i implies j ¼ Oðia1 Þ, for ! 2 j ðÞ we have jSji j as j 1 1=2
! 0:
jas1 1=2 j!1
½ðk þ 1Þa þ i ¼ ½ð½ðk þ 1Þa þ iÞ . The total number of these terms is ½ðk þ 1Þa þ i þ 1 ¼ Oðka Þ. Let us denote by Uk the minimal subset of fuj , j ¼ 1, 2, . . .g such that v k þi is the linear combination of elements from Uk. On the other hand, the neighboring terms in f vk , k 2 Aj ðiÞg are with subscripts kþ1 þ i and k þ i, and ð kþ1 þ iÞ ð k þ iÞ ¼ Oðka1 Þ. This means that
589
Identification for Wiener Systems
starting from some large enough k0 fUk , k ¼ k0 , k0 þ 1, . . .g are disjoint subsets of fuj , j ¼ 1, 2, . . .g whenever a < a 1. Therefore, f vk : k 2 Aj ðiÞg is a sequence of mutually independent random variables with possible exception of finite number of terms. Thus, the conditions of Lemma 6 are satisfied. Lemma 8. Assume C(z) is stable and H1 holds. Then the following series converge a.s. 1 X ak ðukj ykþ1 Eukj ykþ1 Þ 8j 0, k¼1 1 X
ð57Þ ak ðwk yk Ewk yk Þ,
k¼1 1 X
k ð0Þ ! k!1
a:s:,
k ðiÞ ! hi k!1
ð61Þ
a:s:,
i ¼ 1, . . . , p þ q with convergence rates j k ð0Þ j ¼ oðk Þ a:s:,
j k ðiÞ hi j
ð62Þ
a:s
ð58Þ ak ðjwk Ewk j Ejwk Ewk jÞ,
in addition, 6¼ 0 and det 6¼ 0, then cik ! ci , k!1
a:s
i ¼ 1, . . . , p, dik ! dj , j ¼ 1, . . . , q with the same k!1
rate as that given in (62), where cik and djk are given by (22)(23), respectively. ak ðwk jyk j E½wk jyk jÞ:
ð59Þ
k¼1
Proof. Let us prove convergence of the first series in (57). The rest can be done by a similar treatment. By (35)(40) and H1 it is clear that 1 X k¼1 1 X
ak ukj ðykþ1 ykþ1 Þ < 1 a:s:,
and
ak ½Eukj ðykþ1 ykþ1 Þ < 1 a:s:
k¼1
Thus, for the first series in (57) it suffices to show 1 X
Theorem 1. Assume H0-H2 hold. Then
where 0 < < 1=2 and is given in (13). Further, if,
ak ðwk Ewk Þ,
k¼1 1 X
We now in a position to establish the strong consistency of estimation algorithms given in Section 2.
¼ oðk Þ a:s:, i ¼ 1, . . . , p þ q,
k¼1 1 X
5. Convergence of Identification Algorithms
ak ðukj ykþ1 Eukj ykþ1 Þ < 1 a:s:
ð60Þ
k¼1
Notice that Eðukj ykþ1 Eukj ykþ1 Þ ¼ 0 and sup Ejukj ykþ1 Eukj ykþ1 j2 < 1: k
Similar to Lemma 7, separating fuj , j ¼ 1, 2, . . .g to the disjoint subsets Uk with s ¼ 1 and 2 ð0, 1=2Þ, we find that fukj ykþ1 Eukj ykþ1 , k 2 Aj ðiÞg is a sequence of mutually independent random variables with possible exception of a finite number of terms. Then (54) follows from Lemma 6. & Remark 3. As shown in Lemma 6, the stepsize ak in (57)-(59) may be replaced by k1s with s > 12.
Proof. Rewrite (15) as 8 k ðiÞak ð k ðiÞhi Þak kþ1 ðiÞ, > > > < if j ðiÞa ð ðiÞh Þa ðiÞj k k k i k kþ1 kþ1 ðiÞ¼ > M ðiÞ,
k > > : 0, otherwise, where kþ1 ðiÞ ¼ uk zkþiþ1 þ hi ,
ð63Þ
i ¼ 0, 1, . . . , p þ q:
ð64Þ In the sequel we will use the general convergence theorem (GCT) and the convergence rate theorem (CRT) for the SA algorithm with expanding truncations mentioned in Section 2.2 and cited in [11] with reference to [3] for detailed proof. From (63) it is seen that by GCT with regression functions gðiÞ ðxÞ ¼ ðx hi Þ for (61) it suffices to show mðn, tÞ 1 X lim lim sup aj jþ1 ðiÞ ¼ 0 T!0 n!1 T j¼n
ð65Þ
8t 2 ½0, T,
i ¼ 0, 1, . . . , r, ( ) m P where mðn, tÞ ¼ max m : aj t . j¼n Noticing kþ1 ðiÞ ¼ ðuk ykþiþ1 E½uk ykþiþ1 Þ ðE½uk ykþiþ1 di Þ uk kþiþ1 , we conclude (65) by Lemma 2 and the convergence of the first series in (57).
590
X.-L. Hu and H.-F. Chen
For (62) it suffices to set ¼ 1, F ¼ 1, and 00 k ¼ 0 in CRT stated in [11]. The assertions concerning ci and dj follow from (22), (23) and the assumptions 6¼ 0 and det 6¼ 0. &
By stability of C and by the convergence Ck ! C, there exist N1 > 0 and 2 ð0, 1Þ such that
Theorem 2. Assume that H0-H2 hold, 6¼ 0 and det() 6¼ 0 where and are given in (13) and (20), respectively. Then for k(v) given by (29)(30)
where
k ðvÞ ! fðvÞ
Sð , kÞ ¼ N2
ð66Þ
a:s:
k!1
x0 x0 jj þ Sð , kÞ, jj x^kþ1 xkþ1 jj N1 kþ1 j j^
kþ1 X
ð69Þ
kjþ1
j¼1
ðjjCj Cjj jjxj1 jj þ jjDj Djj jjUj jjÞ
Proof. We rewrite (29) as pffiffiffi 8 k ðvÞ ak ðvÞðk ðvÞ fðvÞÞ þ ak ek ðvÞ, > > pffiffiffi > < if jk ðvÞ ak ðvÞðk ðvÞ fðvÞÞ kþ1 ðvÞ ¼ > þak ek ðvÞj Mvk ðvÞ, ð67Þ > > : 0, otherwise, where ek ðvÞ ¼
pffiffiffi ^k ðk ðvÞ zk Þ: ðvÞðk ðvÞ fðvÞÞ w
Byffiffiffi noticing that f(v) is the root of the linear function p ðx fðvÞÞ, by GCT for the theorem it suffices to show mðki , Ti Þ 1 X ðsÞ aj ej ðvÞ ¼ 0, ð68Þ lim lim sup T!0 i!1 T j¼k i
s ¼ 1, 2, 3, 4
8Ti 2 ½0, T
for any convergent subsequence ki ðvÞ ! ðvÞ, where ð1Þ
ð2Þ
ð3Þ
ð4Þ
ek ðvÞ þ ek ðvÞ þ ek ðvÞ þ ek ðvÞ ¼ ek ðvÞ, ð1Þ
^k Þðk ðvÞ zk Þ, ek ðvÞ ¼ ðwk w pffiffiffi ð2Þ ek ðvÞ ¼ ð ðvÞ wk Þk ðuÞ, pffiffiffi ð3Þ ek ðvÞ ¼ wk yk ðvÞfðvÞ, ð4Þ
ð4Þ
ek ðvÞ ¼ ek ¼ wk k : By the second limit in (31) and the convergence of the second series in (57) it follows that (68) holds for s ¼ 3. Further, by Lemma 3 (the first series in (34)) (68) holds for s ¼ 4. We now show (68) for s ¼ 1. For this we estimate ^k wk j. Since 12 3 > 0, there is a $ > 0 jw such that 12 3 3$ > 0. By (24)(26) it follows that x^kþ1 xkþ1 ¼ Ckþ1 ð^ xk xk Þ þ ðCkþ1 CÞxk þ ðDkþ1 DÞUkþ1 :
and
N2 > 0:
By Theorem 1 jjCk Cjj and jjDk Djj are of order oðk1=2$ Þ, and by Lemma 4, a:s:
xk ! 0. k$ k!1
Uk a:s: ! k$ k!1
0,
Then we have
Sð , kÞ ¼
kþ1 X
kjþ1 o
j¼1
1
j1=22$
X kþ1
1=22$ k 1=22$ j k j¼1
X
1=22$ k 1 k j ¼ o 1=22$
kjþ1 k j¼0 1 0 k
X k 2 X C 1 B ¼ o 1=22$ @ þ A j
k j¼0 k
¼o
1
kjþ1
1=22$
j¼ 2 þ1
k kjþ1
k 1 1 2 þ1 þ 2 ¼ o 1=22$ k2 k 1 k
1 ¼ o 1=22$ : k Noticing that the first term on the right-hand side of (69) exponentially decays, we find that
jj x^k xk jj ¼ o
1
k1=22$
a:s:
Since j v^k vk j and jj x^k xk jj are of the same order, we further have
1 1 ^k wk j ¼ o 3 j v^k vk j ¼ o 3 jj x^k xk jj jw b bk
k 1 ð70Þ ¼ o 1=23 2$ : k
591
Identification for Wiener Systems
By the selection of Mk defined in (14) from (29) it follows that jk ðvÞj k , which together with (70) implies ! 1 1 X X 1 ^k wk Þk ðvÞj ¼ o ak jðw : k3=23 2$ k¼1 k¼1 ð71Þ By Lemma 4 and (70) we have 1 1 X X ^k wk Þyk j ¼ ak jðw
1 yk ^ ð w w Þ k k k$ k1$ k¼1 ! 1 X 1 ¼o < 1: k3=23 3$ k¼1
k¼1
ð72Þ Again by (70) it is seen that 1 X
and s, ki ¼ 1 þ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ
ð77Þ
as i ! 1 and T ! 0. By (76) and the convergence of the second series in (34) it follows that X s s, jþ1 aj wj j ¼ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ j¼k i
as i ! 1 and T ! 0. By the last limit in (31) and (59)(77) it is seen that X s s, jþ1 aj wj yj j¼ki
s X
s, jþ1 aj ðwj jyj j E½wj jyj j þ E½wj jyj jÞ
j¼ki
¼ OðTÞ
^k wk Þk j ak jðw
8s 2 ½ki , . . . , mðki , TÞ. From the recursion (74) it follows that
k¼1
¼
1 X k¼1 1 X k¼1
^k wk Þj jk j ak jðw
o
ð78Þ
ðjk j Ejk jÞ þ Ejk j k3=23 2$
< 1:
ð73Þ
Combining (71)–(73) leads to (68) for s ¼ 1. Finally, we prove (68) for s ¼ 2. For this we first show that there exists an 0 with P0 ¼ 1 such that for any ! 2 0 if ki ðvÞ is a convergent subsequence, then ^s ðs ðvÞ zs Þ, sþ1 ðvÞ ¼ s ðvÞ as w
ð74Þ
sþ1 ðvÞ ¼ s ðvÞ as ws ðs ðvÞ zs Þ ^s Þðs ðvÞ zs Þ þ as ðws w s X ¼ s, ki ki ðvÞ þ s, jþ1 aj wj zj j¼ki s X
þ
^j Þðj ðvÞ zj Þ, s, jþ1 aj ðwj w
j¼ki
which incorporating with (71)–(73) yields sþ1 ðvÞ ¼ ki ðvÞ þ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ: ð79Þ
and jjsþ1 ðvÞ ki ðvÞjj cT, s ¼ ki ,ki þ 1, ...,mðki ,TÞ ð75Þ provided i is large enough and T > 0 is sufficiently small, where c > 0 does not depend on ki . Consider the recursion (74) starting from ki ðvÞ. Let
i, j ¼ ð1 ai wi Þ . . . ð1 aj wj Þ,
i j,
j, jþ1 ¼ 1:
By the first limit in (31) and the convergence of the first series in (58) it is clear that s X
aj wj ¼ OðTÞ 8s 2 ½ki , . . . , mðki , TÞ,
j¼ki
which implies log s, ki ¼ O
s X j¼ki
! aj wj
ð76Þ
This means that the algorithm (29)(30) has no truncation for s 2 ½ki , . . . , mðki , TÞ, when i is large enough and T > 0 is sufficiently small. Since Lemmas 1, 3, 4 and 8 are all proved a.s., the possibly exceptional set is with probability zero. This means that 0 may be taken with P0 ¼ 1 such that for all ! 2 0 (74)(75) are verified for large i and small T > 0. Let ki ðvÞ ! ðvÞ. Then we have mðki , Ti Þ 1 X ð2Þ aj ej ðvÞ T j¼k i mðki , Ti Þ 1 X aj ðj ðvÞ ðvÞÞðwj Ewj Þ T j¼k i mðk i , Ti Þ X 1 aj ðwj Ewj Þ þ TðvÞ j¼k i mðk i , Ti Þ X pffiffiffi 1 aj j ðvÞðEwj Þ: ð80Þ þ T j¼k i
592
X.-L. Hu and H.-F. Chen
By the convergencepofffiffiffi the second series in (58), (75), and limk!1 Ewk ¼ ðvÞ (Lemma 1) we see that the first term on the right-hand side of (80) tends to zero as i ! 1 and T ! 0, while its second term tends to zero as i ! 1 by the convergence of the first pffiffiffi series in (58). Finally, by (75) and limk!1 Ewk ¼ ðvÞ the last term in (80) also tends to zero as i ! 1 and T ! 0. Thus, (68) for s ¼ 2 is proved and meanwhile the proof of Theorem 2 is completed. & Remark 4. Although the estimation algorithms consisting of (15),(16),(26)–(30) do not apply to the feedback control system, they are not complicated and suitable for implementation. Fig. 2.
6.
Numerical Experiments
In this section we give some examples with and without complete satisfaction of conditions required in theorems. For a better comparison in all simulations the linear subsystem is taken to be the same: vkþ1 0:3vk þ 0:6vk1 ¼ uk þ 0:27uk1 0:5uk2 þ 0:79uk3 with noise k 2 N ð0, 0:22 Þ and with iid input fuk g, uk 2 N ð0, 1Þ. Take Mk ¼ 0:1 þ k0:1 , 0 ¼ 2, and ak ¼ 1=k. The inputs for all simulations are taken from the same totality consisting of 6000 data. The difference between examples to be given consists in the nonlinear function fðÞ and in the selection of in bk ¼ 1=k . In all figures shown below the solid lines denote the estimates while the dotted lines the true values. Example 1. Let the nonlinear function be 8 < v 13 v3 , jvj 1, fðvÞ ¼ 2=3, v > 1, : 2=3, v < 1: Take bk ¼ k0:16 , which means ¼ 0.16. Obviously, all conditions required in Theorems 1 and 2 are satisfied. Figures 2 and 3 demonstrate how the estimates given in Section 2 approach to the true values. Figure 2 is for the coefficients of the linear subsystem, and Fig. 3 is for the nonlinear function. The estimate for fðÞ is obtained by separating the interval ½-1:5,1:5 into 100 subintervals with equal length and the values of fðÞ at endpoints of the subintervals are estimated. As shown in the figures, the simulation is consistent with the theoretical analysis.
Fig. 3.
An interesting question is that what will happen if some conditions required in theorems are not met. We show this by examples. Example 2. In (14) it is required that 3 þ < 1=2, which is violated if ¼ 0:6 or ¼ 0:9. We take the same system as that discussed in Example 1 but with different . It is clear the change of does not effect the estimates for coefficients of the linear subsystem. Figures 4 and 5 demonstrate the estimates for the nonlinear function for ¼ 0:6 and ¼ 0:9, respectively. From the figures it is seen that the larger the worse the estimates. As a matter of fact, from (27) we observe that bk ¼ 1=k characterizes the concentrating rate of vk on v. Further, from (28) we see that in order v^k be also concentrated on v, bk ¼ 1=k should not decrease faster than vk v^k , which by Theorem 1 tends to zero not faster than 1=k0:5 . This explains why
should be small.
593
Identification for Wiener Systems
Fig. 4.
Fig. 6.
Fig. 5.
Fig. 7.
Example 3. We now consider the situation where fðÞ is discontinuous, i.e., H1 fails to hold. Let 1, v 0, fðvÞ ¼ 1, v < 0, and let the linear subsystem be the same as that in Example 1. Figure 6 corresponds to the case where the parameters of algorithms completely coincide with those in Example 1, while in Fig. 7 the only change is that ¼ 0:6 rather than 0.16. From these figures we find that the algorithms proposed in the paper work well even H1 does not hold, if is appropriately small.
7.
Concluding Remarks
The paper gives the recursive and strongly consistent estimates for the Wiener systems, when its linear
subsystem is an RTF. From the paper we see that in comparison with the FIR type linear subsystem the infinite impulse response makes a great complication in analysis even under a local Lipschitz-type condition H1 on the nonlinearity. Since Gaussian distribution is even, defined by (13) is zero for any even functions. Thus, Theorem 1 does not apply to identifying Wiener systems with even nonlinearity. Therefore, for further research it is important to solve the even function problem. It is also of interest to weaken the conditions on fðÞ and to seek the conditions guaranteeing nondegeneracy of . To accelerate the algorithms used in the paper is of interest too. It is also worth considering the possible extension to the MIMO case.
References 1. Bendat JS. Nonlinear system analysis and identification from random data. New York Wiley 1999
594 2. Billings SA, Fakhouri FY. Identification of nonlinear systems using the Wiener model. Electron Lett 1978; 13: 502–504 3. Chen HF. Stochastic approximation and its applications. Kluwer, Dordrecht 2002 4. Chen HF. Pathwise convergence of recursive identification algorithms for Hammerstein systems. IEEE Trans Autom Control 2004; 49(10): 1641–1649 5. Chen HF. Recursive identification for Wiener model with nonlinearity being discontinuous piece-wise linear function. IEEE Trans Autom Control 2006; 51(3): 390–400 6. Chow YS, Teicher H. Probability theory. Springer, New York 1978 7. Greblicki W. Nonparametric approach to Wiener system identification. IEEE Trans Circuits Syst-I: Fundam Theory Appl 1997; 44(6): 538–545 8. Greblicki W. Recursive identification of Wiener system. Int J Appl Math Comp Sci 2001; 11(4): 977–991 9. Hasiewicz Z. Identification of a linear system observed though zero-memory nonlinearity. Int J Syst Sci 1987; 18: 1595–1607 10. Hu XL, Chen HF. Strong consistency of recursive identification for Wiener systems. Automatica 2005; 41(11): 1905–1916 11. Hu XL, Chen HF. Recursive identification for Wiener systems using Gaussian signals. to appear in Asian J Control 2007.
X.-L. Hu and H.-F. Chen
12. Hunter IW, Korenberg MJ. The identification of nonlinear biological systems: Wiener and Hammerstein cascade models. Biol Cybern 1986; 55: 135–144 13. Krzyzak A, Sasiadek JZ, Kegl B. Non-parametric identification of dynamic non-linear systems by a Hermite Series Approach. Int J Syst Sci 2001; 32(10): 1261–1285 14. Pajunen GA. Adaptive control of Wiener type nonlinear systems. Automatica 1992; 28: 781–785 15. Nordsjo¨ AE, Zetterberg LH. Identification of certain time-varying nonlinear Wiener and Hammerstein systems. IEEE Trans Signal Proc. 2001; 49: 577–592 16. Verhaegen M, Westwick D. Identifying MIMO Wiener systems in the context of subspace model identification methods. Int J Control 1996; 63(2): 331–349 17. Vo¨ro¨s J. Parameter identification of Wiener systems with discontinuous nonlinearities. Syst Control Lett 2001; 44: 363–372 18. Westwick DT, Kearney RE. A new algorithm for the identification of multiple input Wiener systems. Biol Cybern 1992; 68: 75–85 19. Wigren T. Recursive prediction error identification using the nonlinear Wiener mode. Automatica 1993; 29: 1011–1025 20. Wigren T. Convergence analysis of recursive algorithms based on the nonlinear Wiener model. IEEE Trans Automat Control 1994; 39: 2191–2206