Auxiliary model-based least-squares identification ... - Semantic Scholar
Systems & Control Letters 56 (2007) 373 – 380 www.elsevier.com/locate/sysconle
Auxiliary model-based least-squares identification methods for Hammerstein output-error systems夡 Feng Ding a,∗,1 , Yang Shi b , Tongwen Chen c a Control Science and Engineering Research Center, Southern Yangtze University, Wuxi 214122, P.R. China b Department of Mechanical Engineering, University of Saskatchewan, Saskatoon, Canada S7N 5A9 c Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Canada T6G 2V4
Received 28 August 2004; received in revised form 26 October 2006; accepted 30 October 2006 Available online 2 January 2007
Abstract The difficulty in identification of a Hammerstein (a linear dynamical block following a memoryless nonlinear block) nonlinear output-error model is that the information vector in the identification model contains unknown variables—the noise-free (true) outputs of the system. In this paper, an auxiliary model-based least-squares identification algorithm is developed. The basic idea is to replace the unknown variables by the output of an auxiliary model. Convergence analysis of the algorithm indicates that the parameter estimation error consistently converges to zero under a generalized persistent excitation condition. The simulation results show the effectiveness of the proposed algorithms. © 2006 Elsevier B.V. All rights reserved. Keywords: Recursive identification; Parameter estimation; Least squares; Multi-innovation identification; Hierarchical identification; Auxiliary model; Convergence properties; Stochastic gradient; Hammerstein models; Wiener models; Martingale convergence theorem
1. Introduction Two typical classes of nonlinear systems—linear timeinvariant blocks following (or followed by) static nonlinear blocks—are Hammerstein and Wiener (H–W) nonlinear systems, which are common in industry, e.g., the valve saturation nonlinearities, dead-zone nonlinearities and linear systems equipped with nonlinear sensors [35]. In general, existing identification approaches for H–W models can be roughly divided into two categories: the iterative and the recursive algorithms. In order to distinguish on-line from off-line calculation, we use iterative for off-line algorithms, and recursive for on-line ones. We imply that a recursive algorithm can be on-line implemented, but an iterative one cannot. For a 夡 This research was supported by the Natural Sciences and Engineering Research Council of Canada and the National Natural Science Foundation of China (No. 60574051 and 60528007). ∗ Corresponding author. Tel.: 86510 88637783; fax: +86510 85910652. E-mail addresses: [email protected], [email protected] (F. Ding), [email protected] (T. Chen). 1 F. Ding is also an Adjunct Professor in the College of Automation at the Nanjing University of Technology, Nanjing 210009, P.R. China.
0167-6911/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.sysconle.2006.10.026
recursive algorithm, new information (input and/or output data) is always used in the algorithm which recursively computes the parameter estimates every step as time increases. Some iterative and/or off-line algorithms of H–W models were discussed in [2,1,4,6,7,22,23,25,28,30,36,11,20] and other recursive and/or on-line algorithms were studied in, e.g., [35,11,20,32,33,3,5,34]. In the identification area of nonlinear systems, Bai reported a two-stage identification algorithm for Hammerstein–Wiener nonlinear systems based on singular value decomposition (SVD) and least squares [1] and studied identification problem of systems with hard input nonlinearities of known structure [2]; Vörös presented a half-substitution algorithm to identify Hammerstein systems with two-segment nonlinearities and with multisegment piecewise-linear characteristics, but no convergence analysis was carried out [32,33]; Cerone and Regruto analyzed parameter error bounds in the Hammerstein models by using the output measurement error bound [6]. Also, Pawlak used the series expansion approach to study the identification of Hammerstein nonlinear outputerror (state-space) models [29]. Recently, an iterative leastsquares and a recursive least-squares identification methods were reported in [11], and an iterative gradient and a recursive
374
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
as follows [6,11,20]:
v(t) u(t)
f (·)
u¯ (t)
B(z) A(z)
x(t)
+
y(t)
u(t) ¯ = f (u(t)) = c1 1 (u(t)) + c2 2 (u(t)) m + · · · + cm m (u(t)) = cj j (u(t)).
(1)
j =1
Fig. 1. The Hammerstein nonlinear output-error system.
Then the Hammerstein nonlinear output-error model in Fig. 1 may be expressed as stochastic gradient algorithms were developed in [20] for nonlinear ARMAX models based only on the available input–output data, and the convergence properties of the recursive algorithms involved were proved. However, most of the identification approaches in the literature assume that the systems under consideration are nonlinear ARX models, or equation-error-like models [35,2,1,7,3]. That is, each element of the information vector consisting of input–output data is measured. In this paper, we focus on the identification problem of a class of Hammerstein output-errortype nonlinear systems and present an auxiliary-model leastsquares (AMLS) algorithm, which is different from the ones mentioned above in that the information vector in our identification model contains unknown variables (namely, unavailable noise-free outputs), and we adopt an auxiliary model or reference model to estimate these unknown variables and further use the outputs of the auxiliary model instead of the unknown noisefree outputs to identify the system parameters. The basic idea is to extend the Landau’s output-error method to study identification problem of nonlinear systems [21,31]. To the best of our knowledge, few publications addressed identification methods of Hammerstein nonlinear output-error systems, especially the convergence problem of the algorithms involved, which are the focus of this work. The objective of this paper is, by means of the auxiliary model identification principle, to derive an algorithm to estimate the system parameters of the nonlinear output-error models based on the available input–output data {u(t), y(t)}, and to study the properties of the algorithm involved. Briefly, the paper is organized as follows. Section 2 describes the identification algorithms related to the Hammerstein systems. Section 3 analyzes the properties of the proposed stochastic algorithm. Section 4 provides an illustrative example to show the effectiveness of the algorithm proposed. Finally, we offer some concluding remarks in Section 5. 2. The algorithm description Consider the Hammerstein output-error system shown in Fig. 1 that consists of a nonlinear memoryless element followed by a linear output-error model [6,29], where the true output (namely, the noise-free output) x(t) and the inner variable u(t) ¯ (namely, the output of the nonlinear block) are unmeasurable, u(t) is the system input, y(t) is the measurement of x(t), v(t) is an additive noise with zero mean. The nonlinear part in the Hammerstein model is a polynomial of a known order in the input [7,22,28], or, more generally, a nonlinear function of a known basis (1 , 2 , . . . , m )
x(t) =
B(z) B(z) u(t) ¯ = [c1 1 (u(t)) + c2 2 (u(t)) A(z) A(z) + · · · + cm m (u(t))],
y(t) = x(t) + v(t). Here, A(z) and B(z) are polynomials in the shift operator z−1 [z−1 y(t) = y(t − 1)] with A(z) = 1 + a1 z−1 + a2 z−2 + · · · + an z−n , B(z) = b1 z−1 + b2 z−2 + b3 z−3 + · · · + bn z−n . Notice that for the Hammerstein model shown in Fig. 1, f (u) and G(z) := B(z)/A(z) are actually not unique. Any pair (f (u), G(z)/) for some nonzero and finite constant would produce identical input and output measurements. In other words, any identification scheme cannot distinguish between (f (u), G(z)) and (f (u), G(z)/). Therefore, to get a unique parameterization, without loss of generality, one of the gains of f (u) and G(z) has to be fixed. There are several ways to normalize the gains [1,6,22]. Here, we adopt the assumption [22,3]: the first coefficient of the function f (·) equals 1; i.e., c1 = 1 [11,20]. Eq. (2) can be rewritten as a recursive form x(t) = − = −
n i=1 n i=1
ai x(t − i) + ai x(t − i) +
n i=1 n
bi u(t ¯ − i) bi
i=1
m
cj j (u(t − i)).
j =1
Define the parameter vector and information vector 0 (t) as ⎡ ⎤ a ⎢ c1 b ⎥ ⎢ ⎥ n0 c b⎥ =⎢ ⎢ 2. ⎥ ∈ R , ⎣ .. ⎦ cm b ⎡ ⎤ −x(t − 1) ⎢ −x(t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 , n0 := (m + 1)n, 0 (t) = ⎢ (2) . ⎢ ⎥ ⎣ −x(t − n) ⎦ (t) a = [a1 , a2 , . . . , an ]T ∈ Rn ,
b = [b1 , b2 , . . . , bn ]T ∈ Rn ,
c = [c2 , c3 , . . . , cm ]T ∈ Rm−1 , (t) = [T1 (t), T2 (t), . . . , Tm (t)]T ∈ Rmn ,
(3)
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
j (t) = [j (u(t − 1)), j (u(t − 2)), n
. . . , j (u(t − n))] ∈ R , T
j = 1, 2, . . . , m.
(4)
Here, the superscript T denotes the matrix transpose. Then we obtain the identification model of the HOE systems: x(t) = T0 (t),
y(t) = T0 (t) + v(t).
(6)
v(t − p + 1) From (5) and (6), we have Y(p, t) = 0 (p, t) + V(p, t).
(7)
Because 0 (p, t) contains unmeasured inner variables x(t − i), i = 1, 2, . . . , n, the standard least-squares method cannot be applied directly to obtain the least-squares estimate [T0 (p, t)0 (p, t)]−1 T0 (p, t)Y(p, t) of the parameter vector by minimizing the cost function Y(p, t) − 0 (p, t)2 even if V(p, t) is a “white” noise vector with zero mean, where X2 = tr[XXT ]. If these unknowns x(t−i) are replaced by the outputs xa (t−i) of an auxiliary (or reference) model, xa (t) = Pa (z)f (u(t)),
According to the auxiliary model identification principle: the unknown variables x(t − i) in 0 (t) are replaced by the output xa (t − i) of the auxiliary model, 0 (t) by (t) and 0 (p, t) by (p, t); then it is easy from (7) to get a recursive pseudo-linear regression least-squares identification algorithm based on the auxiliary model as follows:
(5)
The systems in (5) by parameterization in (2) contain more parameters than actually needed, but the benefit to do this is that a linear regression form can be obtained. However, the disadvantage is that it requires extra computation, which can be approached using the hierarchical estimation methods based on the hierarchical identification principle [8,12–15,18]. From (5), we can see that the information vector 0 (t) consists of unknown noise-free outputs x(t − i) and system inputs. Let p be the data length (p?n0 ), and define ⎡ ⎤ y(t) ⎢ y(t − 1) ⎥ ⎥, Y(p, t) := ⎢ .. ⎣ ⎦ . y(t − p + 1) ⎤ ⎡ T0 (t) ⎢ T (t − 1) ⎥ 0 ⎥ ⎢ 0 (p, t) := ⎢ ⎥, .. ⎦ ⎣ . T (t − p + 1) ⎡ 0 ⎤ v(t) ⎢ v(t − 1) ⎥ ⎥. V(p, t) := ⎢ .. ⎣ ⎦ .
375
or xa (t) = T (t)a ,
then the identification problem of can be solved using xa (t) instead of x(t). Here, Pa (z), (t) and a are the transfer function, information vector and parameter vector of the auxiliary model, respectively. If we use the estimate of B(z)/A(z) as an ˆ auxiliary model Pa (z), namely, take a to be the estimate (t) of , and (t) to be the regressive vector of xa (t) and u(t), and use (t) as 0 (t), then the identification algorithms based on this idea are called the auxiliary (or reference) model identification method or output-error method [21,31]. Of course, there are other ways to choose auxiliary models, e.g., using the finite impulse response model [9,10].
ˆ = (t ˆ − 1) + P(t)T (p, t)[Y(p, t) − (p, t)(t ˆ − 1)], (t) (8) P−1 (t) = P−1 (t − 1) + T (p, t)(p, t), ˆ xa (t − i) = T (t − i)(t), ⎡ ⎢ (p, t) = ⎢ ⎣
T (t) T (t − 1) .. .
P(0) = p0 I,
i = 0, 1, 2, . . . , n,
(9) (10)
⎤ ⎥ ⎥, ⎦
T (t − p + 1) ⎤ −xa (t − 1) ⎢ −xa (t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 . (t) = ⎢ . ⎢ ⎥ ⎣ −x (t − n) ⎦ a (t) ⎡
(11)
Here, p0 is a large positive number, e.g., p0 =106 , and we refer to xa (t) as the estimate of x(t). Eqs. (8)) to (11) is also called the AMLS identification algorithm of estimating in (7) of the Hammerstein output-error model, HOE-AMLS algorithm for short. Note that the algorithm is combined in the sense the parameters and true outputs are estimated simultaneously, and can be implemented on-line. Note c1 = 1, the estimates aˆ = [aˆ 1 , aˆ 2 , . . . , aˆ n ]T and bˆ = [bˆ1 , bˆ2 , . . . , bˆn ]T of a and b can be read from the first and ˆ respectively. Let ˆ i be the ith element second n entries of , ˆ referring to the definition of , then the estimates of cj , of , j = 2, 3, . . . , m, may be computed by cˆj = ˆ j n+i /bˆi ,
j = 2, 3, . . . , m; i = 1, 2, . . . , n.
Since we do not need such n estimates cˆj , we may calculate the c-parameters by a least-squares fit or take their average like in [11,20] as the estimate of cj , i.e., cˆj =
n 1 ˆ j n+i , j = 2, 3, . . . , m. n bˆi i=1
3. The performance analysis Let us introduce some notation first. The symbol I stands for an identity matrix of appropriate sizes; |X| = det[X] represents the determinant of the square matrix X; max [X] and min [X] represent the maximum and minimum eigenvalues of the symmetric matrix X, respectively; for g(t) 0, we write f (t) = O(g(t)) if there exists positive constants 1 and t0 such that |f (t)| 1 g(t) for t t0 .
376
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
Dividing (13) by [ln |P−1 (t)] and summing for t give
Define P−1 0 (t) =
t i=1 −1
r(t) = tr[P
T0 (p, i)0 (p, i) + (t)];
∞ T (t − i)P(t)(t − i)
1 I; p0
Hence, we easily get (12)
Lemma 1. For each i (i = 0, 1, . . . , p − 1), the following inequalities hold:
2
(j − i)P(j )(j − i) ln |P T (t − i)P(t)(t − i)
t=1
[ln |P−1 (t)|]
< ∞,
−1
(t)| + n0 ln p0 , a.s.,
a.s., for any > 1.
P−1 (t − 1) = P−1 (t) − T (p, t)(p, t) P−1 (t) − (t − i)T (t − i) = P−1 (t)[I − P(t)(t − i)T (t − i)]. Taking determinants on both sides and using the formula det[I+ DE] = det[I + ED] yields |P−1 (t − 1)| |P−1 (t)||I − P(t)(t − i)T (t − i)|
˜ = (t) ˆ − , (t)
(14) (15)
Lemma 2. For the system in (7) and the HOE-AMLS algorithm in (8)–(11), assume the noise sequence {v(t)} with zero mean and bounded time-varying variance satisfies [24]:
where {v(t), Ft } is a martingale sequence defined on a probability space {, F, P } and {Ft } is the algebra sequence generated by {v(t)}. Then the following inequality holds: E[W (t) + S(t)|Ft−1 ] W (t − 1) + S(t − 1) p−1
T (t − i)P(t)(t − i)¯ 2v ,
a.s.,
i=0
Hence
where
T (t − i)P(t)(t − i)
|P−1 (t)| − |P−1 (t − 1)| |P−1 (t)|
.
(13)
Replacing t with j and summing for j from 1 to t yields (noting that |P−1 (t)| is a nondecreasing function of t)
− 1)|
|P−1 (j )|
j =1 t |P−1 (j )|
−1 j =1 |P (j −1)|
a.s.,
(16)
i=1
˜ = Y(t)
1 2
˜ (p, t)(t) (17) (18)
Here, (A3) guarantees that S(t) 0. The proof can be done in a similar way in [11,19] and is omitted here. Next, we give prove the main results of this paper.
dx |P−1 (j )|
|P−1 (t)| dx = ln |P−1 (t)| − ln |P−1 (0)| −1 |P (0)| x = ln |P−1 (t)| + n0 ln p0 ,
˜ T (i)Y(i), ˜ U
˜ ˜ U(t) = −(p, t)(t).
j =1
|P−1 (j )| − |P−1 (j
S(t) = 2
t
ˆ − V(p, t)], + [Y(p, t) − (p, t)(t)
T (j − i)P(j )(j − i) t
˜ W (t) = ˜ (t)P−1 (t)(t).
+ 2p
= |P−1 (t)|[1 − T (t − i)P(t)(t − i)].
=
˜ and a nonDefine the parameter estimation error vector (t) negative definite function W (t) as
(A1) E[v(t)|Ft−1 ] = 0, a.s., (A2) E[v 2 (t)|Ft−1 ] = 2v (t) ¯ 2v < ∞, a.s., 1 − 21 is strictly positive real, (A3) H (z) = A(z)
Proof. From the definition of P(t) in (9), we have
|P−1 (∞)|
−1 1
= −1 −1 −1 − 1 x(ln x) (ln x) |P (0)| |P (0)| 1 1 1 − = − 1 [ln |P−1 (0)|]−1 [ln |P−1 (∞)|]−1 < ∞, a.s. dx
T
T
j =1 ∞
t
|P−1 (t)|[ln |P−1 (t)|]
|P−1 (∞)|
The stochastic process theory [27,16] and the stochastic martingale theory [11,10,17,24] are the main tool of analyzing performance of identification algorithms. Here, we prove the main convergence results of the HOE-AMLS algorithm by using the martingale convergence theorem. To do this, some mathematical preliminaries are required.
1
∞ |P−1 (t)| − |P−1 (t − 1)| t=1
|P−1 (t)|r n0 (t); r(t) n0 max [P−1 (t)]; ln |P−1 (t)| = O(ln r(t)).
t
[ln |P−1 (t)|]
t=1
r0 (t) = tr[P−1 0 (t)].
a.s.
Theorem 1. For the system in (7), assume that (A1)–(A3) hold, and A(z) is stable, i.e., all zeros of A(z) are inside the unit circle. Then for any > 1, the parameter estimation error by
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
the HOE-AMLS algorithm in (8)–(11) satisfies:
[ln r0 (t)] 2 ˆ , a.s. (t) − = O min [P−1 0 (t)]
(22), we have
min [P−1 0 (t)] + [ln r0 (t)] [ln r0 (t)] =O , a.s. for any > 1. min [P−1 0 (t)]
(19)
Let [ln |P−1 (t)|]
.
Since ln |P−1 (t)|] is nondecreasing, according to Lemma 2, we have E[Z(t)|Ft−1 ] Z(t − 1) + 2p
p−1
T (t − i)P(t)(t − i)
i=0
[ln |P−1 (t)|]
¯ 2v ,
a.s. (20)
Using Lemma 1 and applying the martingale convergence theorem (Lemma D.5.3 in [24]) to (20), we conclude that Z(t) converges a.s. to a finite random variable, say, Z0 ; i.e., Z(t) =
W (t) + S(t) [ln |P−1 (t)|]
→ Z0 < ∞,
a.s.,
W (t) = O([ln |P−1 (t)|] ), a.s., S(t) = O([ln |P−1 (t)|] ), a.s.
(21)
Since H (z) is a strictly positive real function, from the definition of S(t) and referring to [11], we have 2 ˜ U(i) = O([ln |P−1 (t)|] ).
i=1
(i) − V(p, i)2 k1
t
ˆ = (t ˆ − 1) + P(t)(t)[y(t) − T (t)(t ˆ − 1)], (t) P(t − 1)(t)T (t)P(t − 1) , 1 + T (t)P(t − 1)(t)
ˆ xa (t) = T (t)(t), ⎡ ⎤ −xa (t − 1) ⎢ −xa (t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 . (t) = ⎢ . ⎢ ⎥ ⎣ −x (t − n) ⎦ a (t)
(24) (25) (26)
(27)
B(z) u(t) ¯ + v(t), A(z) A(z) = 1 + a1 z−1 + a2 z−2 = 1 − 1.60z−1 + 0.80z−2 , y(t) =
(22)
2 ˜ U(i) + k2
i=1
= O([ln |P−1 (t)|] ) = O([ln r(t)] ).
From Theorem 1, we can see that if [ln r0 (t)] = o(min [P−1 0 (t)]) (namely, the generalized persistent excitation (PE) condition [10], the weak PE condition [10] or the strong PE condition [16] holds), then the estimation error (t) − converges to zero. This requires that G(z) is stable, u(t) ¯ is a PE signal, and so is u(t) = f −1 (u(t)). ¯ The positive real condition in (A3) depends on the unknown model parameters and is a more restrictive assumption than A(z) being stable [27]. ˆ − 1) ∈ R1 is called the As p = 1, e(t) := y(t) − T (t)(t ˆ − 1)] ∈ innovation [26], and E(t) := [Y(p, t) − (p, t)(t Rp is referred as to the innovation vector. p here may be also known as the innovation length, the algorithm in (8)–(11) is also called the auxiliary model multi-innovation least-squares identification algorithm for Hammerstein output-error systems. When p = 1, we get a simple HOE-AMLS algorithm:
An example is given to demonstrate the effectiveness of the proposed algorithms. Consider the following system:
Like in [11,19], it follows that there exist positive constants k1 and k2 such that t
4. Example
i=1
From (19), (21) and (12), we have
−1 (t)|] [ln |P 2 ˜ (t) =O min [P−1 (t)]
[ln r(t)] =O , a.s. for any > 1. min [P−1 (t)]
This proves Theorem 1.
P(t) = P(t − 1) −
or
t
T ˜ ˜ (t)P−1 (t)(t) W (t) 2 ˜ (t) = . −1 min [P (t)] min [P−1 (t)]
W (t) + S(t)
[ln r0 (t)]
ˆ − 2 = O (t)
Proof. From the definition of W (t), we have
Z(t) =
377
(23)
Using the similar way in [11,19], it is easy to get r(t)=O(r0 (t)), min [P−1 (t)] = O(min [P−1 0 (t)] + O([ln r0 (t)] ). Thus, from
B(z) = b1 z−1 + b2 z−2 = 0.85z−1 + 0.65z−2 , u(t) ¯ = f (u(t)) = c1 u(t) + c2 u2 (t) + c3 u3 (t) = u(t) + 0.5u2 (t) + 0.25u3 (t), s = [a1 , a2 , b1 , b2 , c2 , c3 ]T . {u(t)} is taken as a PE signal sequence with zero mean and unit variance 2u = 1.002 , and {v(t)} as a white noise sequence with zero mean and constant variance 2v = 0.502 and 2v = 2.002 . Apply the HOE-AMLS algorithm with p = 1 to estimate the parameters of this system, the parameter estimates = [ 1 , 2 , . . . , 8 ]T and s and their errors with different
378
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
Table 1 The estimates of (2v = 0.502 ) t
1
2
3
4
5
6
7
8
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.56901 −1.57996 −1.57992 −1.59002 −1.58361 −1.59141 −1.59375 −1.59131 −1.59376
0.77630 0.78122 0.78325 0.78999 0.78623 0.79217 0.79512 0.79218 0.79477
0.68233 0.79989 0.83494 0.87510 0.85321 0.83909 0.86519 0.85338 0.85300
1.04526 0.91609 0.85387 0.79772 0.71349 0.69136 0.69088 0.68343 0.66855
0.37656 0.35609 0.36760 0.37839 0.37916 0.37660 0.38268 0.38296 0.38739
0.40992 0.37159 0.39116 0.39249 0.39854 0.38461 0.38439 0.38079 0.37451
0.21864 0.22293 0.20882 0.19782 0.20294 0.21061 0.20713 0.21037 0.21052
0.13271 0.15234 0.15149 0.15553 0.17042 0.17118 0.16963 0.16826 0.16821
20.43326 13.14011 10.34029 7.93938 5.08542 4.10887 3.95682 3.62598 3.03148
True values
−1.60000
0.80000
0.85000
0.65000
0.42500
0.32500
0.21250
0.16250
Table 2 The estimates of (ai , bi , ci ) (2v = 0.502 ) t
a1
a2
b1
b2
c2
c3
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.56901 −1.57996 −1.57992 −1.59002 −1.58361 −1.59141 −1.59375 −1.59131 −1.59376
0.77630 0.78122 0.78325 0.78999 0.78623 0.79217 0.79512 0.79218 0.79477
0.68233 0.79989 0.83494 0.87510 0.85321 0.83909 0.86519 0.85338 0.85300
1.04526 0.91609 0.85387 0.79772 0.71349 0.69136 0.69088 0.68343 0.66855
0.47202 0.42540 0.44918 0.46220 0.50149 0.50256 0.49934 0.50296 0.50717
0.22370 0.22250 0.21376 0.21051 0.23835 0.24930 0.24247 0.24636 0.24920
20.05581 13.13781 9.97809 7.41968 3.15539 2.05765 2.08352 1.66260 1.00617
True values
−1.60000
0.80000
0.85000
0.65000
0.50000
0.25000
Table 3 The estimates of (2v = 2.002 ) t
1
2
3
4
5
6
7
8
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.54443 −1.55876 −1.55951 −1.58706 −1.56785 −1.58573 −1.59032 −1.58242 −1.59036
0.75677 0.76487 0.76570 0.78914 0.77375 0.78720 0.79314 0.78315 0.79298
0.35049 0.67966 0.77079 0.85335 0.81059 0.77961 0.87038 0.83949 0.84170
1.37101 1.14544 1.03161 0.93939 0.78027 0.76794 0.77592 0.76061 0.71035
0.41764 0.33138 0.37081 0.35883 0.33976 0.32739 0.34143 0.33394 0.34846
0.43136 0.36260 0.42585 0.43890 0.47286 0.43761 0.44251 0.43979 0.42285
0.25470 0.26717 0.21899 0.18798 0.20004 0.22060 0.20796 0.21446 0.21475
0.10459 0.13651 0.12295 0.13346 0.17329 0.16684 0.16391 0.15866 0.16187
41.00730 24.88620 18.98097 14.78786 10.27459 9.39654 8.89719 8.55583 6.40464
True values
−1.60000
0.80000
0.85000
0.65000
0.42500
0.32500
0.21250
0.16250
Table 4 The estimates of (ai , bi , ci ) (2v = 2.002 ) t
a1
a2
b1
b2
c2
c3
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.54443 −1.55876 −1.55951 −1.58706 −1.56785 −1.58573 −1.59032 −1.58242 −1.59036
0.75677 0.76487 0.76570 0.78914 0.77375 0.78720 0.79314 0.78315 0.79298
0.35049 0.67966 0.77079 0.85335 0.81059 0.77961 0.87038 0.83949 0.84170
1.37101 1.14544 1.03161 0.93939 0.78027 0.76794 0.77592 0.76061 0.71035
0.75310 0.40206 0.44694 0.44386 0.51259 0.49490 0.48129 0.48800 0.50463
0.40149 0.25613 0.20165 0.18118 0.23443 0.25011 0.22509 0.23203 0.24150
43.00408 24.82545 18.52710 14.04942 6.65812 6.43023 6.10921 5.36474 2.91094
True values
−1.60000
0.80000
0.85000
0.65000
0.50000
0.25000
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
379
1
0.8
0.6 2
v = 0.502
0.4
0.2
2
v = 2.002
0 0
500
1000
1500
2000
2500
3000
2500
3000
t Fig. 2. The parameter estimation errors vs. t.
1
0.8
s
0.6 2
= 0.502 v
0.4
0.2
0
2v = 2.002 0
500
1000
1500
2000
t Fig. 3. The parameter estimation errors s vs. t.
noise variances are shown in Tables 1–4, and the parameter estimation errors and s vs. t are shown in Figs. 2 and 3, ˆ =(t)−/, s =ˆ s (t)−s /, ˆ s (t), is the estimate of s . From Tables 1–4 and Figs. 2 and 3, we can draw the following conclusions:
Hammerstein output-error models. The analysis using the martingale convergence theorem indicates that the proposed algorithm can give consistent parameter estimation.
• A high noise level results in a slow rate of convergence of the parameter estimates to the true parameters. • It is clear that the errors and s are becoming smaller (in general) as t increases. This confirms the proposed theorem.
[1] E.W. Bai, An optimal two-stage identification algorithm for Hammerstein-Wiener nonlinear systems, Automatica 34 (3) (1998) 333–338. [2] E.W. Bai, Identification of linear systems with hard input nonlinearities of known structure, Automatica 38 (5) (2002) 853–860. [3] E.W. Bai, A blind approach to the Hammerstein–Wiener model identification, Automatica 38 (6) (2002) 967–979. [4] E.W. Bai, A random least-trimmed-squares identification algorithm, Automatica 39 (9) (2003) 1561–1659. [5] M. Boutayeb, M. Darouach, Recursive identification method for MISO Wiener–Hammerstein model, IEEE Trans. Automat. Control 40 (2) (1995) 287–291.
5. Conclusions A recursive AMLS algorithm based on replacing unavailable variables (noise-free outputs) by their estimates is derived for
References
380
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
[6] V. Cerone, D. Regruto, Parameter bounds for discrete-time Hammerstein models with bounded output errors, IEEE Trans. Automat. Control 48 (10) (2003) 1855–1860. [7] F. Chang, R. Luus, A noniterative method for identification using Hammerstein model, IEEE Trans. Automat. Control 16 (5) (1971) 464–468. [8] F. Ding, T. Chen, Hierarchical gradient-based identification of multivariable discrete-time systems, Automatica 41 (2) (2004) 315–325. [9] F. Ding, T. Chen, Identification of dual-rate systems based on finite impulse response models, Internat. J. Adapt. Control Signal Process. 18 (7) (2004) 589–598. [10] F. Ding, T. Chen, Combined parameter and output estimation of dualrate systems using an auxiliary model, Automatica 40 (10) (2004) 1739–1748. [11] F. Ding, T. Chen, Identification of Hammerstein nonlinear ARMAX systems, Automatica 41 (9) (2005) 1479–1489. [12] F. Ding, T. Chen, Hierarchical least squares identification methods for multivariable systems, IEEE Trans. Automat. Control 50 (3) (2005) 397–402. [13] F. Ding, T. Chen, Hierarchical identification of lifted state-space models for general dual-rate systems, IEEE Trans. Circuits Systems I. Regular Papers 52 (6) (2005) 1179–1187. [14] F. Ding, T. Chen, Iterative least squares solutions of coupled Sylvester matrix equations, Systems Control Lett. 54 (2) (2005) 95–107. [15] F. Ding, T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Trans. Automat. Control 50 (8) (2005) 1216–1221. [16] F. Ding, T. Chen, Performance bounds of forgetting factor least squares algorithm for time-varying systems with finite measurement data, IEEE Trans. Circuits Systems I. Regular Papers 52 (3) (2005) 555–566. [17] F. Ding, T. Chen, Parameter estimation of dual-rate stochastic systems by using an output error method, IEEE Trans. Automat. Control 50 (9) (2005) 1436–1441. [18] F. Ding, T. Chen, On iterative solutions of general coupled matrix equations, SIAM J. Control Optim. 44 (6) (2006) 2269–2284. [19] F. Ding, H.B. Chen, M. Li, Multi-innovation least squares identification methods based on the auxiliary model for MISO systems, Appl. Math. Comput., 2007, in press. [20] F. Ding, Y. Shi, T. Chen, Gradient-based identification algorithms for Hammerstein nonlinear ARMAX models, Nonlinear Dynam. 45 (1–2) (2006) 31–43.
[21] L. Dugard, I.D. Landau, Recursive output error identification algorithms theory and evaluation, Automatica 16 (5) (1980) 443–462. [22] P.G. Gallman, A comparison of two Hammerstein model identification algorithms, IEEE Trans. Automat. Control 21 (1) (1976) 124–126. [23] F. Giri, F.Z. Chaoui, Y. Rochdi, Parameter identification of a class of Hammerstein plants, Automatica 37 (5) (2001) 749–756. [24] G.C. Goodwin, K.S. Sin, Adaptive Filtering, Prediction and Control, Prentice-Hall, Englewood Cliffs, NJ, 1984. [25] N.D. Haist, F. Chang, R. Luus, Nonlinear identification in the presence of correlated noise using a Hammerstein model, IEEE Trans. Automat. Control 18 (5) (1973) 553–555. [26] L. Ljung, System Identification: Theory for the User, second ed., Prentice-Hall, Englewood Cliffs, NJ, 1999. [27] L. Ljung, T. Söderström, Theory and Practice of Recursive Identification, MIT Press, Cambridge, MA, 1983. [28] K.S. Narendra, P.G. Gallman, An iterative method for the identification of nonlinear systems using a Hammerstein model, IEEE Trans. Automat. Control 11 (3) (1966) 546–550. [29] M. Pawlak, On the series expansion approach to the identification of Hammerstein system, IEEE Trans. Automat. Control 36 (6) (1991) 763–767. [30] P. Stoica, On the convergence of an iterative algorithm used for Hammerstein system identification, IEEE Trans. Automat. Control 26 (4) (1981) 967–969. [31] P. Stoica, T. Söderström, Analysis of an output error identification algorithm, Automatica 17 (6) (1981) 861–863. [32] J. Vörös, Iterative algorithm for parameter identification of Hammerstein systems with two-segment nonlinearities, IEEE Trans. Automat. Control 44 (11) (1999) 2145–2149. [33] J. Vörös, Modeling and parameter identification of systems with multisegment piecewise-linear characteristics, IEEE Trans. Automat. Control 47 (1) (2002) 184–189. [34] T. Wigren, Convergence analysis of recursive identification algorithms based on the nonlinear Wiener model, IEEE Trans. Automat. Control 39 (11) (1994) 2191–2206. [35] T. Wigren, A.E. Nordsjö, Compensation of the RLS algorithm for output nonlinearities, IEEE Trans. Automat. Control 44 (10) (1999) 1913–1918. [36] Y. Zhu, Estimation of an N-L-N Hammerstein–Wiener model, Automatica 38 (9) (2002) 1607–1614.
Auxiliary model-based least-squares identification methods for Hammerstein output-error systems夡 Feng Ding a,∗,1 , Yang Shi b , Tongwen Chen c a Control Science and Engineering Research Center, Southern Yangtze University, Wuxi 214122, P.R. China b Department of Mechanical Engineering, University of Saskatchewan, Saskatoon, Canada S7N 5A9 c Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Canada T6G 2V4
Received 28 August 2004; received in revised form 26 October 2006; accepted 30 October 2006 Available online 2 January 2007
Abstract The difficulty in identification of a Hammerstein (a linear dynamical block following a memoryless nonlinear block) nonlinear output-error model is that the information vector in the identification model contains unknown variables—the noise-free (true) outputs of the system. In this paper, an auxiliary model-based least-squares identification algorithm is developed. The basic idea is to replace the unknown variables by the output of an auxiliary model. Convergence analysis of the algorithm indicates that the parameter estimation error consistently converges to zero under a generalized persistent excitation condition. The simulation results show the effectiveness of the proposed algorithms. © 2006 Elsevier B.V. All rights reserved. Keywords: Recursive identification; Parameter estimation; Least squares; Multi-innovation identification; Hierarchical identification; Auxiliary model; Convergence properties; Stochastic gradient; Hammerstein models; Wiener models; Martingale convergence theorem
1. Introduction Two typical classes of nonlinear systems—linear timeinvariant blocks following (or followed by) static nonlinear blocks—are Hammerstein and Wiener (H–W) nonlinear systems, which are common in industry, e.g., the valve saturation nonlinearities, dead-zone nonlinearities and linear systems equipped with nonlinear sensors [35]. In general, existing identification approaches for H–W models can be roughly divided into two categories: the iterative and the recursive algorithms. In order to distinguish on-line from off-line calculation, we use iterative for off-line algorithms, and recursive for on-line ones. We imply that a recursive algorithm can be on-line implemented, but an iterative one cannot. For a 夡 This research was supported by the Natural Sciences and Engineering Research Council of Canada and the National Natural Science Foundation of China (No. 60574051 and 60528007). ∗ Corresponding author. Tel.: 86510 88637783; fax: +86510 85910652. E-mail addresses: [email protected], [email protected] (F. Ding), [email protected] (T. Chen). 1 F. Ding is also an Adjunct Professor in the College of Automation at the Nanjing University of Technology, Nanjing 210009, P.R. China.
0167-6911/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.sysconle.2006.10.026
recursive algorithm, new information (input and/or output data) is always used in the algorithm which recursively computes the parameter estimates every step as time increases. Some iterative and/or off-line algorithms of H–W models were discussed in [2,1,4,6,7,22,23,25,28,30,36,11,20] and other recursive and/or on-line algorithms were studied in, e.g., [35,11,20,32,33,3,5,34]. In the identification area of nonlinear systems, Bai reported a two-stage identification algorithm for Hammerstein–Wiener nonlinear systems based on singular value decomposition (SVD) and least squares [1] and studied identification problem of systems with hard input nonlinearities of known structure [2]; Vörös presented a half-substitution algorithm to identify Hammerstein systems with two-segment nonlinearities and with multisegment piecewise-linear characteristics, but no convergence analysis was carried out [32,33]; Cerone and Regruto analyzed parameter error bounds in the Hammerstein models by using the output measurement error bound [6]. Also, Pawlak used the series expansion approach to study the identification of Hammerstein nonlinear outputerror (state-space) models [29]. Recently, an iterative leastsquares and a recursive least-squares identification methods were reported in [11], and an iterative gradient and a recursive
374
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
as follows [6,11,20]:
v(t) u(t)
f (·)
u¯ (t)
B(z) A(z)
x(t)
+
y(t)
u(t) ¯ = f (u(t)) = c1 1 (u(t)) + c2 2 (u(t)) m + · · · + cm m (u(t)) = cj j (u(t)).
(1)
j =1
Fig. 1. The Hammerstein nonlinear output-error system.
Then the Hammerstein nonlinear output-error model in Fig. 1 may be expressed as stochastic gradient algorithms were developed in [20] for nonlinear ARMAX models based only on the available input–output data, and the convergence properties of the recursive algorithms involved were proved. However, most of the identification approaches in the literature assume that the systems under consideration are nonlinear ARX models, or equation-error-like models [35,2,1,7,3]. That is, each element of the information vector consisting of input–output data is measured. In this paper, we focus on the identification problem of a class of Hammerstein output-errortype nonlinear systems and present an auxiliary-model leastsquares (AMLS) algorithm, which is different from the ones mentioned above in that the information vector in our identification model contains unknown variables (namely, unavailable noise-free outputs), and we adopt an auxiliary model or reference model to estimate these unknown variables and further use the outputs of the auxiliary model instead of the unknown noisefree outputs to identify the system parameters. The basic idea is to extend the Landau’s output-error method to study identification problem of nonlinear systems [21,31]. To the best of our knowledge, few publications addressed identification methods of Hammerstein nonlinear output-error systems, especially the convergence problem of the algorithms involved, which are the focus of this work. The objective of this paper is, by means of the auxiliary model identification principle, to derive an algorithm to estimate the system parameters of the nonlinear output-error models based on the available input–output data {u(t), y(t)}, and to study the properties of the algorithm involved. Briefly, the paper is organized as follows. Section 2 describes the identification algorithms related to the Hammerstein systems. Section 3 analyzes the properties of the proposed stochastic algorithm. Section 4 provides an illustrative example to show the effectiveness of the algorithm proposed. Finally, we offer some concluding remarks in Section 5. 2. The algorithm description Consider the Hammerstein output-error system shown in Fig. 1 that consists of a nonlinear memoryless element followed by a linear output-error model [6,29], where the true output (namely, the noise-free output) x(t) and the inner variable u(t) ¯ (namely, the output of the nonlinear block) are unmeasurable, u(t) is the system input, y(t) is the measurement of x(t), v(t) is an additive noise with zero mean. The nonlinear part in the Hammerstein model is a polynomial of a known order in the input [7,22,28], or, more generally, a nonlinear function of a known basis (1 , 2 , . . . , m )
x(t) =
B(z) B(z) u(t) ¯ = [c1 1 (u(t)) + c2 2 (u(t)) A(z) A(z) + · · · + cm m (u(t))],
y(t) = x(t) + v(t). Here, A(z) and B(z) are polynomials in the shift operator z−1 [z−1 y(t) = y(t − 1)] with A(z) = 1 + a1 z−1 + a2 z−2 + · · · + an z−n , B(z) = b1 z−1 + b2 z−2 + b3 z−3 + · · · + bn z−n . Notice that for the Hammerstein model shown in Fig. 1, f (u) and G(z) := B(z)/A(z) are actually not unique. Any pair (f (u), G(z)/) for some nonzero and finite constant would produce identical input and output measurements. In other words, any identification scheme cannot distinguish between (f (u), G(z)) and (f (u), G(z)/). Therefore, to get a unique parameterization, without loss of generality, one of the gains of f (u) and G(z) has to be fixed. There are several ways to normalize the gains [1,6,22]. Here, we adopt the assumption [22,3]: the first coefficient of the function f (·) equals 1; i.e., c1 = 1 [11,20]. Eq. (2) can be rewritten as a recursive form x(t) = − = −
n i=1 n i=1
ai x(t − i) + ai x(t − i) +
n i=1 n
bi u(t ¯ − i) bi
i=1
m
cj j (u(t − i)).
j =1
Define the parameter vector and information vector 0 (t) as ⎡ ⎤ a ⎢ c1 b ⎥ ⎢ ⎥ n0 c b⎥ =⎢ ⎢ 2. ⎥ ∈ R , ⎣ .. ⎦ cm b ⎡ ⎤ −x(t − 1) ⎢ −x(t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 , n0 := (m + 1)n, 0 (t) = ⎢ (2) . ⎢ ⎥ ⎣ −x(t − n) ⎦ (t) a = [a1 , a2 , . . . , an ]T ∈ Rn ,
b = [b1 , b2 , . . . , bn ]T ∈ Rn ,
c = [c2 , c3 , . . . , cm ]T ∈ Rm−1 , (t) = [T1 (t), T2 (t), . . . , Tm (t)]T ∈ Rmn ,
(3)
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
j (t) = [j (u(t − 1)), j (u(t − 2)), n
. . . , j (u(t − n))] ∈ R , T
j = 1, 2, . . . , m.
(4)
Here, the superscript T denotes the matrix transpose. Then we obtain the identification model of the HOE systems: x(t) = T0 (t),
y(t) = T0 (t) + v(t).
(6)
v(t − p + 1) From (5) and (6), we have Y(p, t) = 0 (p, t) + V(p, t).
(7)
Because 0 (p, t) contains unmeasured inner variables x(t − i), i = 1, 2, . . . , n, the standard least-squares method cannot be applied directly to obtain the least-squares estimate [T0 (p, t)0 (p, t)]−1 T0 (p, t)Y(p, t) of the parameter vector by minimizing the cost function Y(p, t) − 0 (p, t)2 even if V(p, t) is a “white” noise vector with zero mean, where X2 = tr[XXT ]. If these unknowns x(t−i) are replaced by the outputs xa (t−i) of an auxiliary (or reference) model, xa (t) = Pa (z)f (u(t)),
According to the auxiliary model identification principle: the unknown variables x(t − i) in 0 (t) are replaced by the output xa (t − i) of the auxiliary model, 0 (t) by (t) and 0 (p, t) by (p, t); then it is easy from (7) to get a recursive pseudo-linear regression least-squares identification algorithm based on the auxiliary model as follows:
(5)
The systems in (5) by parameterization in (2) contain more parameters than actually needed, but the benefit to do this is that a linear regression form can be obtained. However, the disadvantage is that it requires extra computation, which can be approached using the hierarchical estimation methods based on the hierarchical identification principle [8,12–15,18]. From (5), we can see that the information vector 0 (t) consists of unknown noise-free outputs x(t − i) and system inputs. Let p be the data length (p?n0 ), and define ⎡ ⎤ y(t) ⎢ y(t − 1) ⎥ ⎥, Y(p, t) := ⎢ .. ⎣ ⎦ . y(t − p + 1) ⎤ ⎡ T0 (t) ⎢ T (t − 1) ⎥ 0 ⎥ ⎢ 0 (p, t) := ⎢ ⎥, .. ⎦ ⎣ . T (t − p + 1) ⎡ 0 ⎤ v(t) ⎢ v(t − 1) ⎥ ⎥. V(p, t) := ⎢ .. ⎣ ⎦ .
375
or xa (t) = T (t)a ,
then the identification problem of can be solved using xa (t) instead of x(t). Here, Pa (z), (t) and a are the transfer function, information vector and parameter vector of the auxiliary model, respectively. If we use the estimate of B(z)/A(z) as an ˆ auxiliary model Pa (z), namely, take a to be the estimate (t) of , and (t) to be the regressive vector of xa (t) and u(t), and use (t) as 0 (t), then the identification algorithms based on this idea are called the auxiliary (or reference) model identification method or output-error method [21,31]. Of course, there are other ways to choose auxiliary models, e.g., using the finite impulse response model [9,10].
ˆ = (t ˆ − 1) + P(t)T (p, t)[Y(p, t) − (p, t)(t ˆ − 1)], (t) (8) P−1 (t) = P−1 (t − 1) + T (p, t)(p, t), ˆ xa (t − i) = T (t − i)(t), ⎡ ⎢ (p, t) = ⎢ ⎣
T (t) T (t − 1) .. .
P(0) = p0 I,
i = 0, 1, 2, . . . , n,
(9) (10)
⎤ ⎥ ⎥, ⎦
T (t − p + 1) ⎤ −xa (t − 1) ⎢ −xa (t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 . (t) = ⎢ . ⎢ ⎥ ⎣ −x (t − n) ⎦ a (t) ⎡
(11)
Here, p0 is a large positive number, e.g., p0 =106 , and we refer to xa (t) as the estimate of x(t). Eqs. (8)) to (11) is also called the AMLS identification algorithm of estimating in (7) of the Hammerstein output-error model, HOE-AMLS algorithm for short. Note that the algorithm is combined in the sense the parameters and true outputs are estimated simultaneously, and can be implemented on-line. Note c1 = 1, the estimates aˆ = [aˆ 1 , aˆ 2 , . . . , aˆ n ]T and bˆ = [bˆ1 , bˆ2 , . . . , bˆn ]T of a and b can be read from the first and ˆ respectively. Let ˆ i be the ith element second n entries of , ˆ referring to the definition of , then the estimates of cj , of , j = 2, 3, . . . , m, may be computed by cˆj = ˆ j n+i /bˆi ,
j = 2, 3, . . . , m; i = 1, 2, . . . , n.
Since we do not need such n estimates cˆj , we may calculate the c-parameters by a least-squares fit or take their average like in [11,20] as the estimate of cj , i.e., cˆj =
n 1 ˆ j n+i , j = 2, 3, . . . , m. n bˆi i=1
3. The performance analysis Let us introduce some notation first. The symbol I stands for an identity matrix of appropriate sizes; |X| = det[X] represents the determinant of the square matrix X; max [X] and min [X] represent the maximum and minimum eigenvalues of the symmetric matrix X, respectively; for g(t) 0, we write f (t) = O(g(t)) if there exists positive constants 1 and t0 such that |f (t)| 1 g(t) for t t0 .
376
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
Dividing (13) by [ln |P−1 (t)] and summing for t give
Define P−1 0 (t) =
t i=1 −1
r(t) = tr[P
T0 (p, i)0 (p, i) + (t)];
∞ T (t − i)P(t)(t − i)
1 I; p0
Hence, we easily get (12)
Lemma 1. For each i (i = 0, 1, . . . , p − 1), the following inequalities hold:
2
(j − i)P(j )(j − i) ln |P T (t − i)P(t)(t − i)
t=1
[ln |P−1 (t)|]
< ∞,
−1
(t)| + n0 ln p0 , a.s.,
a.s., for any > 1.
P−1 (t − 1) = P−1 (t) − T (p, t)(p, t) P−1 (t) − (t − i)T (t − i) = P−1 (t)[I − P(t)(t − i)T (t − i)]. Taking determinants on both sides and using the formula det[I+ DE] = det[I + ED] yields |P−1 (t − 1)| |P−1 (t)||I − P(t)(t − i)T (t − i)|
˜ = (t) ˆ − , (t)
(14) (15)
Lemma 2. For the system in (7) and the HOE-AMLS algorithm in (8)–(11), assume the noise sequence {v(t)} with zero mean and bounded time-varying variance satisfies [24]:
where {v(t), Ft } is a martingale sequence defined on a probability space {, F, P } and {Ft } is the algebra sequence generated by {v(t)}. Then the following inequality holds: E[W (t) + S(t)|Ft−1 ] W (t − 1) + S(t − 1) p−1
T (t − i)P(t)(t − i)¯ 2v ,
a.s.,
i=0
Hence
where
T (t − i)P(t)(t − i)
|P−1 (t)| − |P−1 (t − 1)| |P−1 (t)|
.
(13)
Replacing t with j and summing for j from 1 to t yields (noting that |P−1 (t)| is a nondecreasing function of t)
− 1)|
|P−1 (j )|
j =1 t |P−1 (j )|
−1 j =1 |P (j −1)|
a.s.,
(16)
i=1
˜ = Y(t)
1 2
˜ (p, t)(t) (17) (18)
Here, (A3) guarantees that S(t) 0. The proof can be done in a similar way in [11,19] and is omitted here. Next, we give prove the main results of this paper.
dx |P−1 (j )|
|P−1 (t)| dx = ln |P−1 (t)| − ln |P−1 (0)| −1 |P (0)| x = ln |P−1 (t)| + n0 ln p0 ,
˜ T (i)Y(i), ˜ U
˜ ˜ U(t) = −(p, t)(t).
j =1
|P−1 (j )| − |P−1 (j
S(t) = 2
t
ˆ − V(p, t)], + [Y(p, t) − (p, t)(t)
T (j − i)P(j )(j − i) t
˜ W (t) = ˜ (t)P−1 (t)(t).
+ 2p
= |P−1 (t)|[1 − T (t − i)P(t)(t − i)].
=
˜ and a nonDefine the parameter estimation error vector (t) negative definite function W (t) as
(A1) E[v(t)|Ft−1 ] = 0, a.s., (A2) E[v 2 (t)|Ft−1 ] = 2v (t) ¯ 2v < ∞, a.s., 1 − 21 is strictly positive real, (A3) H (z) = A(z)
Proof. From the definition of P(t) in (9), we have
|P−1 (∞)|
−1 1
= −1 −1 −1 − 1 x(ln x) (ln x) |P (0)| |P (0)| 1 1 1 − = − 1 [ln |P−1 (0)|]−1 [ln |P−1 (∞)|]−1 < ∞, a.s. dx
T
T
j =1 ∞
t
|P−1 (t)|[ln |P−1 (t)|]
|P−1 (∞)|
The stochastic process theory [27,16] and the stochastic martingale theory [11,10,17,24] are the main tool of analyzing performance of identification algorithms. Here, we prove the main convergence results of the HOE-AMLS algorithm by using the martingale convergence theorem. To do this, some mathematical preliminaries are required.
1
∞ |P−1 (t)| − |P−1 (t − 1)| t=1
|P−1 (t)|r n0 (t); r(t) n0 max [P−1 (t)]; ln |P−1 (t)| = O(ln r(t)).
t
[ln |P−1 (t)|]
t=1
r0 (t) = tr[P−1 0 (t)].
a.s.
Theorem 1. For the system in (7), assume that (A1)–(A3) hold, and A(z) is stable, i.e., all zeros of A(z) are inside the unit circle. Then for any > 1, the parameter estimation error by
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
the HOE-AMLS algorithm in (8)–(11) satisfies:
[ln r0 (t)] 2 ˆ , a.s. (t) − = O min [P−1 0 (t)]
(22), we have
min [P−1 0 (t)] + [ln r0 (t)] [ln r0 (t)] =O , a.s. for any > 1. min [P−1 0 (t)]
(19)
Let [ln |P−1 (t)|]
.
Since ln |P−1 (t)|] is nondecreasing, according to Lemma 2, we have E[Z(t)|Ft−1 ] Z(t − 1) + 2p
p−1
T (t − i)P(t)(t − i)
i=0
[ln |P−1 (t)|]
¯ 2v ,
a.s. (20)
Using Lemma 1 and applying the martingale convergence theorem (Lemma D.5.3 in [24]) to (20), we conclude that Z(t) converges a.s. to a finite random variable, say, Z0 ; i.e., Z(t) =
W (t) + S(t) [ln |P−1 (t)|]
→ Z0 < ∞,
a.s.,
W (t) = O([ln |P−1 (t)|] ), a.s., S(t) = O([ln |P−1 (t)|] ), a.s.
(21)
Since H (z) is a strictly positive real function, from the definition of S(t) and referring to [11], we have 2 ˜ U(i) = O([ln |P−1 (t)|] ).
i=1
(i) − V(p, i)2 k1
t
ˆ = (t ˆ − 1) + P(t)(t)[y(t) − T (t)(t ˆ − 1)], (t) P(t − 1)(t)T (t)P(t − 1) , 1 + T (t)P(t − 1)(t)
ˆ xa (t) = T (t)(t), ⎡ ⎤ −xa (t − 1) ⎢ −xa (t − 2) ⎥ ⎢ ⎥ .. ⎥ ∈ Rn0 . (t) = ⎢ . ⎢ ⎥ ⎣ −x (t − n) ⎦ a (t)
(24) (25) (26)
(27)
B(z) u(t) ¯ + v(t), A(z) A(z) = 1 + a1 z−1 + a2 z−2 = 1 − 1.60z−1 + 0.80z−2 , y(t) =
(22)
2 ˜ U(i) + k2
i=1
= O([ln |P−1 (t)|] ) = O([ln r(t)] ).
From Theorem 1, we can see that if [ln r0 (t)] = o(min [P−1 0 (t)]) (namely, the generalized persistent excitation (PE) condition [10], the weak PE condition [10] or the strong PE condition [16] holds), then the estimation error (t) − converges to zero. This requires that G(z) is stable, u(t) ¯ is a PE signal, and so is u(t) = f −1 (u(t)). ¯ The positive real condition in (A3) depends on the unknown model parameters and is a more restrictive assumption than A(z) being stable [27]. ˆ − 1) ∈ R1 is called the As p = 1, e(t) := y(t) − T (t)(t ˆ − 1)] ∈ innovation [26], and E(t) := [Y(p, t) − (p, t)(t Rp is referred as to the innovation vector. p here may be also known as the innovation length, the algorithm in (8)–(11) is also called the auxiliary model multi-innovation least-squares identification algorithm for Hammerstein output-error systems. When p = 1, we get a simple HOE-AMLS algorithm:
An example is given to demonstrate the effectiveness of the proposed algorithms. Consider the following system:
Like in [11,19], it follows that there exist positive constants k1 and k2 such that t
4. Example
i=1
From (19), (21) and (12), we have
−1 (t)|] [ln |P 2 ˜ (t) =O min [P−1 (t)]
[ln r(t)] =O , a.s. for any > 1. min [P−1 (t)]
This proves Theorem 1.
P(t) = P(t − 1) −
or
t
T ˜ ˜ (t)P−1 (t)(t) W (t) 2 ˜ (t) = . −1 min [P (t)] min [P−1 (t)]
W (t) + S(t)
[ln r0 (t)]
ˆ − 2 = O (t)
Proof. From the definition of W (t), we have
Z(t) =
377
(23)
Using the similar way in [11,19], it is easy to get r(t)=O(r0 (t)), min [P−1 (t)] = O(min [P−1 0 (t)] + O([ln r0 (t)] ). Thus, from
B(z) = b1 z−1 + b2 z−2 = 0.85z−1 + 0.65z−2 , u(t) ¯ = f (u(t)) = c1 u(t) + c2 u2 (t) + c3 u3 (t) = u(t) + 0.5u2 (t) + 0.25u3 (t), s = [a1 , a2 , b1 , b2 , c2 , c3 ]T . {u(t)} is taken as a PE signal sequence with zero mean and unit variance 2u = 1.002 , and {v(t)} as a white noise sequence with zero mean and constant variance 2v = 0.502 and 2v = 2.002 . Apply the HOE-AMLS algorithm with p = 1 to estimate the parameters of this system, the parameter estimates = [ 1 , 2 , . . . , 8 ]T and s and their errors with different
378
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
Table 1 The estimates of (2v = 0.502 ) t
1
2
3
4
5
6
7
8
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.56901 −1.57996 −1.57992 −1.59002 −1.58361 −1.59141 −1.59375 −1.59131 −1.59376
0.77630 0.78122 0.78325 0.78999 0.78623 0.79217 0.79512 0.79218 0.79477
0.68233 0.79989 0.83494 0.87510 0.85321 0.83909 0.86519 0.85338 0.85300
1.04526 0.91609 0.85387 0.79772 0.71349 0.69136 0.69088 0.68343 0.66855
0.37656 0.35609 0.36760 0.37839 0.37916 0.37660 0.38268 0.38296 0.38739
0.40992 0.37159 0.39116 0.39249 0.39854 0.38461 0.38439 0.38079 0.37451
0.21864 0.22293 0.20882 0.19782 0.20294 0.21061 0.20713 0.21037 0.21052
0.13271 0.15234 0.15149 0.15553 0.17042 0.17118 0.16963 0.16826 0.16821
20.43326 13.14011 10.34029 7.93938 5.08542 4.10887 3.95682 3.62598 3.03148
True values
−1.60000
0.80000
0.85000
0.65000
0.42500
0.32500
0.21250
0.16250
Table 2 The estimates of (ai , bi , ci ) (2v = 0.502 ) t
a1
a2
b1
b2
c2
c3
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.56901 −1.57996 −1.57992 −1.59002 −1.58361 −1.59141 −1.59375 −1.59131 −1.59376
0.77630 0.78122 0.78325 0.78999 0.78623 0.79217 0.79512 0.79218 0.79477
0.68233 0.79989 0.83494 0.87510 0.85321 0.83909 0.86519 0.85338 0.85300
1.04526 0.91609 0.85387 0.79772 0.71349 0.69136 0.69088 0.68343 0.66855
0.47202 0.42540 0.44918 0.46220 0.50149 0.50256 0.49934 0.50296 0.50717
0.22370 0.22250 0.21376 0.21051 0.23835 0.24930 0.24247 0.24636 0.24920
20.05581 13.13781 9.97809 7.41968 3.15539 2.05765 2.08352 1.66260 1.00617
True values
−1.60000
0.80000
0.85000
0.65000
0.50000
0.25000
Table 3 The estimates of (2v = 2.002 ) t
1
2
3
4
5
6
7
8
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.54443 −1.55876 −1.55951 −1.58706 −1.56785 −1.58573 −1.59032 −1.58242 −1.59036
0.75677 0.76487 0.76570 0.78914 0.77375 0.78720 0.79314 0.78315 0.79298
0.35049 0.67966 0.77079 0.85335 0.81059 0.77961 0.87038 0.83949 0.84170
1.37101 1.14544 1.03161 0.93939 0.78027 0.76794 0.77592 0.76061 0.71035
0.41764 0.33138 0.37081 0.35883 0.33976 0.32739 0.34143 0.33394 0.34846
0.43136 0.36260 0.42585 0.43890 0.47286 0.43761 0.44251 0.43979 0.42285
0.25470 0.26717 0.21899 0.18798 0.20004 0.22060 0.20796 0.21446 0.21475
0.10459 0.13651 0.12295 0.13346 0.17329 0.16684 0.16391 0.15866 0.16187
41.00730 24.88620 18.98097 14.78786 10.27459 9.39654 8.89719 8.55583 6.40464
True values
−1.60000
0.80000
0.85000
0.65000
0.42500
0.32500
0.21250
0.16250
Table 4 The estimates of (ai , bi , ci ) (2v = 2.002 ) t
a1
a2
b1
b2
c2
c3
(%)
100 200 300 500 1000 1500 2000 2500 3000
−1.54443 −1.55876 −1.55951 −1.58706 −1.56785 −1.58573 −1.59032 −1.58242 −1.59036
0.75677 0.76487 0.76570 0.78914 0.77375 0.78720 0.79314 0.78315 0.79298
0.35049 0.67966 0.77079 0.85335 0.81059 0.77961 0.87038 0.83949 0.84170
1.37101 1.14544 1.03161 0.93939 0.78027 0.76794 0.77592 0.76061 0.71035
0.75310 0.40206 0.44694 0.44386 0.51259 0.49490 0.48129 0.48800 0.50463
0.40149 0.25613 0.20165 0.18118 0.23443 0.25011 0.22509 0.23203 0.24150
43.00408 24.82545 18.52710 14.04942 6.65812 6.43023 6.10921 5.36474 2.91094
True values
−1.60000
0.80000
0.85000
0.65000
0.50000
0.25000
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
379
1
0.8
0.6 2
v = 0.502
0.4
0.2
2
v = 2.002
0 0
500
1000
1500
2000
2500
3000
2500
3000
t Fig. 2. The parameter estimation errors vs. t.
1
0.8
s
0.6 2
= 0.502 v
0.4
0.2
0
2v = 2.002 0
500
1000
1500
2000
t Fig. 3. The parameter estimation errors s vs. t.
noise variances are shown in Tables 1–4, and the parameter estimation errors and s vs. t are shown in Figs. 2 and 3, ˆ =(t)−/, s =ˆ s (t)−s /, ˆ s (t), is the estimate of s . From Tables 1–4 and Figs. 2 and 3, we can draw the following conclusions:
Hammerstein output-error models. The analysis using the martingale convergence theorem indicates that the proposed algorithm can give consistent parameter estimation.
• A high noise level results in a slow rate of convergence of the parameter estimates to the true parameters. • It is clear that the errors and s are becoming smaller (in general) as t increases. This confirms the proposed theorem.
[1] E.W. Bai, An optimal two-stage identification algorithm for Hammerstein-Wiener nonlinear systems, Automatica 34 (3) (1998) 333–338. [2] E.W. Bai, Identification of linear systems with hard input nonlinearities of known structure, Automatica 38 (5) (2002) 853–860. [3] E.W. Bai, A blind approach to the Hammerstein–Wiener model identification, Automatica 38 (6) (2002) 967–979. [4] E.W. Bai, A random least-trimmed-squares identification algorithm, Automatica 39 (9) (2003) 1561–1659. [5] M. Boutayeb, M. Darouach, Recursive identification method for MISO Wiener–Hammerstein model, IEEE Trans. Automat. Control 40 (2) (1995) 287–291.
5. Conclusions A recursive AMLS algorithm based on replacing unavailable variables (noise-free outputs) by their estimates is derived for
References
380
F. Ding et al. / Systems & Control Letters 56 (2007) 373 – 380
[6] V. Cerone, D. Regruto, Parameter bounds for discrete-time Hammerstein models with bounded output errors, IEEE Trans. Automat. Control 48 (10) (2003) 1855–1860. [7] F. Chang, R. Luus, A noniterative method for identification using Hammerstein model, IEEE Trans. Automat. Control 16 (5) (1971) 464–468. [8] F. Ding, T. Chen, Hierarchical gradient-based identification of multivariable discrete-time systems, Automatica 41 (2) (2004) 315–325. [9] F. Ding, T. Chen, Identification of dual-rate systems based on finite impulse response models, Internat. J. Adapt. Control Signal Process. 18 (7) (2004) 589–598. [10] F. Ding, T. Chen, Combined parameter and output estimation of dualrate systems using an auxiliary model, Automatica 40 (10) (2004) 1739–1748. [11] F. Ding, T. Chen, Identification of Hammerstein nonlinear ARMAX systems, Automatica 41 (9) (2005) 1479–1489. [12] F. Ding, T. Chen, Hierarchical least squares identification methods for multivariable systems, IEEE Trans. Automat. Control 50 (3) (2005) 397–402. [13] F. Ding, T. Chen, Hierarchical identification of lifted state-space models for general dual-rate systems, IEEE Trans. Circuits Systems I. Regular Papers 52 (6) (2005) 1179–1187. [14] F. Ding, T. Chen, Iterative least squares solutions of coupled Sylvester matrix equations, Systems Control Lett. 54 (2) (2005) 95–107. [15] F. Ding, T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Trans. Automat. Control 50 (8) (2005) 1216–1221. [16] F. Ding, T. Chen, Performance bounds of forgetting factor least squares algorithm for time-varying systems with finite measurement data, IEEE Trans. Circuits Systems I. Regular Papers 52 (3) (2005) 555–566. [17] F. Ding, T. Chen, Parameter estimation of dual-rate stochastic systems by using an output error method, IEEE Trans. Automat. Control 50 (9) (2005) 1436–1441. [18] F. Ding, T. Chen, On iterative solutions of general coupled matrix equations, SIAM J. Control Optim. 44 (6) (2006) 2269–2284. [19] F. Ding, H.B. Chen, M. Li, Multi-innovation least squares identification methods based on the auxiliary model for MISO systems, Appl. Math. Comput., 2007, in press. [20] F. Ding, Y. Shi, T. Chen, Gradient-based identification algorithms for Hammerstein nonlinear ARMAX models, Nonlinear Dynam. 45 (1–2) (2006) 31–43.
[21] L. Dugard, I.D. Landau, Recursive output error identification algorithms theory and evaluation, Automatica 16 (5) (1980) 443–462. [22] P.G. Gallman, A comparison of two Hammerstein model identification algorithms, IEEE Trans. Automat. Control 21 (1) (1976) 124–126. [23] F. Giri, F.Z. Chaoui, Y. Rochdi, Parameter identification of a class of Hammerstein plants, Automatica 37 (5) (2001) 749–756. [24] G.C. Goodwin, K.S. Sin, Adaptive Filtering, Prediction and Control, Prentice-Hall, Englewood Cliffs, NJ, 1984. [25] N.D. Haist, F. Chang, R. Luus, Nonlinear identification in the presence of correlated noise using a Hammerstein model, IEEE Trans. Automat. Control 18 (5) (1973) 553–555. [26] L. Ljung, System Identification: Theory for the User, second ed., Prentice-Hall, Englewood Cliffs, NJ, 1999. [27] L. Ljung, T. Söderström, Theory and Practice of Recursive Identification, MIT Press, Cambridge, MA, 1983. [28] K.S. Narendra, P.G. Gallman, An iterative method for the identification of nonlinear systems using a Hammerstein model, IEEE Trans. Automat. Control 11 (3) (1966) 546–550. [29] M. Pawlak, On the series expansion approach to the identification of Hammerstein system, IEEE Trans. Automat. Control 36 (6) (1991) 763–767. [30] P. Stoica, On the convergence of an iterative algorithm used for Hammerstein system identification, IEEE Trans. Automat. Control 26 (4) (1981) 967–969. [31] P. Stoica, T. Söderström, Analysis of an output error identification algorithm, Automatica 17 (6) (1981) 861–863. [32] J. Vörös, Iterative algorithm for parameter identification of Hammerstein systems with two-segment nonlinearities, IEEE Trans. Automat. Control 44 (11) (1999) 2145–2149. [33] J. Vörös, Modeling and parameter identification of systems with multisegment piecewise-linear characteristics, IEEE Trans. Automat. Control 47 (1) (2002) 184–189. [34] T. Wigren, Convergence analysis of recursive identification algorithms based on the nonlinear Wiener model, IEEE Trans. Automat. Control 39 (11) (1994) 2191–2206. [35] T. Wigren, A.E. Nordsjö, Compensation of the RLS algorithm for output nonlinearities, IEEE Trans. Automat. Control 44 (10) (1999) 1913–1918. [36] Y. Zhu, Estimation of an N-L-N Hammerstein–Wiener model, Automatica 38 (9) (2002) 1607–1614.
