Nonlinear systems widely exist in practical applications, like communication systems, chemical processes, biomedical systems and so on. Therefore, nonlinear systems identificati Jingfan Liu on is quite significant both in theory and application. This thesis presents the identification algorithms for a class of nonlinear systems based on the Youth Project of Central University. Considering the identification of the input nonlinear systems with the colored noise, An extended Newton recursive algorithm are derived for comparison. In the simulation, the results show that the Newton recursive algorithm can get better accurate parameter estimates, The simulation results show the effectiveness of the proposed algorithms.
linear and nonlinear subsystems of the hysteresis Hammerstein system respectively using least squares [23]. Xiang Wei and Zonghai Chen proposed a new identification method for Volterra sequence (Laguerre function), a specific dynamic model of the nonlinear system. Many new methods are also used to identify nonlinear systems such as neural networks and genetic algorithms [24,25]. 2
Hammerstein model nonlinear system is the input nonlinear system, it consists of a static nonlinear segment and a dynamic linear segment(Figure 1) [26], where y(t) is output, u(t) is input, u(t) is the output of the nonlinear part and the noise v(t) is assumed to be i.i.d. random sequences with zero mean, A(z) , B(z) and D(z) are polynomials in the unit backward shift operator z−1[z−1 y(t) = y(t − 1)] , with A(z) =1+α1z−1 +α2z−2 ++α n z−na ,
a B(z) = β1z−1 + β 2 z−2 + β 3z−3 + + β nb z−nb , D(z) = 1 + d1z−1 + d2 z−2 + + dn z−nd .
d
The intermediate variables μ (t) , x(t) and h(t) are immesurable, and g(⋅) is static nonlinear function of state. The nonlinear part is an unknown polynomial that can be expressed as [27] : = g(μ (t))c We can write the Hammerstein-CARMA model (Figure 1) into the following formula: μ (t) = g(μ (t)) = c1g1(μ (t)) + c2 g2 (μ (t)) + + cnc gnc (μ (t))
A(z) y(t) = B(z)μ (t) + D(z)v(t) (1) (2)
The unknown parameters needed to be estimated are: the linear subsystem parameters α , β , d and the nonlinear part parameters c .Let the superscript T represent the matrix transpose. α := [α1,α 2 ,,α n ]T ∈ na ,
a d := [d1, d2 ,, dn ]T ∈ nd，
d α γ := ∈ na +nd , d γ θ := β ∈ n , n := na + nb + nc + nd .
c β := [β1,β 2 ,,β n ]T ∈ nb ,
b cˆt γˆt γ ˆ at time t .Define the information vector φ(t) as: Let θˆ(t)：= βt denote the estimate of θ := β c The output sequence can be written as: φ(t) := [− y(t −1), − y(t − 2),, − y(t − na ), v(t −1), v(t − 2),, v(t − nd )]T ∈ na +nd
y(t) = ζ T (t)J + v(t) Thereinto J is the parameterized vector, ζ (t) is the information vector, and they are defined as: α J := ∈ na +nbnc +nd ,
β ⊗ c ζ (t) := [φT (t), g(μ (t −1)), g(μ (t − 2)),
g(μ (t − nb ))]T ∈ na +nbnc +nd .
β ⊗ c is the Kronecker product of β and c . ln most of the existing papers, the combined parameter β ⊗ c is identified, and the combined parameter needs to be decomposed after the identification result is obtained [28], which increases the computational burden. The goal of this paper is to identify the parameters through the extended Newton recursive algorithm, obtain the parameter vectors α , β , c and d . 3 3.1
In this section, Newton method will be used to derive the augmented Newton recursive identification algorithm based on Hammerstein-CARMA model, its basic idea is to introduce stacking output vector and stacking information matrix. Define the input information:
We define a quadratic criterion function as follows: g(μ (t −1)) G(t) := g(μ (t − 2) ∈ nb×nc g(μ (t − nb )) y(t) = φT (t)γ + β TG(t)c + v(t)
J1(θ) = J1(γ, β, c) = [ y(t) − φT (t)γ − β TG(t)c]2
Since the Hessian matrix H[J1(γ, β, c)] of the criterion function J1 is singular, it is useful to introduce the stacked data, so Newton algorithm is used to solve the identification optimization problem.
φ(t)φT (t) φ(t)cTGT (t) φ(t) βTG(t) H[J1(γ, β, c)]=2 G(t)cφT (t) G(t)cφT (t)GT (t) h23(γ, β, c,t) ∈ n×n
GT (t) βφT (t) h32 (γ, β, c,t) GT (t) ββTG(t)
Where: h23 (γ, β, c, t)：= − ∂ ∂c = − G(t)[ y(t) − φT (t)γ + βTG(t)c] {G(t)c[ y(t) − φT (t)γ + βTG(t)c]} h32 (γ, β, c,t)：= − ∂ {GT (t) β[ y(t) − φT (t)γ + βTG(t)c]}
∂b = − GT (t)[ y(t) − φT (t)γ + βTG(t)c]+GT (t) βcTGT (t) +G(t)cβTG(t) ∈ nb ×nc
=h23T (γ, β, c, t) ∈ nc×nb (3) (4) (5) (6) (7)
Consider the newest p data and define stacked output vector Y ( p,t) and stacked matrices Φ0 ( p, t) , Φ(c, t) and Ψ ( β, t) .
y(t) Y ( p, t)：= y(t −1) ∈ p y(t − p + 1) φT (t) Φ0 ( p, t)：= φT (t-1) ∈ p×(na +nd ) φT (t − p + 1) Φ(c, t)：= cTG T (t)
cTG T (t-1)
∈ p×nb cTG T (t − p + 1)
β TG(t) Ψ ( β, t)：= β TG(t-1) ∈ p×nc β TG(t − p + 1) ：= Y ( p, t) − Φ0 ( p, t)γ −Ψ ( β, t)c 2
= Y ( p, t) − Φ0 ( p, t)γ − Φ( β, t)γ 2 Then define a new criterion function:
J2 (θ) = J2 (γ, β, c) (8)
Eq.(8) is equivalent to the following criterion function constructed from the data in a dynamical window with length p .
t J3 (γ, β, c) = [ y(i) − φT (i)γ − β TG(i)c]2 (9)
i=t− p+1
That is J2 (θ) = J3 (γ, β, c) . If we take t = N and p = N ( N is the data length), then Eq.(8) and Eq.(9) are the least squares criterion functions [29].
Computing the gradient of J2 (γ, β, c) gives: Φ0T ( p,t) gradθ [J2 (γ, β, c)] = −2 ΦT (c, t) [Y ( p,t) − Φ0 ( p,t)γ −Ψ ( β,t)c] Ψ T ( β,t) Φ0T ( p,t) = − 2 ΦT (c, t) [Y ( p,t) − Φ0 ( p,t)γ − Φ(c,t) β]
Ψ T ( β,t)
Define the extended generalized information matrix Ξ (t) and expanding innovation into innovation vector E( p, t) as: Φ0T ( p, t) Ξ (t)：= ΦT (cˆt−1, t) ∈ n× p Ψ T ( βˆt−1, t)
= Y ( p,t) −Φ0 ( p,t)γˆt−1 −Φ(cˆt−1,t)βˆt−1 ∈p Thus, we have: Φ0T( p,t) gradθ [J2(γ, β,c)] = −2ΦT(cˆt−1,t) Y( p,t) −Φ0( p,t)γˆt−1 −Ψ(βˆt−1,t)cˆt−1 ΨT(βˆt−1,t) Φ0T( p,t) = − 2ΦT(cˆt−1,t) Y(p,t) −Φ0(p,t)γˆt−1 −Φ(cˆt−1,t)βˆt−1
ΨT(βˆt−1,t) = − 2Ξ(t)E( p,t) Computing the Hessian matrix of the criterion function J2 (γ, β, c) .
H[J2(γ, β,c)]= ∂gardθ [J2(γ, β,c)]
∂θ T Φ0T ( p,t)Φ0( p,t) Φ0T ( p,t)Φ(c,t) Φ0T ( p,t)Ψ (β,t) = 2 ΦT (c,t)Φ0( p,t) ΦT (c,t)Φ(c,t) H23(γ, β,c,t) ∈ n×n
Ψ T (β,t)Φ0( p,t) H23T (γ, β,c,t) Ψ T (β,t)Ψ (β,t) H23(γ, β,c,t)：= − ∂∂c {ΦT (c,t)[Y ( p,t) −Φ0( p,t)γ −Ψ (β,t)c]}
= − ∂∂c {[G(t)c,G(t-1)c，,G(t-p +1)c][Y ( p,t) −Φ0( p,t)γ −Ψ (β,t)c]}
We can summarize the Newton extended recursive algorithm (the H-ENR algorithm) for the Hammerstein-CARMA models as follows: θ (t) = θ (t −1) + Πˆ −1(t)Ξ (t) E ( p, t) (13) v(t) = y(t) − φT (t)γ − β TG(t)c v(t) = y(t) − φT (t)γ − β TG(t)c = − ∂∂c ip=−01G(t − i)c y(t − i) − φT (t − i)γ − βTG(t − i)c
p−1 = p−−1i=0 {G(t − i) y(t − i) − φT (t − i)γ − βTG(t − i)c −G(t − i)cβTG(t − i)} ={G(t − i) −y(t − i) + φT (t − i)γ + βTG(t − i)c + G(t − i)cβTG(t − i)} i=0 p−1 ={G(t − i) −y(t − i) + φT (t − i)γ + βTG(t − i)c} + ΦT(c,t)Ψ (β,t)∈nb×nc i=0
Using the Newton method to minimize J2 (θ) , we can obtain the following recursive relation of computing θ (t) :
= θ(t −1) + 2{H[J2 (ηˆt−1, bˆt−1,cˆt−1)]} θ(t) = θ(t −1) −{H[J2(ηˆt−1,bˆt−1,cˆt−1)]}−1 gradθ[J2(ηt−1, bt−1,ct−1)] −1 Ξ(t)E( p,t)
On the right side of the equation(12) containing the unknown Hessian matrix H[J2 (γˆt−1, βˆt−1, cˆt−1)] , extended generalized information matrix Ξ (t) and innovation vector E( p, t) , and φ(t) contains the unpredictable noise v(t − i), i = 1, 2,, nd .In order to solve these difficulties, according to the principle of recursive identification, replacing H[J2 (γˆt−1, βˆt−1, cˆt−1)] , Ξ (t) and E( p, t) in the above Eq.(12) with H[J2 (γˆt−1, βˆt−1, cˆt−1)] , Ξˆ (t) and Eˆ ( p, t) , let vˆ(t − i) denote the estimate of v(t − i) to define the estimate of φ(t) as follows: φˆ (t) := [− y(t −1), − y(t − 2),, − y(t − na ),
vˆ(t −1), vˆ(t − 2),, vˆ(t − nd )]T From Eq.(6), v(t) can be written as: (10) (11) (12) (14) Πˆ (t)= 12{Hˆ [J2(γˆt−1, βˆt−1,cˆt−1)]} Φˆ0T( p,t)Φˆ0( p,t) Φˆ0T( p,t)Φ(ct−1,t) Φˆ0T( p,t)Ψ(βt−1,t) =ΦT(ct−1,t)Φˆ0( p,t) ΦT(ct−1,t)Φ(ct−1,t) Πˆ 23(t)
ΨT(βt−1,t)Φˆ0( p,t) Πˆ 23T(t) ΨT(βt−1,t)Ψ(βt−1,t) Πˆ 23(t)={G(t − i) −y(t − i) + φˆT (t − i)γˆt−1 + βt−1TG(t − i)cˆt−1} p−1 i=0 + ΦT (cˆt−1,t)Ψ(βt−1,t) Ξˆ (t)= Φˆ 0T ( p, t) ΦT (cˆt−1, t) Ψ T ( βˆt−1, t)
T Eˆ ( p, t) = Y ( p, t) − Φˆ 0 ( p, t)γˆt−1 − Φ(cˆt−1, t) βˆt−1 Y ( p, t)= [ y(t)，y(t −1)，，y(t − p + 1)] Φˆ 0 ( p, t)= [φˆ (t)，φˆ (t-1)，，φˆ (t − p + 1)]
T T T Ψ(βˆt−1,t)= GT (t)βˆt−1,GT (t −1)βˆt−1,GT(t − p +1)βˆt−1 Φ(cˆt−1,t)=[G(t)cˆt−1,G(t −1)cˆt−1,,G(t − p +1)cˆt−1]
T φˆ (t) = [− y(t −1), − y(t − 2),, − y(t − na ),
vˆ(t −1), vˆ(t − 2),, vˆ(t − nd )]T cˆt = sgn[θˆna +nb +nd +1 (t)] t
[θˆ(t)](na + nb + nd + 1 : n) [θˆ(t)](na + nb + nd + 1 : n) [θˆ(t)](na + nb + nd + 1: n) = cˆtt (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25)
The inverse matrix Ωˆ −1 (t) in Eq. (13), for all t ,the stacked data length p in the nonsingular matrix Ωˆ (t) should be large enough to invert the nonsingular matrix. The process of computing θˆ(t) by the H-ENR algorithm is summarised as follows:
(1) Choose the stacked data length p and initialize: let t = 1, θ（ˆ0）be an arbitrary real vector with cˆ0 = 1 .
(2) Collect the measured data u(t) and y(t) ,form stacked vector Y ( p, t) by Eq.(19) , G(t) by Eq.(5) the information vector φˆ(t) by Eq.(23) and Φˆ 0 ( p, t) by Eq.(20).
(3) Compute and form Ψ (bˆt−1, t) by Eq. (21) and Φ(cˆt−1, t) by Eq. (22). (4) Form information matrix Ξˆ (t) by Eq. (17) and compute innovation vector Eˆ ( p, t) by Eq. (18). (5) Compute Πˆ 23 (t) by Eq. (16) and Πˆ (t) by Eq. (15). (6) Update the parameter estimation θ（ˆt）by Eq. (13). (7) Normalize cˆt by Eq. (24) and Eq. (25) with the first positive element.
(8) Increase t by 1 and go to Step 2. 3.2
Consider the following Hammerstein nonlinear system: A( z) y(t) = B( z)μ (t) + D( z)v(t) A( z) = 1 +α 1z −1 +α 2 z −2 = 1 −1.07 z −1 + 0.675z −2 ,
B( z) = β1z −1 + β 2 z −2 = 1.55z −1 + 1.20z −2 ,
In simulation, the input {u(t)} is taken as a persistent excitation signal sequence,the noise {v(t)} is taken as a white noise sequence with zero mean and and variance σ 2 = 0.302 , and the data length is taken as p=100 and p=160 .Adopting the the Newton extended recursive algorithm(H-ENR) to estimate the paramelers of this Hammerstein-CARMA process,the corresponding noise-to-signal ratio is δns = 8.76%，where the noise-to-signal ratio δns is as follows( h(t) and x(t) in Figure 1): δ ns = h(t)=
D(z) A(z) var [h(t)] var [ x(t)] v(t) , x(t)=
μ (t) . ×100% ,
B(z) A(z)
The parameter estimates and their errors are shown in Table 1 , the estimation errors versus t are shown in Figure 2.
From Table 1 and Figure 2, we can get the following conclusions: (1) For the confirmed Hammerstein-CARMA model( v(t) =0 orσ 2 =0 ),the Newton extended recursive algorithm can converge to the true value faster than the extended projection algorithm. For the stochastic Hammerstein-CARMA model ( σ 2 ≠ 0,σ 2 =0.302 ),the parameter estimation of Newton extended recursive algorithm fluctuates greatly, especially for the small stacked data length p , and its estimation errors cannot converge to zero even if the data length t tends to infinity. The reason is that the increment of the extended recursive algorithm does not approach zero. However, when the length of stacked data length p increases, the parameter estimation will become getting more stationary, as shown Table 1 and Figure 2 with p=100 and p=160 . (2) The parameter estimation errors rapidly converges to a small constant as the data length t increases. As the data length t goes to infinity, this constant is going to get very small and close to zero. This shows that the extended Newton recursive algorithm is effective. 4
This paper studies the parameter estimation methods for the nonlinear Hammerstein-CARMA model. A extended Newton recursive(H-ENR) algorithms are derived based on the Newton method. Aiming at the difficulty that the information vector of Hammerstein-CARMA model contains unmeasured noiseterms, the principle of recursion identification is applied. The unknown noise terms contained in the information vector are replaced by its estimated value, and the estimated value is calculated by the parameterestimated value of the previous time or the previous time. Compared with the extended stochastic gradient algorithms the H-ENR algorithm has improved parameter estimation accuracy.The numerical example shows that the parameter estimates for the proposed H-ENR algorithm converge to their true values.At present, there is a lot of work tobe done in the study of nonlinear systems. ln this paper, only single-input single-output nonlinearsystems are studied. How to extend it to multi-inputmulti-output nonlinearsystems and apply it in the field is the next problem to be considered. p
This work is supported by the Fundamental Research Funds for the Central Universities under Grant 31920210075. [1] L. Ljung, System Identification: Theory for the User, 2nd ed., Prentice-Hall, Englewood Cliffs, NJ, 1999. [2] F. Ding, System Identification - New Theory and Methods, Science Press, Beijing, 2013. [3] Y.Cao, P. Li，Y. Zhang, Parallel processing algorithm for railway signal fault diagnosis data