Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm

ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

Research article

Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm Feng Yu n, Zhizhong Mao, Ping Yuan, Dakuo He, Mingxing Jia Department of Control Theory and Control Engineering, Northeastern University, No.11, Lane 3, WenHua Road, HePing District, Shenyang, China

art ic l e i nf o

a b s t r a c t

Article history: Received 15 May 2016 Received in revised form 5 January 2017 Accepted 21 May 2017

This paper focuses on the recursive parameter estimation for the single input single output Hammerstein-Wiener system model, and the study is then extended to a rarely mentioned multiple input single output Hammerstein-Wiener system. Inspired by the extended Kalman filter algorithm, two basic recursive algorithms are derived from the first and the second order Taylor approximation. Based on the form of the first order approximation algorithm, a modified algorithm with larger parameter convergence domain is proposed to cope with the problem of small parameter convergence domain of the first order one and the application limit of the second order one. The validity of the modification on the expansion of convergence domain is shown from the convergence analysis and is demonstrated with two simulation cases. & 2017 ISA. Published by Elsevier Ltd. All rights reserved.

Keywords: Hammerstein-Wiener system Recursive parameter estimation Extended Kalman filter algorithm Multiple input single output Uniform convergence

1. Introduction Identification of a class of nonlinear dynamical systems, which are so-called block oriented systems, gradually becomes the focus of research on the nonlinear system identification. This is because these kinds of systems widely exist in the industrial areas, such as chemical process, electronic circuits, biologic systems, etc. Two simplest examples in this class of systems are Hammerstein (H) systems, which consist of a static nonlinear block followed by a linear dynamic one; and Wiener (W) systems, which are composed by the inverse configuration of H systems. Many efforts have been devoted to identify these two kinds of nonlinear block oriented systems [1–9]. With identification algorithms for H and W systems gradually reaching maturity, research on identifying the more complex Hammerstein-Wiener (H-W) systems attracts more attention [10]. The structure of single input single output (SISO) H-W systems is composed with two static nonlinear blocks and a linear dynamic block in between. The first study focused on identifying the H-W systems may be carried out by [11]. In the paper, an approximate model was used to express the H-W systems and a two-stage identification algorithm based on singular value decomposition and least squares method was proposed to identify the unknown n

Corresponding author. E-mail addresses: [email protected] (F. Yu), [email protected] (Z. Mao), [email protected] (P. Yuan), [email protected] (D. He), [email protected] (M. Jia).

parameters iteratively. Strictly speaking, as some important terms were omitted, the approximate model used in [11] cannot completely reflect the original system characteristics. After this, other methods were proposed to identify the H-W systems, for example, the relaxed iterative method [12], frequency domain method [13], maximum likelihood estimation method [14,15], iterative least squares method [16,17], and so on. However, all these kinds of methods cannot be used for online identification because of their heavy computation load. Therefore, some recursive methods were proposed. [18] and [19] extended the work [11] to the recursive identification of noisy H-W systems. However, the similar approximate model was still used in these works. Also based on [11,20] proposed a blind method to identify H-W systems. This algorithm works only under special input signals. The author also didn’t discuss the convergence of the algorithm and the interference of noise. [21] proposed a hierarchical least squares estimation algorithm to estimate the parameters within each block of the H-W systems, respectively. [22] worked on recursive identification of the H-W system with a dead-zone input block. However, [21] didn’t analyze the convergence of the algorithm and [22] didn’t consider the effect of noise on the algorithm. In addition, all these algorithms include two defects. Firstly, these methods only focus on the SISO systems. Even for the multiple input single output (MISO) case, these algorithms cannot work. Secondly, the assumption of irreversibility for the system output nonlinearity is necessary. That means even if the output nonlinearity possesses a simple quadratic function characteristic, the algorithm is no longer applicable. [23] used an extended Kalman filter (EKF) kind

http://dx.doi.org/10.1016/j.isatra.2017.05.012 0019-0578/& 2017 ISA. Published by Elsevier Ltd. All rights reserved.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

Nomenclature a b c e f m p q u v w x y A B C G I K O P R2 U V W

ε

model parameter of linear part model parameter of input nonlinear part model parameter of output nonlinear part model estimation error nonlinear characteristic order of FIR filter the basis function number of input nonlinear part the basis function number of output nonlinear part system input output of input nonlinear part input of output nonlinear part true system output system output measurement vector form of a vector form of b vector form of c linear characteristic identity matrix update gain Lagrange remainder term covariance matrix of the parameter estimation error estimate of δ2 vector form of u vector form of v vector form of w measurement noise

algorithm to identify H-W system recursively. Although the usually used irreversible assumption was not necessary, the initial parameter estimates close to the true values were required. In summary, the research on recursive identification of H-W systems is still in the primary stage and much work remains to be done for this. The aim of this paper is to propose a recursive algorithm for online estimating the parameters of the H-W systems with measurement noise and extend it to the MISO case. Based on our previous work [23], a modified EKF-kind algorithm(MEKF) is proposed, which is derived from the two basic algorithms with the first and the second order Taylor approximation, respectively. The convergence of the proposed algorithm is studied and a sufficient condition for achieving the uniform convergence of the parameter estimation is obtained from the convergence analysis. The validity of this algorithm is demonstrated with two simulation examples, including a practical wind tunnel system case. The contributions of this paper are summarized as follows: 1. The parameterized model studied in the paper is not an approximate one. Any items in the model are not omitted. Therefore, the parameterize model we used can completely reflect the dynamic characteristics of the system. 2. The proposed algorithm can be used for the MISO H-W systems, which, to the best of our knowledge, have been rarely mentioned in the previous literatures. 3. The usually used irreversible assumption for the output nonlinearity is no longer necessary for the proposed algorithm. 4. The main modification, namely add the parameter λ into the traditional algorithm, is very simple, but it greatly enlarges the convergence domain of the algorithm so that the algorithm can tolerate a larger initial parameter estimation error. In addition, the introduction of λ also strengthens the robustness of the

λ λ(
) φ(
) θ δ Λ Φ 1

algorithm variable basis function of input nonlinear part basis function of output nonlinear part model parameter vector variance vector form of λ vector form of φ positive infinity

Operators tr E ℱ -

trace of a matrix expectation s function sequence generated by measurable y approach to

Subscripts t I O

sampling time input part output part

Superscripts T ∧ 

transposition estimate estimation error

algorithm on disturbance.

the

unmodeled

dynamic

and

unmodeled

The outline of the paper is as follows. A nomenclature is shown at the end of the introduction section. The H-W system and its parameterized representation are presented in Section 2. The derivation of the MEKF is given in Section 3. The convergence property is analyzed in Section 4. Section 5 discusses 3 essential problems closely related to the algorithms. Section 6 shows the results of the two simulation cases. Finally, the conclusions are drawn and further research direction is discussion in Section 7.

2. Problem formulation 2.1. Parameterization of SISO H-W systems Consider the SISO H-W systems shown in Fig. 1, where ut is the system input; the internal variables vt and wt are unmeasurable; xt is the true system output, namely the noise-free system output; yt is the measurement of xt; εt is the measurement noise; G(vt) is the linear block transfer function; fI(ut) and fO(wt) are the functions of input and output nonlinear blocks, respectively. The input and output nonlinear static blocks are modeled as a linear combination of the known basis functions λi and φi,

Fig. 1. SISO H-W system with output errors.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

respectively: p

vt = fI (ut ) =

∑ biλi( ut ) = B¯T Λ( ut ) i=1

q

xt = fO (wt ) =

∑ ciφi( wt ) = C¯

T

Φ( wt )

i=1

(1)

(2)

where p and q are the numbers of the basis functions; T T B¯ = ⎡⎣ b1, b2 , ⋯ , bp⎤⎦ and C¯ = ⎡⎣ c1, c2, ⋯ , cq⎤⎦ are coefficient vectors; T T Λ( ⋅) = ⎡⎣ λ1( ⋅) , λ2( ⋅) , ⋯ , λ p( ⋅)⎤⎦ ; Φ( ⋅) = ⎡⎣ φ1( ⋅) , φ2( ⋅) , ⋯ , φq( ⋅)⎤⎦ . The linear dynamic block is modeled by a FIR filter:

m

wt = G(vt ) =

∑ aivt − i = AT Vt − 1 i=1

(3)

T where m is the order of the FIR filter; A = ⎡⎣ a1, a2, ⋯ , am⎤⎦ is the T coefficient vector of the linear block; Vt − 1 = ⎡⎣ vt − 1, vt − 2, ⋯ , vt − m⎤⎦ . With (1), (2), and (3), the parameterized model of the SISO H-W systems shown in Fig. 1 can be represented as:

yt = f ( Ut − 1, θ¯ ) + εt

(4)

T T T⎤ T ⎡ where Ut − 1 = ⎡⎣ ut − 1, ut − 2 , ⋯ , ut − m⎤⎦ ; θ¯ = ⎢⎣ AT , B¯ , C¯ ⎥⎦ .

3

vectors Ai, B¯i , and C¯ i are not unique due to the unmeasurable intermediate signals vi,t and wi,t. Therefore, arbitrary gains may be distributed among the linear block and the nonlinear blocks. A simple but effective solution to this problem is setting bi,1 and ci,1 to be a constant value of 1, not participating in parameter updating [20,24]. As a result, the parameter vectors B¯i and C¯ i can be reT T and written as B¯i = ⎡⎣ 1, bi,2 /bi,1, ⋯ , bi, p /bi,1⎤⎦ = ⎡⎣ 1, BiT ⎤⎦ i T T C¯ i = ⎡⎣ 1, ci,2/ci,1, ⋯ , ci, q /ci,1⎤⎦ = ⎡⎣ 1, CiT ⎤⎦ . It is obvious that only the i T updating of the parameter vector θ = ⎡⎣ A1T , B1T , C1T , ⋯ArT , BrT , CrT ⎤⎦ is required. Furthermore, the following assumptions are required for the derivation of the algorithm.

Assumption 1. The orders of the transfer function for the linear blocks mi and the numbers of the basis functions for the nonlinear blocks pi and qi are known. This implies that the estimation of system structure is not needed. Otherwise the identification will become a problem of structural estimation [20]. Assumption 2. The second derivative of the output nonlinearity basis functions exists. Assumption 3. The linear block of the system is stable so that it can be represented as a FIR filter. Assumption 4. The sequence {εt} is a stochastic noise with zero mean and bounded variance of δ2, namely:

2.2. Parameterization of MISO H-W systems The general structure of the MISO H-W systems with r inputs considered in this paper is shown in Fig. 2. Based on the parameterization model of SISO H-W systems, we add the subscript i ( i = 1, ⋯ , r ) to describe the symbols pertaining to the ith subsystem of the MISO H-W system. The output of the MISO H-W system equals to the sum of the outputs of all subsystems, namely: r

xt = f ( Ut − 1, θ¯ ) =

∑ x i, t i=1

(5)

yt is the measurement of xt and is disturbed by the measurement noise εt. Therefore, the parameterized model of the MISO H-W system can also be expressed by (4), while θ¯ is redefined as T T ⎤ T ⎡ T ¯ T ¯T T ¯ T ¯T T ⎢⎣ A1 , B1 , C1 , A2 , B2 , C2 , ⋯ , Ar , B¯r , C¯ r ,⎥⎦ and Ut-1 is redefined as

⎡ U T , U T , ⋯ , U T ⎤T . ⎣ 1, t − 1 2, t − 1 r , t − 1⎦ The objective of this paper is to establish a recursive identification algorithm to estimate the parameter vector θ¯ in the parameterized H-W system model (4) online by using the available input-output data Ut and yt, and to study the properties of the algorithms involved.

3. Recursive estimation algorithm Before the derivation of the algorithm, the problem of parameter estimation consistency should be clarified. That is, the

E⎡⎣εt t − 1⎤⎦ = 0 and E⎡⎣εt2 t − 1⎤⎦ = δ 2 < ∞ where {ℱt  1} is the s function sequence generated by measurable {yt  1} obtained from the deterministic input sequence {ut  1}. Remark 1. According to Assumption 2, only the differentiability of the basis function for the output block is required, which means the input block can be represented by a linear combination of nondifferentiable basis functions. Therefore, the proposed algorithm can be used to identify the H-W systems with special input blocks, such as the input block containing preload nonlinearity, piecewise-linear nonlinearity or a linear combination of other nondifferentiable basis functions. Remark 2. The assumption on the invertibility of output nonlinearity, which is commonly used in the most of the studies on the H-W system identification, is no longer necessary. This is because the linear block function Gi(vi,t) is assumed to be represented as a FIR filter so that the internal variables wi,t does not appear in the right-hand side of the linear block function (3). Otherwise, the inverse function of output nonlinearity is required to estimate internal variable wi,t. This kind of representation of the linear block is usually utilized in the H or W systems identification [25–28]. Based on the assumptions above, we can give the deviation of the algorithm as below. Define the structure of algorithm with the recursive form of EKF:

θ^t = θ^t − 1 + Ktet

(

(6)

Fig. 2. MISO H-W system with output errors.

Pt = Pt − 1 − KtPθTx

−1

)

Kt = Pθx Pxx + R2

(7)

(8)

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

where θ^t is the estimate of θ ; Kt is the update gain; R2 is the estimate of δ2; et is the model estimation error; Pt is the covariance matrix of the parameter estimation errors. We have:

approach is the key to tolerate a bad θ^0 . A straightforward approach is to use the second order Taylor expansion to approximate the system parameterization model, namely:

et = yt − y^t

T 1 f ( θ ) ≈ f θ^t − 1 + f t′−T 1 θ − θ^t − 1 + θ − θ^t − 1 f t"− 1 θ − θ^t − 1 2

(9)

(10)

where y^t is the system estimation; θ˜t is the parameter estimation error, defined as:

θ˜t = θ − θ^t

(11)

Similar with the derivation of the EKF algorithm, we have the following equations under the Gaussian approximation assumption:

∫R

n

f ( τ )φ( τ )dτ

∫R

)

(20)

2

2

t−1

t−1

Therefore, we have:

∫R

≈

n

T ⎡ ^ ⎤ 1 ^ ′T τ − θ^t − 1 f t″− 1 τ − θ^t − 1 ⎥ ⎢ f θt − 1 + f t − 1 τ − θt − 1 + ⎣ ⎦ 2

( )

(

) (

) (

)

φ( τ )dτ 1 =f θ^t − 1 + tr f t″− 1 Pt − 1 2

( )

(

)

(21)

(12)

( τ − θ^ )⎡⎣⎢ f ( τ) − f ( θ^ )⎤⎦⎥φ( τ)dτ t−1

n

) (

y^t,2 = x^t,2

Pθx,2 ≈ Pθx =

( ) ( = ∂ f ( θ^ )/∂θ^ .

where f t″− 1

⎡ ~ ~T ⎤ Pt = E⎢θtθt t − 1⎥ ⎣ ⎦

y^t = x^t =

( )

t−1

(13)

∫R

⎡

n

T

⎤

(τ − θ^ )⎢⎣ f (τ − θ^ ) + 21 (τ − θ^ ) f (τ − θ^ )⎥⎦ t−1

′T t−1

t−1

t−1

″ t−1

t−1

φ( τ )dτ =f t′− 1 Pt − 1 = Pθx,1

Pxx =

∫R

⎡ ⎤ ^ ⎢⎣ f ( τ ) − f θt − 1 ⎥⎦ φ( τ )dτ

( )

n

(14)

where n is the number of parameters to be identified; φ( τ ) is defined as:

φ( τ ) =

) (

)

⎡ 2π ndetP ⎤1/2 t − 1⎦ ⎣( )

(15)

Here we need to derive the recursive expressions of y^t , Pθx and Pxx . Considering the complicacy in the nonlinearity of the system model, we use the first order Taylor expansion to approximate the model, namely:

( )

( ) ( ) = ∂f ( θ^ )/∂θ^ .

f ( θ ) ≈ f θ^t − 1 + f t′−T 1 θ − θ^t − 1 = f θ^t − 1 + f t′−T 1θ˜t where f t′− 1

t−1

(16)

t−1

Substituting (16) into (12)–(14), we obtain:

y^t,1 = x^t,1 ≈

Pθx,1 ≈

∫R

n

∫R

n

⎡ ^ ⎤ ^ ^ ′T ⎣⎢ f θt − 1 + f t − 1 τ − θt − 1 ⎥⎦φ( τ )dτ = f θt − 1

( )

(

)

( τ − θ^ )f ( τ − θ^ )φ( τ)dτ = f t−1

′T t−1

t−1

Pxx,2 ≈

∫R

⎡ ′T ⎤2 1 ^ ^ T f ″ τ − θ^ f τ − θ + τ − θ ⎢ t−1 t−1 t − 1 ⎥ φ( τ )dτ t−1 n ⎣ t−1 ⎦ 2

(

=f t′−T 1Pt − 1f t′− 1 +

T ⎡ 1 ⎤ exp⎢ − 2 τ − θ^t − 1 Pt−−11 τ − θ^t − 1 ⎥ ⎣ ⎦

(

(22)

2

( )

′ P t−1 t−1

(17)

(18)

) (

) (

)

1 tr f t″− 1Pt − 1f t″− 1 Pt − 1 2

(

)

(23)

where subscript 2 denotes the expression derived from the second order Taylor approximation. Theoretically, (21)–(23) with higher order Taylor approximation should make the algorithm (6)–(8) more robust on a bad θ^ . In 0

practical, however, that is difficult to achieve. First, the biased Pt obstructs the realization of this superiority of the second approximation algorithm. This deviation disturbs the estimates of yt,2 especially: as a simple example, if θ^ = θ while P ≠ 0, (21) cannot t

yield a reasonable estimate. On this condition, the formulas (21)– (23) may result in an unstable estimate even if θ^t is close to the true value sufficiently. Unfortunately, the true value of Pt is quite difficult to be obtained either in initialization or during identification process. In addition, the necessity of calculating f ″ in the second order approximation algorithm increases the computational complexity exceedingly, especially for θ with higher dimension. Therefore, it is considerable to modify the algorithm to be more robust on θ^ simply and effectively. 0

According to (9), (17) and (21), we have:

Pxx,1 ≈

∫R

⎡ ′T ⎤2 ^ ′T ′ n⎢ ⎣ f t − 1 τ − θt − 1 ⎥⎦ φ( τ )dτ = f t − 1Pt − 1f t − 1

(

)

( )

et,1 = yt − f θ^t − 1

(19)

where subscript 1 denotes the expression derived from the first order Taylor approximation. Remark 3. It is rational to use formulas (17)–(19) to implement algorithm (6)–(8) provided that the parameter estimates are sufficiently close to the true values. This is because the first order Taylor approximation is applicable only if the parameter estimation error is small enough; otherwise, the approximation is fallacious and can even lead to an unstable parameter estimates. Therefore, θ^0 close to the true value is needed for the algorithm. According to the Remark 3, a more efficacious approximation

(24)

1 et,2 = yt − f θ^t − 1 − tr f t″− 1Pt − 1 2

( )

(

)

(25)

Comparing (17)–(19) and (21)–(23), we have:

et,2 = λ t, eet,1

(26)

(

)

Pxx,2 + R2 = λ t, P Pxx,1 + R2

(27)

where λt,e and λt,P are coefficients; From (26) and (27) as well as (6) and (7), it is obtained:

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

θ^t = θ^t − 1 + et,2Pθx,2 Pxx,2 + R2

(

−1

= θ^t − 1 + λ t, eλ t−, P1et,1Pθx,1 Pxx,1 + R2

)

(

= θ^t − 1 + λ t−1et,1Pθx,1 Pxx,1 + R2

(

λt−1

−1

)

St−−41 ft′−T 1Pt ft′− 1 rt

−1

)

≤

where = Considering the more efficacious approximation effect of the second order Taylor expansion provided that θ^ is closed to the true value sufficiently, we should have 0 < λt , e ≤ 1 when excluding the interference of measurement noise. In addition, from (19) and (23), we have λt , P ≥ 1, then λt ≥ 1 follows. According to the analysis above and (28), we modify the algorithm by introducing a coefficient λ ≥ 1 into Kt:

Kt, m

=

)

(29)

where Kt,m is the modified update gain. By introducing a constant coefficient λ into the update gain, we modify the first order approximation algorithm to achieve a better tolerance to a bad θ^ similar to the second order one. As λ ≥ 1 is

′ St−−21 ft′−T 1PP t t − 1 ft − 1

f t′−T 1Pt − 1 ft′− 1 + St2− 1

St−−21 ft′−T 1Pt − 1Pt − 1 ft′− 1 f t′−T 1Pt − 1 ft′− 1 + St2− 1

(

) =S

St−−21tr f t′−T 1Pt − 1Pt − 1 ft′− 1 f t′−T 1Pt − 1 ft′− 1 + St2− 1

−2 t − 1tr

( Pt − 1 − Pt )

≤ λ−1R−2( trPt − 1 − trPt )

(37)

Therefore: ∞

∑

−1 = Pθx,1⎡⎣ λ Pxx,1 + R2 ⎤⎦

(

≤ St−−41 ft′−T 1Pt2 ft′− 1 =

(28)

λt , eλt−, P1.

5

δ 2St−−41 ft′−T 1Pt ft′− 1 rt

t=1

≤ δ 2λ−1R−2( trP0 − trP∞) ≤ δ 2λ−1R−2trP0

(38)

□

This completes the proof.

Assumption 5. There exist positive real numbers η1 and finite positive integer N so that for all t ≥ N , we have:

η2 and a

0

chosen from imprecise analysis, however, a strict mathematical analysis on the convergence of the algorithm, as well as the benefits of the improvement on the performance of the algorithm, will be given in the following section. To conclude the section, the recursive steps of the MEKF are summed up as below.

⎡ ⎤ θ^t = θ^t − 1 + Kt, met = θ^t − 1 + Kt, m⎣⎢ yt − f θ^t − 1 ⎥⎦

(30)

−1 Kt, m = − Pt − 1 ft′− 1⎡⎣ λ f t′−T 1Pt − 1 ft′− 1 + R2 ⎤⎦

(31)

( )

(

)

−1 Pt = Pt − 1 − Pt − 1 ft′− 1 ft′−T 1Pt − 1⎡⎣ λ f t′−T 1Pt − 1 ft′− 1 + R2 ⎤⎦

(

)

(32)

The modified algorithm degrades into the first order approximation algorithm if λ = 1.

t

f t′− 1 ft′−T 1 ≤ η2I

∑

η1I ≤

(39)

i=t−N+1

Lemma 2. Under the Assumptions 1–5, there exists a real number 2 2 S¯ > 0 so that St2− 1 ≤ S¯ . Proof. From Assumption 5 and (35), we have:

P0 ≥ P1 ≥ ⋯ ≥ Pt

(40)

f t′ f t′T ≤ η2I

(41)

Therefore:

St2− 1

(

(

)

≤ ( λ − 1)tr η2P0

4. Convergence analysis

2 + λR = S¯ 2

(42) □

This completes the proof. The following definition and three lemmas are required prior to the main conclusion. Define the variable St2− 1:

St2− 1 = ( λ − 1) ft′−T 1Pt − 1 ft′− 1 + λR2

)

≤ ( λ − 1) ft′−T 1P0 ft′− 1 + λR2 = ( λ − 1)tr P0 ft′− 1 ft′−T 1 + λR2

Lemma 3. When applying Algorithm (30)–(32) to System (4) on Assumptions 1–5, the following equality holds:

(33) ∞

rt − rt − 1 =∞ rt − 1 t=1

According to (32), we have:

(

∑ −1

)

Pt ft′− 1 = St2− 1Pt − 1 ft′− 1 f t′−T 1Pt − 1 ft′− 1 + St2− 1

(43)

(34) Proof. With Assumption 5 and (35), we obtain the following two inequalities:

Pt−1 = Pt−−11 + St−−21f t′− 1f t′−T 1

(35)

Lemma 1. When applying Algorithm (30)–(32) to System (4) on Assumptions 1–4, the following inequality holds:

∞

∑ t=1

δ 2St−−41 ft′−T 1Pt ft′− 1 rt

<∞

(36)

( )

where the real number rt is defined as rt = tr Pt−1 . Proof. From definition (33) and formula (32) and (34), we have:

⎛ lim t /rt = lim t /tr ⎜⎜ P0−1 + t →∞ t →∞ ⎝

t−1

⎞

i=0

⎠

∑ Si−2 fi′ fi′T ⎟⎟

t−1 ⎛ ⎞ −2 ≤ lim t /tr ⎜⎜ P0−1 + S¯ ∑ f i′ fi′T ⎟⎟ t →∞ ⎝ ⎠ i=0 k−1 ⎛ N−1 ⎡ ⎞⎤ −2 ′T ⎟⎥ = lim kN /⎢ tr P0−1 + S¯ ∑ tr ⎜⎜ ∑ f ′jN + i f jN + i⎟⎥ ⎢⎣ k →∞ ⎝ i=0 ⎠⎦ j=0

( )

−2

( )

≤ lim kN /kS¯ tr η1I < ∞ k →∞

(44)

and Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

∞

N

∞

rt − rt − 1 = rt t=1

∑

∑∑ k=0 i=1 ∞

∑

=

k=0

rkN + N − rkN = rkN + N

∑

tr

∑∑ k=0 i=1

(

∑

Assuming that

(

)

(

(

rkN + N

k →∞

) ≥ l S¯

(

∑

(

rkN + N

−2

1

( )

tr η1I > 0

)

−1 −2 ⎡ ′T ′ ⎤ × ⎢2(1 − μt − 1) − 1 − (1 − μt − 1) St−−21 ft − 1 Pt ft − 1 ⎥ , a . s. ⎣ ⎦

(46)

(

(47)

et,1 = xt + εt − x^t,1 = x˜ t,1 + εt ≈ f ′tT θ˜t − 1 + εt

(56)

− 1 ′T f t − 1θ˜t − 1

≥0

)

(48)

μt − 1 ≤

1 − St−−21 ft′−T 1Pt ft′− 1

(58)

⎛L ⎞ L ⎛ 1 St−−41δ 2f t′−T 1Pt f t′− 1 1⎞ E⎜⎜ t t − 1⎟⎟ ≤ t − 1 − L t − 1⎜ − ⎟+ , a. s . rt − 1 rt ⎠ rt ⎝ rt − 1 ⎝ rt ⎠

(60)

Considering the conclusion (36) of Lemma 1, and applying the martingale convergence Theorem [31] to (60), we deduce that Lt /rt convergences a.s. to a finite random variable and obtain:

⎛ 1 1⎞ − ⎟ < ∞ , a. s. rt ⎠ ⎝ rt − 1

∞

(50)

(59)

Substituting (33) and (34) into (59), (51) follows. From (58), we have:

∑ Lt − 1⎜ + εt

(57)

Solving (57), we have:

(49)

)

− 2 − 2 ′T St − 1 ft − 1Pt ft′− 1

)

holds, we have:

(

Therefore we have:

(

(

− 1 − 1 − μt − 1

E L t t − 1 ≤ L t − 1 + St−−41δ 2f t′−T 1Pt f t′− 1 , a. s .

where x˜ t is the system estimation error. As the Lagrange remainder term Ot  1 is neglected in (16), f t′T θ˜t − 1 is the approximation of the true model estimation error x˜ t . To amend this approximation, we introduce an unknown coefficient μt, as defined in the following formula, to take those residues into account [29,30].

μt − 1 = Ot − 1/x˜ t

−1

)

2 1 − μt − 1

) =∞

With (45) and (47), (43) follows. This completes the proof. □ Substituting (16) into (24), we obtain:

et = 1 − μt − 1

(55)

T Noticing that St−−21 ft′−T 1θ˜t − 1θ˜t − 1 ft′− 1 ≥ 0, if

−2 N−1 ′ S¯ ⋅tr ∑i = 0 f kN f ′T + i kN + i

k=1

⎤2 + εt ⎥ ⎦

−2 ′T ~ ⎤ ~T ′ ′T ′ ⎡ + St−−41 ft − 1 Pt ft − 1 ⎢(1 − μt − 1) ft − 1 θt − 1θt − 1 ft − 1 + δ 2⎥ , a . s. ⎣ ⎦ ~T ′ ′T ~ ′T ′ = L t − 1 + St−−41δ 2ft − 1 Pt ft − 1 − St−−21 ft − 1 θt − 1θt − 1 ft − 1

Therefore: ∞

− 1 ′T f t − 1θ˜t − 1

)

−1 ~T ′ ~T ′ ′T ~ ′T ~ E L t t − 1 = L t − 1 − 2(1 − μt − 1) St−−21θt − 1 ft − 1 ft − 1 θt − 1 + St−−21θt − 1 ft − 1 ft − 1 θt − 1

= l1, where l1 is a finite real number and

−2 N−1 ′ S¯ tr ∑i = 0 f kN f ′T + i kN + i

⎤ + εt ⎥ ⎦

Therefore:

(45)

considering (45), we have:

lim ( kN + N )

− 1 ′T f t − 1θ˜t − 1

)

⎡ +St−−41 ft′−T 1Pt ft′− 1⎢ 1 − μt − 1 ⎣

)

rkN + N t lim t →∞ rt

⎡ T −2St−−21θ˜t − 1 ft′− 1⎢ 1 − μt − 1 ⎣

rkN + N

k=0

(

k=0

T L t = L t − 1 + St−−21θ˜t − 1 ft′− 1 ft′−T 1θ˜t − 1

rkN + i − rkN + i − 1 rkN + N

N − 1 −2 ′ ′T ∑i = 0 SkN + i f kN + i f kN + i

−2 N−1 ′ S¯ tr ∑i = 0 f kN f ′T + i kN + i

∞

≥

∞

N

∞

rkN + i − rkN + i − 1 ≥ rkN + i

t=1

(61)

From (61), we have:

The main conclusion is summarized as below:

⎛ 1 1⎞ − ⎟= rt ⎠ ⎝ rt − 1

∞

Theorem 1. When applying algorithm (30)–(32) to System (4) on Assumptions 1–5, if.

∞

∑ Lt − 1⎜

∑

t=1

t=1

L t − 1 rt − rt − 1 ⋅ < ∞ , a. s. rt rt − 1

(62)

Considering the conclusion (43) of Lemma 3, we have Lt − 1/rt → 0, a.s., namely: −1

μt − 1 ≤

T P θ˜t − 1 t − 1 θ˜t − 1 → 0 , a. s. rt

( λ − 1)f t′−T 1Pt − 1 ft′− 1 + λR2 λf t′−T 1Pt − 1 ft′− 1 + λR2

(51)

It is obvious that:

holds, the algorithm can ensure that

θ˜t → 0 , a. s.

(52)

lim

Pt−−11

t →∞ rt − 1

>0

(64)

Therefore we obtain the conclusion (52). This completes the proof. □

Proof. Using the Lyapunov function.

T L t = θ˜t Pt−1θ˜t

(53)

Considering (30), (31), (34) and (35), the Lyapunov function becomes: T

(63)

T

L t = L t − 1 + St−−21θ˜t − 1 ft′− 1 ft′−T 1θ˜t − 1 − 2St−−21θ˜t − 1 ft′− 1et + St−−41 ft′−T 1Pt ft′− 1et2 Substituting (50) into (54), we have:

(54)

Remark 4. (51) can be considered as a sufficient condition of parameter convergence. This is the condition restricted on the value of μt, which reflects the proportional relationship between the Lagrange remainder term Ot  1 and the true model estimation error x˜ t . Obviously, a λ 4 1 expends this domain. Remark 5. (39) can be regarded as the persistent excitation condition of uniform convergence in parameter estimates since f ′t − 1 is related to the system input data. Actually, this condition is a transformation of the uniform observability condition in EKF algorithm [32,33] if we consider the parameter as a special kind of

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

state, namely

θt þ 1 ¼ θt.

7

(

et+ = yt − f Ut − 1, Θ^t

Remark 6. Although R2 is defined as an estimate of the noise variance δ2, the unbiased estimation of this value is not necessary. It is because that according to the proof of Theorem, no assumption on the unbiased estimation of δ2 is required. A “small” value of R2 is suggested in the MEKF.

)

(69)

θ^t = θ^t − 1 + Kt, met I⎡ e + − e ⎣ t

⎤

t ≤ l1⎦

(70)

−1 Pt = Pt − 1 − Pt − 1 ft′− 1 ft′−T 1Pt − 1⎡⎣ λ f t′−T 1Pt − 1 ft′− 1 + R2 ⎤⎦ I⎡ e + − e ≤ l ⎤ t 1⎦ ⎣ t

(71)

λ t = λ t − 1 + l2I⎡ e + − e

(72)

(

5. Some related problems 5.1. The selection of

λ

⎣ t

The most important modification of the algorithm is the introduction of parameter λ. The main function of λ is to enhance the convergence of the proposed algorithm so that to tolerate a worse initial condition. From Remark 4, it is obtained that (51) define the convergence domain of the algorithm. For the traditional EKF algorithm, named λ ¼ 1, this domain is enslaved to the relationship between f t′−T 1Pt − 1 ft′− 1 and R2. If f t′−T 1Pt − 1 f t′− 1 ≫ R2, the convergence domain defined by (51) becomes quite small so that the algorithm maybe unstable even if the initial parameter estimates are close to the true values. For the MEKF, λ plays an important role for the extension of the convergence domain: if f t′−T 1Pt − 1 f t′− 1 ≫ R2, we have

μt − 1 ≤

1 − λ−1 approximately; if

f t′−T 1Pt − 1 f t′− 1 ≪ R2, we have

μt − 1 ≤ 1 approximately. That means the convergence domain of the algorithm is controlled by λ. On the basis of the above analysis, if λ is greater than 1, the convergence domain will be extended greatly. Therefore, a large λ is suggested. However, the stability of the algorithm can’t be guaranteed in any case, especially if the prior knowledge is absent. How to set an appropriate value of λ is a significant issue to be discussed. As the true values of the system parameters are unknown, it is difficult to give an appropriate value of λ before the algorithm is running. Adaptive adjusting λ during the identification process is a feasible method. Here we suggest a practical method for adjusting λ. During the identification process, it is impossible to check whether the parameter estimates are close to the true values directly. But it may be indirectly reflected from the system estimation error: if it becomes larger, it is believed that the parameter estimating is carrying out in the wrong direction. Based on this, a further modified algorithm is given. Firstly, define the model posterior estimation error et+ as the following:

(

^ et+ = yt − f Ut − 1, θ¯t

)

(65)

et+

If − et ≥ l1, where l1 ≥ 0, it is believed that the identification diverges as the convergence domain is not large enough. Therefore, the updating of Pt and θ^t is invalid at this sampling time and a larger value of λ is given for the next identification process. Based on this updating mechanism, the recursive steps of the MEKF are as follows:

(

)

^ et = yt − f Ut − 1, θ¯t − 1

(66)

−1 Kt, m = − Pt − 1 ft′− 1⎡⎣ λ f t′−T 1Pt − 1 ft′− 1 + R2 ⎤⎦

(

)

(67)

)

⎤

t > l1⎦

where λ 0 ≥ 1; l2 > 0; I[ ⋅] is the indicator function, i.e., IA( ω) = 1 if ω ∈ A ; IA( ω) = 0, otherwise, Θ^t is the parameter interval estimation vector in the recursive steps. By introducing the mechanism of adaptive adjusting λ, the adaptive enlarging of the convergence domain is realized. 5.2. For the unstable system In the Assumption 3, the stable linear block is required. For the unstable system, the algorithm can also work with some modification. The modification focuses on the parameterized representation of the system. If the linear block is unstable, the FIR filter cannot be used as the linear part model. Here, we use the following IIR filter model as a substitute:

wt + d1wt − 1 + ⋯ + dmwt − n = a1vt − 1 + ⋯ + a nvt − m

(73)

Notice that the internal variable wt and vt is immeasurable, the parameterized representation of the nonlinear part are also required. The representation of the input nonlinear part is with the same way as formula (1). The focus is on the output nonlinear part. A feasible method is to use the measurement of system output yt to express variable wt by the following formula [11,19,20]: q

ϕi ( yt ) ( ) ∑ ⌢ci⌢

wt = fO−1 yt =

(74)

i=1

⌢ where, ϕi and ⌢ ci are the new basis function and parameter, respectively. It is obviously that the irreversible assumption for the output nonlinear part is necessary in this method. In order to eliminate this deficiency, another method for expressing wt is as follows [23]:

wt = AVt − 1 − DWt − 1

(75) T

T

where Wt − 1 = ⎡⎣ wt − 1, wt − 2, ⋯ , wt − n⎤⎦ , D = ⎡⎣ d1, d2, ⋯ , dm⎤⎦ . The representation of the output nonlinear part is with the same way as formula (2). When the algorithm is running, the estimate of wt is ^ ^ obtained by At − 1, Dt − 1 and the previous estimates of vt and wt, namely:

^ ^ ^ = A^ V^ − D w t t−1 t−1 t − 1Wt − 1

(76)

Based on the two parameterized method, we can easily acquire the whole system model. The MEKF can also work well on the estimation of the unknown parameters in the two system models. For the parameterized representation (2), it is noticed that the initial estimate error of parameter set di should not be too large even if a large λ is selected. 5.3. The robustness properties of the algorithm

Θ^t = θ^t − 1 + Kt, met

(68)

The robustness is an important property for an identification

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

algorithm. The purpose of introducing λ is to enhance the robustness to the initial parameter estimation error. In Section 3, we assume that the system model structure is known and the disturbance sequence is a stochastic noise with zero mean and bounded variance. These assumptions ensure the uniform convergence of the parameter estimates. Actually, the introduction of λ also enhances the robustness of the algorithm to the unmodeled dynamic and unmodeled disturbance. Generally speaking, the effect of the traditional KF will be worse if the model structure is mismatched or the system subjects to a nonzero mean noise with no priori knowledge. Therefore, many new algorithms are proposed to overcome these shortcomings [34–37], i.e., H1 filter. Compare with KF, H1 filter gives a larger Pt and Kt. A larger Kt means the algorithm will be more focused on the system measured values so that it will give stronger robustness on the unmodeled dynamic and unmodeled disturbance [38]. The introduction of λ has similar functions. From (32), it is obviously that the introduction of λ gives a larger Pt than traditional algorithm. It seems that we obtain a smaller Kt from (32). However, with the difference between Pts obtained by the two algorithms increasing gradually, the Kt given by the MEKF will quickly become larger than the traditional one. Therefore, the MEKF is more robustness to EKF. For this, it will be further illustrated in the simulation examples. However, as the Assumptions 1 and 4 are not met, the algorithm will only give a bounded parameter estimates if the condition (51) holds. For identification and tracking the time-varying parameters in the time-varying H-W system, further modification for the algorithms is required. Here, we give a shallow discussion on this issue. In order to achieve the identification and tracking the timevarying parameters in the system, we add a term Q into the recursive steps (31) and (32), namely:

Kt, m = − ( Pt − 1 + Q )f t′− 1 λ⎡⎣ f t′−T 1( Pt − 1 + Q )f t′− 1 + R2⎤⎦

{

−1

}

(77)

Pt = Pt − 1 + Q − ( Pt − 1 + Q )f t′− 1f t′−T 1( Pt − 1 + Q ) × λ⎡⎣ f t′−T 1( Pt − 1 + Q ) ft′− 1 + R2⎤⎦

{

−1

}

(78)

Bounded parameter tracking results will be obtained by the above recursive steps.

6. Numerical simulations

w2, t = 0.14v2, t − 1 + 0.23v2, t − 2 + 0.28v2, t − 3 + 0.25v2, t − 4 + 0.14v2, t − 5 − 0.1v2, t − 6 − 0.12v2, t − 7 − 0.1v2, t − 8 − 0.05v2, t − 9 + 0.04v2, t − 10 + 0.07v2, t − 11 + 0.05v2, t − 12 + 0.02v2, t − 13 while the output blocks are also represented as polynomial functions:

x1, t = w1, t + 1.2w1,2 t , x2, t = w2, t − 0.5w2,2 t The example contains two simulation parts. In part one, we verify the convergence of the algorithm and the effect of λ on the extension of the convergence domain. In part two, we show the robustness of the algorithm on the unmodeled dynamic and unmodeled disturbance. Part one. In this simulation part, the noise is generated by the zero mean WGN and its signal-to-noise ratio (SNR) is 20 dB. The subsystem inputs are both WGN sequences with mean zero and variance 1. The identification model is the same as the system model. P0 ¼108I. In order to test the performance of the algorithm, an ill initialization condition, namely all the initial parameter estimates are far from the true values, is considered as θ^ = 50 × (1, …, 1)T . 0

Three kinds of algorithms are used for comparison in this simulation: 1. The traditional EKF, namely λ ¼1; 2. The MEKF with a fixed λ, set λ ¼5; 3. The MEKF with a variational λ, set λt = λt − 1 + l2I⎡ ⎣ λ0 ¼5; l1 ¼1; l2 ¼0.5.

et+ − et > l1⎤⎦ ,

where

In order to verify the stability of the algorithm, we carry out 100 times independent simulation tests. In each test, the simulation goes through 10,000 samples. The norms of the parameter estimation error ||θ˜t||2 at the end of simulation in each simulation experiments are shown in Fig. 3. The final values of λt of algorithm 3 in each experiment are shown in Fig. 4. All the identification processes of the algorithm 1 are diverging, because the convergence domain of traditional EKF is not sufficiently large. Therefore, Fig. 3 doesn’t show these results. From Fig. 3, it can be seen that algorithm 2 gives divergent results only in 2 tests and algorithm 3 obtains satisfactory results in all tests. In Fig. 4, it is noticed that the mechanism of adaptive adjusting λt plays a role in some experiments. Therefore, the algorithm 3 gives a better result than algorithm 2. The simulation results show the function of λ in extension of the algorithm convergence domain and the practicability of the mechanism of adaptive adjusting λt.

In this section, we use two simulation examples to show the validity of the proposed algorithm. Example 1. In this simulation example, we consider a MISO H-W system composed of two subsystems with the input nonlinear blocks represented as polynomial functions:

v1, t = u1, t + 0.4u1,2 t − 0.2u1,3 t , v2, t = u2, t + 0.6u2,2 t − 0.6u2,3 t and the linear blocks represented as: w1, t = 0.15v1, t − 1 + 0.2v1, t − 2 + 0.28v1, t − 3 + 0.22v1, t − 4 + 0.13v1, t − 5 + 0.1v1, t − 6 + 0.08v1, t − 7 − 0.08v1, t − 8 − 0.09v1, t − 9 − 0.11v1, t − 10 − 0.07v1, t − 11

, Fig. 3. The norms of the parameter estimation error.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 4. The final values of λt of algorithm 3.

9

Fig. 7. The average system estimation error.

Part two. After demonstrating the robustness of the MEKF on the initial parameter estimation error, we verify the robustness on the unmodeled dynamic and unmodeled disturbance in this simulation part. The noise term is generated by the GN sequences with variance 0.05 and the mean value μ of GN is generated by μ = 0.1 sin( 0.01t ). Reduced-order models for the two linear parts are used in the simulation. We set m1 ¼ 10 and m2 ¼ 12 for the identification model rather than m1 ¼ 11 and m2 ¼ 13 in the original system model. In order to facilitate the study of the robustness on the unmodeled dynamic and unmodeled disturbance of the algorithm, an initial parameter estimate closed to the “true value” is used in this simulation part. We set θ^ = 0.01 × (1, … , 1)T . 0

Two kinds of algorithms are used for comparison in the simulation: Fig. 5. The norms of the parameter estimation error.

1. The traditional EKF, namely λ ¼ 1; 2. The MEKF with a fixed value of λ, set

λ ¼ 5;

M The index about the system estimation error x¯˜t = ∑i = 1 x˜ t , i /M at each sampling point is used to show the identification results in Fig. 7, where M is the number of the tests. The average norms of M the update gain ∑i = 1 ‖Kt , i‖2 /M of 100 tests at each sampling point are shown in Fig. 8. There are two divergent results obtained from the algorithm 1. Therefore, the results shown in Fig. 7 are calculated by the rest 98 results. It is obviously that the MEKF shows stronger robustness on the unmodeled dynamic and unmodeled disturbance than the traditional one as its x¯˜t is much smaller. Fig. 8 verifies the discussion in 5.3, namely the Kt given by the proposed algorithm will

Fig. 6. The values of λt at the end of simulation of algorithm 3.

In order to further illustrate the effectiveness of the algorithm, a much worse initialization condition θ^0 = 100 × (1, … , 1)T is used in the simulation. The results are shown in Figs. 5 and 6. From Fig. 5, it can be seen that more divergent results are obtained from algorithm 2 while algorithm 3 still gives convergent results in all tests. This is because the convergence domain is no longer large enough so that the mechanism of adaptive adjusting λt plays a greater role in the simulation processes. This point can be seen from Fig. 6, where more values of λt in the experiment are changed by the adaptive adjusting mechanism. These results show the practicability of the mechanism of adaptive adjusting λt further.

Fig. 8. The average norms of the update gain.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

⎡ ⎛ ⎞2/7 ⎤ P 5⎢ ⎜ o ⎟ − 1⎥ ⎢ ⎝ Ps ⎠ ⎥ ⎣ ⎦

(81)

Ps = − 0.5622w1, t + 0.1095w1,2 t

(82)

Po = 0.7021w1, t − 0.0974w1,2 t

(83)

Ma =

Fig. 9. Wind tunnel system.

quickly become larger than the traditional one with the difference between Pts obtained by the two algorithms increasing gradually. This is the reason why the proposed algorithm has stronger robustness. Example 2. In this example, we use a practical wind tunnel system simulation case to show the practicability of the proposed algorithm. The wind tunnel system is used primarily for the purpose of flight sciences. One of the key concepts behind wind tunnels is the Mach number (MA). Having an accurate grasp of the dynamic properties of the wind tunnel system on-line will have direct manipulation over the MA.

where Ps and Po are the static pressure and the total pressure in the stilling chamber room, respectively. Before the system identification, we should parameterize the system model. The following expressions are used to parameterize the backlash nonlinearity in the choke finger system and the main exhaust system:

vt = ( ut − c )mh1(ut ) + ( ut + c )mh2(ut ) + vt − 1⎡⎣ 1 − h1(ut )⎤⎦⎡⎣ 1 − h2(ut )⎤⎦

(84)

where The wind tunnel system we studied is in the China Aerodynamics Research and Development Center. The system block diagram is shown in Fig. 9 [39], where the control signal u1,t and u2,t drive the choke finger system and the main exhaust system to adjust the system opening w1,t and w1,t, respectively. These two variables manipulate the system output xt, namely MA. yt is the measurements of xt. Both of the choke finger system and the main exhaust system are operated by the hydraulic systems with self regulating function. Generally speaking, the inputs and outputs characteristic of this kind of hydraulic systems can be represented by a linear dynamic model. However, due to mechanic friction and other interference, the overall hydraulic system experiences a backlash nonlinear characteristic. Therefore, we can use an H model with a backlash nonlinearity to express the choke finger system and the main exhaust system. This kind of nonlinearity can be formulated by:

⎧ ( ut − c )m ut ≥ z r ⎪ vt = ⎨ ( ut + c )m ut ≤ zl ⎪ vt − 1 zl < ut < z r ⎩

(79)

where zr = vt − 1/m + c ; zl = vt − 1/m − c . For the choke finger system, m1 ¼ 1, c1 ¼ 0.3; For the main exhaust system, m2 ¼ 1, c2 ¼ 0.2. The linear part of the H model can be expressed by the following FIR filter model:

wt = AT Vt − 1

(80)

⎧ 1 ( ut − c )m − vt − 1 > 0 h1(ut ) = ⎨ ⎩ 0 ( ut − c )m − vt − 1 ≤ 0

(85)

⎧ 1 ( ut + c )m − vt − 1 < 0 h2(ut ) = ⎨ ⎩ 0 ( ut + c )m − vt − 1 ≥ 0

(86)

⎪

⎪

⎪

⎪

And the formula (80) is used as the model for the linear part of the choke finger system and the main exhaust system. The original system orders are used in the parameterized model. It is very difficult to obtain the model structure of the flow field. Even the formulas (81) to (83) in the simulation are approximate. Therefore, the following polynomial function is used to approximate the static characteristic of the flow field: 3

xt = f1, O (w1, t ) + f2, O (w2, t ) =

i=1

i=1

(87)

Summarizing the parameterizing process mentioned above, we obtain the system model structure shown in Fig. 10. Recursive identification of the parameters in this model is exactly what we're studying in the paper. The simulation comprises 5000 samples. In order to ensure the consistency of parameter estimation, the estimates of a11, a21, m1 and m2 are set to be 1. The noise is generated by the zero mean WGN and its SNR is 30db. The system inputs u1,t and u2,t are 214 þ Δ and 350 þ Δ, respectively, where Δ is generated by an uniform distribution sequence on [  5,5]. P and θ^ are initialized as 108I and (0.01, … , 0.01)T , re0

For the choke finger system, A1 ¼ [0.3668, 0.1899, 0.1681, 0.1539, 0.0632, 0.0553]T; For the main exhaust system, A2 ¼ [0.5426, 0.2476, 0.1120, 0.0693, 0.0187, 0.0135, 0.0039]T. The characteristic of the flow field can be represented by the following formulas [40,41]:

3

∑ c1, iw1,i t + ∑ c2, iw2,i t

0

spectively. R2 is set to be 0.01. An adaptive adjusting λt is used and we set λt = λt − 1 + I⎡ e+ − et > 0.05⎤, where λ0 ¼ 10; 100 independent ⎣ t

⎦

simulation tests are carried out. In each experiment, the simulation goes through 5000 samples. The average estimation error of each parameter and system output at each sampling time are

Fig. 10. Wind tunnel system model block diagram.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 11. Estimates of A1.

Fig. 12. Estimates of A2.

11

Fig. 14. Estimate error of Ma.

the study is extended to a rarely mentioned multiple input single output Hammerstein-Wiener systems. Two basic recursive algorithms are proposed based on the extended Kalman filter algorithm and the Taylor approximation of the nonlinear system model. By introducing a coefficient λ into the first order approximation algorithm, a modified algorithm is proposed, which copes with the problem of small parameter convergence domain of the first order one and the application limit of the second order one. We give a convergence analysis on the modified algorithm and obtain a sufficient condition for achieving the uniform convergence of the parameter estimation. Based on this condition, the validity of the modification on the expansion of convergence domain is verified. Furthermore, we suggest a practical method for adjusting λ, which can adaptive enlarge the convergence domain during the identification process. We also discuss the application of the algorithm to the unstable system and the robustness on the unmodeled dynamic and unmodeled disturbance. Simulation examples show satisfactory performance of the modified algorithm. For further research it is of interest to consider identification multiple input multiple output Hammerstein-Wiener system, unstable Hammerstein-Wiener system or time-varying Hammerstein-Wiener system. The algorithm which is worth considering is to estimate the parameters of the system linear part using parameter identification method and to estimate the system nonlinear input and output relations using the non-parameter identification algorithm.

References

Fig. 13. Estimates of C.

shown in Figs. 11–14. Figs. 11–13 show that all the parameter estimates converge to their true values in each simulation test. Fig. 14 shows that the estimate error of Ma in each test converges to zero. The results of the practical system simulation further illustrate the validity of the MEKF.

7. Conclusions In this paper, we present a study on the recursive identification of single input single output Hammerstein-Wiener systems, and

[1] Wang DQ. Hierarchical parameter estimation for a class of MIMO Hammerstein systems based on the reframed models. Appl Math Lett 2016;57(7):13–9. [2] Jafari M, Salimifard M, Dehghani M. Identification of multivariable nonlinear systems in the presence of colored noises using iterative hierarchical least squares algorithm. ISA Trans 2014;53(4):1243–52. [3] Mehta U, Majhi S. Identification of a class of Wiener and Hammerstein-type nonlinear processes with monotonic static gains. ISA Trans 2010;49(4):501–9. [4] Ding F, Liu XG, Chu J. Gradient-based and least-squares-based iterative algorithms for Hammerstein systems using the hierarchical identification principle. IET Control Theory Appl 2013;7(2):176–84. [5] Vörös J. Parameter identification of Wiener systems with multisegment piecewise-linear nonlinearities. Syst Control Lett 2007;56(2):99–105. [6] Vörös J. Compound operator decomposition and its application to hammerstein and wiener systems. In: Giri F, Bai E-W, editors. Block-oriented nonlinear system identification. London: Springer; 2010. [7] Gu Y, Ding F. Parameter estimation for an input nonlinear state space system with time delay. J Frankl Inst 2014;351(12):5326–39. [8] Ikhouane F, Giri F. A unified approach for the parametric identification of SISO/ MIMO Wiener and Hammerstein systems. J Frankl Inst 2014;351(12):1717–27. [9] Zhang B, Mao ZZ Adaptive control of stochastic Hammerstein systems with dead-zone input non-linearity. Transactions of the Institute of Measurement and Control DOI: http://dx.doi.org/10.1177/0142331214546521(2014).

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i

12

F. Yu et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[10] Ławryńczuk M. Nonlinear predictive control for Hammerstein–Wiener systems. ISA Trans 2015;55(3):49–62. [11] Bai EW. An optimal two-stage identification algorithm for Hammerstein-Wiener nonlinear systems. Automatica 1998;34(3):333–8. [12] Zhu YC. Estimation of an N-L-N Hammerstein-Wiener model. Automatica 2002;38(9):1607–14. [13] Crama P, Schoukens J. Hammerstein-Wiener system estimator initialization. Automatica 2004;40(9):1543–50. [14] Wills A, Ninness B. Generalised Hammerstein–Wiener system estimation and a benchmark application. Control Eng Pract 2012;20(11):1097–108. [15] Wills A, et al. Identification of Hammerstein-Wiener models. Automatica 2013;49(1):70–81. [16] Vörös J. Iterative identification of nonlinear dynamic systems with output backlash using three-block cascade models. Nonlinear Dyn 2015;79(3):2187–95. [17] Vörös J. Identification of nonlinear dynamic systems with input saturation and output backlash using three-block cascade models. J Frankl Inst 2014;351 (12):5455–66. [18] Bolkvadze GR. The Hammerstein-wiener model for identification of stochastic systems. Autom Remote Control 2003;64(9):1418–31. [19] Wang DQ, Ding F. Extended stochastic gradient identification algorithms for Hammerstein-Wiener ARMAX systems. Comput Math Appl 2008;56 (12):3157–64. [20] Bai EW. A blind approach to the Hammerstein-Wiener model identification. Automatica 2002;38(6):967–79. [21] Wang DQ, Ding F. Hierarchical least squares estimation algorithm for Hammerstein-Wiener systems. IEEE Signal Process Lett 2012;19(12):825–8. [22] Yu F, Mao ZZ, Jia M. Recursive identification for Hammerstein–Wiener systems with dead-zone input nonlinearity. J Process Control 2013;23(8):1108–15. [23] Yu F, Mao ZZ, Jia M, Yuan P. Recursive parameter identification of Hammerstein–Wiener systems with measurement noise. Signal Process 2014;105 (12):137–47. [24] Zhang B, Hong HC, Mao ZZ Adaptive control of Hammerstein–Wiener nonlinear systems. International Journal of Systems Science DOI: http://dx.doi. org/10.1080/00207721.2014.971089(2014). [25] Cai ZJ, Bai EW. How nonlinear parametric Wiener system identification is under gaussian inputs? IEEE Trans Autom Control 2012;57(3):738–42. [26] Chen HF. Recursive identification for Wiener model with discontinuous piece-

wise linear function. IEEE Trans Autom Control 2006;51(3):390–400. [27] Li GQ, Wen CY. Identification of Wiener systems with clipped observations. IEEE Trans Signal Process 2012;60(7):3845–52. [28] Wang Z, Wang Y, Ji Z Stochastic gradient algorithm for multi-input multioutput Hammerstein FIR-MA-like systems using the data filtering. Journal of the Franklin Institute DOI: http://dx.doi.org/10.1016/j.jfranklin.2015.01.015 (2015). [29] Boutayeb M, Aubry D. A strong tracking extended Kalman observer for nonlinear discrete-time systems. IEEE Trans Autom Control 1999;44(8):1550–6. [30] Boutayeb M, Rafaralahy H, Darouach M. Convergence analysis of the extended Kalman filter used as an observer for nonlinear deterministic discrete-time systems. IEEE Trans Autom Control 1997;42(4):581–6. [31] Goodwin GC, Sin KS. Adaptive filtering prediction and control.Mineola, New York: Courier Dover Publications; 2013. [32] Reif K, Gunther S, Yaz E, Unbehauen R. Stochastic stability of the discrete-time extended Kalman filter. IEEE Trans Autom Control 1999;44(4):714–28. [33] Song Y, Grizzle JW. The extended Kalman filter as a local asymptotic observer for discrete-time nonlinear systems. J Math Syst, Estim Control 1995;5(1):59– 78. [34] Zhang D, Cai W, Xie L, et al. Nonfragile distributed filtering for T–S fuzzy systems in sensor networks. IEEE Trans Fuzzy Syst 2015;23(5):1883–90. [35] Zhang D, Shi P, Zhang WA, et al. Energy-efficient distributed filtering in sensor networks: a unified switched system approach. IEEE Trans Cybern 2016:1–12. [36] Zhang D, Song H, Yu L. Robust fuzzy-model-based filtering for nonlinear cyber-physical systems with multiple stochastic incomplete measurements, 2016, 1–13. [37] Zhang D, Cai W, Wang QG. Energy-efficient H1 filtering for networked systems with stochastic signal transmissions. Signal Process 2014;101(8):134–41. [38] Simon D. Optimal state estimation: Kalman, H1, and nonlinear approaches. Wiley-Interscience; 2006. [39] Wang X, Li S, Cai W, et al. Multi-model direct adaptive decoupling control with application to the wind tunnel system. ISA Trans 2005;44(1):131–43. [40] Soeterboek RAM, Pels AF, Verbruggen HB, et al. A predictive controller for the Mach number in a transonic wind tunnel. IEEE Control Syst 1991;11(1):63–72. [41] Nott CR, Ölçmen SM, Lewis DR, et al. Supersonic, variable-throat, blow-down wind tunnel control using genetic algorithms, neural networks, and gain scheduled PID. Appl Intell 2008;29(1):79–89.

Please cite this article as: Yu F, et al. Recursive parameter estimation for Hammerstein-Wiener systems using modified EKF algorithm. ISA Transactions (2017), http://dx.doi.org/10.1016/j.isatra.2017.05.012i