# Instrumental Variable-Based OMP Identification Algorithm for Hammerstein Systems

Authors: Shuo Zhang; Dongqing Wang; Yaru Yan
Journal: Complexity
 (2018)

---

## Abstract

Hammerstein systems are formed by a static nonlinear block followed by a dynamic linear block. To solve the parameterizing difficulty caused by parameter coupling between the nonlinear part and the linear part in a Hammerstein system, an instrumental variable method is studied to parameterize the Hammerstein system. To achieve in simultaneously identifying parameters and orders of the Hammerstein system and to promote the computational efficiency of the identification algorithm, a sparsity-seeking orthogonal matching pursuit (OMP) optimization method of compressive sensing is extended to identify parameters and orders of the Hammerstein system. The idea is, by the filtering technique and the instrumental variable method, to transform the Hammerstein system into a simple form with a separated nonlinear expression and to parameterize the system into an autoregressive model, then to perform an instrumental variable-based orthogonal matching pursuit (IV-OMP) identification method for the Hammerstein system. Simulation results illustrate that the investigated method is effective and has advantages of simplicity and efficiency.


---

## Body

## 1. Introduction

Nonlinear system modeling and identification are very important in theory and application [1–6]; block-oriented nonlinear systems, which combine nonlinear and linear blocks in various styles, are the typical representation of nonlinear systems. Hammerstein systems, a static nonlinear block plus a dynamic linear block, are a type of commonly used block-oriented nonlinear systems, and a lot of work has contributed in identification methods for Hammerstein systems. The traditional identification methods for Hammerstein systems mainly include the overparameterization model-based methods [7–10] and the recursive/iterative identification methods [11–14]. The highly efficient identification methods include the key-term separation principle-based identification methods [15–17], the hierarchical identification methods [18, 19], the filtering technique based-identification methods [20–22], the maximum likelihood estimation methods [23–25], and the evolution optimization methods [26–28]. Recently, Chen et al. investigated a particle swarm optimization method [29], Krishnanathan et al. discussed a continuous-time nonlinear systems using approximate Bayesian computation [30], and Wang et al. studied a model recovery for Hammerstein systems using the auxiliary model-based orthogonal matching pursuit method [31] and so on.It is difficult to parameterize Hammerstein systems due to existing parameter coupling between the nonlinear part and the linear part of Hammerstein systems. In order to solve this problem, an instrumental variable-based method is studied to parameterize the Hammerstein systems. Further, to achieve in simultaneously identifying parameters and orders and to promote the computational efficiency of the estimated method, a sparsity-seeking orthogonal matching pursuit optimization method of compressive sensing is extended to identify parameters and orders of the Hammerstein systems. The idea is, by the filtering technique and the instrumental variable method, to transform the Hammerstein system into a simple form with a separated nonlinear expression, so as to easily parameterize the system into autoregressive form, and to perform the instrumental variable-based orthogonal matching pursuit (IV-OMP) identification method for the Hammerstein systems.In the compressive sensing theory [32–34], there mainly exist two sparse signal recovery methods: the orthogonal matching pursuit algorithms and the basis pursuit algorithms. The orthogonal matching pursuit (OMP) algorithm is a kind of greedy parameter recovering method [35–37], it selects the best fitting column of the measurement matrix and the corresponding sparse signal in each selected step. Due to the selection being orthogonal, the OMP algorithm has a lower computational complexity compared with the basis pursuit algorithms [38–40].Recently, Mao et al. investigated parameter estimation algorithms for Hammerstein time-delay systems based on the OMP scheme to estimate the system parameters and the time delay [41] by means of the compressed sensing recovery theory and the auxiliary model identification idea. In this paper, an instrumental variable-based orthogonal matching pursuit (IV-OMP) algorithm is investigated to simultaneously estimate the orders and parameters of a Hammerstein system. The contributions lie in four aspects:
(i)
To solve the parameterizing difficulty caused by the coupling between the nonlinear part and the linear part in Hammerstein systems, a filtering technique-based instrumental variable method is studied to parameterize the Hammerstein systems.
(ii)
To achieve in simultaneously identifying parameters and orders and to promote the computational efficiency of the system, an instrumental variable-based orthogonal matching pursuit (IV-OMP) optimization method of compressive sensing is extended to identify parameters and orders of the systems.
(iii)
The advantage of the proposed IV-OMP method over the traditional methods is that it is not necessary to collect a lot of data and invest a lot of power into the parameter identification based on the sparse principle.
(iv)
Simulation results illustrate that the investigated method is effective and has advantages of simplicity and efficiency.The rest of the paper is organized as follows. Section2 demonstrates the identification problem of a Hammerstein system. Section 3 studies the parameterizing method of the Hammerstein system. Section 4 presents the IV-OMP identification algorithm by using the instrumental variables. Section 5 provides a numerical example for the proposed method. Finally, the concluding remarks are involved in Section 6.


## 2. The Problem Formulation

For the narrative convenience, we define some notation. “M≕N” represents “M is defined as N”; z represents a unit forward shift operator: zyt=yt+1 and z−1yt=yt−1;x^t stand for the estimate of x at time t. The input nonlinear and output linear functions of a Hammerstein CARMA system in Figure 1 are expressed as
(1)xt=fut=∑k=1ncckfkut,(2)Azyt=Bzxt+Dzvt,where ut and yt are the system input and output, xt is an internal variable, and vt is stochastic white noise with zero mean; the input nonlinearity f is modeled as a linear combination of basis functions fk; nc is the number of the basis functions; the linear block is a CARMA model; Az, Bz, and Dz are polynomials in the unit backward shift operator z−1z−1yt=yt−1 and defined by
(3)Az≔1+a1z−1+a2z−2+⋯+anaz−na,Bz≔1+b1z−1+b2z−2+⋯+bnbz−nb,Dz≔1+d1z−1+d2z−2+⋯+dndz−nd.Figure 1
A system described by the Hammerstein CARMA model.If the orderna is known and the orders nb, nc, and nd are unknown, then we face two parameterizing problems for the Hammerstein system:
(i)
The coupling between the nonlinear part and the linear part in a Hammerstein system makes the parametrization of the Hammerstein system more difficult.
(ii)
Simultaneous and efficient identification of system parameters and orders is difficult to be achieved, due to existing unknown orders.


## 3. The Parametrization of The Hammerstein CARMA System

Divide both sides of (2) by the filtering function Bz, we get
(4)AzBzyt=xt+DzBzvt.Define the instrumental variableszt and wt as
(5)zt≔AzBzyt,(6)wt≔DzBzvt.Then, (4) can be written as
(7)zt=xt+wt.Replacing the expressions ofAz and Bz and Bz and Dz into (5) and (6), respectively, we can get
(8)zt=−∑i=1nbbizt−i+yt+∑i=1naaiyt−i,(9)wt=−∑i=1nbbiwt−i+vt+∑i=1nddivt−i.Substitutingzt and wt in (8) and (9) and xt in (1) into (7) gives
(10)yt=−∑i=1naaiyt−i+∑i=1nccifiut+∑i=1nbbizt−i−wt−i+∑i=1nddivt−i+vt.Since the ordersnb, nc, and nd are unknown, we assume a sufficient length of order L for nb, nc, and nd (L>nb,L>nc,L>nd). Then, (10) can be rewritten as
(11)yt=−∑i=1naaiyt−i+∑i=1Lcifiut+∑i=1Lbizt−i−wt−i+∑i=1Ldivt−i+vt.Define the information vectorsφt, φut, and φvt as
(12)φt≔φutφvt∈ℝN,N≔3L+na,φut≔−yt−1,−yt−2,…,−yt−na,f1ut,f2ut,…,fLutT∈ℝL+na,φvt≔zt−1−wt−1,zt−2−wt−2,…,zt−L−wt−L,vt−1,vt−2,…,vt−LT∈ℝ2Land the parameter vectors Θ, θu and θv, as
(13)Θ≔θuθv∈ℝN,θu≔a1,a2,…,ana,c1,c2,…,cnc,0,0,…,0⏟L−ncT∈ℝL+na,θv≔b1,b2,…,bnb,0,0,…,0⏟L−nb,d1,d2,…,dnd,0,0,…,0⏟L−ndT∈ℝ2L.Then, (11) can be written as
(14)yt=φuTtθu+φvTtθv+vt=φTtΘ+vt.Samplingm sets of data t=1,2,…,m and substituting them into (14) get
(15)y1=φT1Θ+v1,y2=φT2Θ+v2,…,ym=φTmΘ+vm.Define the accumulated information vectors and the matrix,(16)Y≔y1y2⋮ym∈ℝm,V≔v1v2⋮vm∈ℝm,Φ≔φT1φT2⋮φTm∈ℝm×N,then (15) can be described as
(17)Y=ΦΘ+V.If there are enough measurements, that is,m reaches several thousands, we can get the least squares estimate of Θ, as follows:
(18)Θ^LS=ΦTΦ−1ΦTY.But with an assumed big order lengthL (≥5) and N, the above least squares algorithm leads a big computational burden and is not suitable for solving system parameters and orders on line. To avoid this problem and get an efficient identification algorithm, we choose a more advantageous method over the traditional least squares method; the aim is to identify the parameter vector Θ with less observations (K<m<n) by using the OMP theory based on the mentioned instrumental variable method.According to the CS theory, we can regard the parameter vectorΘ as a sparse signal. Let Θ0=K be the number of the nonzero entries in Θ, that is, the sparsity level of Θ, then the identification problem can be described as an orthogonal matching pursuit (OMP) method:
(19)Θ^=argminΘ0,s.t.y−φTtΘ2<ε,where Θ^ is the estimate of Θ and εε>0 is the error tolerance.


## 4. The IV-OMP Identification Algorithm

For the sparse parameter vectorΘ, if K=Θ0K<N, there will be K nonzero scalar components in Θ. What we should do in this identification method is to pick out and recover valid data in Θ from Φ by using the OMP theory.Defineϕi as the ith column vector of Φ (a ϕi is also called an atom), and θi as the ith element of Θi=1,2,…,N, that is,
(20)Φ=ϕ1,ϕ2,…,ϕN,Θ=θ1,θ2,…,θNT.Then, the vectorY in (17) can be written as
(21)Y=∑i=1Nϕiθi+V.Obviously, the output vectorY contains a linear combination of all atoms ϕii=1,2,…,N. The main idea is to find the K nonzero items and corresponding valid atoms at the right-hand side of (21).To describe the IV-OMP recovery algorithm, we need to give some notations. Letk=1,2,… be the iterative number and λk be the index of the solution support ϕλk in Φ at the kth iteration; Λk is a set composed of λi, i=1,2,…,k; ΦΛk is the submeasurement matrix composed by the k columns of Φ indexed by Λk; rk denotes the residual at the kth iteration; and Θ^k is the parameter estimation at the kth iteration. The algorithm is initialized as r0=Y,Λ0=∅ and Θ^Λ0=0.Define a cost function at thekth iteration
(22)Jθi=rk−1−ϕiθi2,i=1,2,…,N.Remark 1.
MinimizingJθi with respect to θi means that the derivation of Jθi with respect to θi is equal to 0, that is,
(23)∂Jθi∂θi=−2ϕiTrk−1−ϕiθi=0,and we get
(24)θi=ϕiTrk−1ϕi2.
Substitute the aboveθi into (22) to get a minimized Jθi as
(25)Jθimin=ϕiϕiTrk−1ϕi2−rk−12=ϕiϕiTrk−1ϕi2−rk−1TϕiϕiTrk−1ϕi2−rk−1=rk−1TϕiϕiTϕi2−rk−1TϕiϕiTrk−1ϕi2−rk−1=rk−1TϕiϕiTϕiϕiTrk−1ϕi2ϕi2−rk−1TϕiϕiTrk−1ϕi2−rk−1TϕiϕiTrk−1ϕi2+rk−1Trk−1=rk−12−ϕiTrk−12ϕi2.Remark 2.
The resultJθimin=rk−12−ϕiTrk−12/ϕi2says that minimizing Jθi is equivalent to the quest for the largest inner product between the residual rk−1 and the normalized column vector ϕi of Φ.The next step is to find the indexi and the column ϕi corresponding to the largest inner product between the residual rk−1 and the normalized column vectors ϕi and assign them into λk and ϕλk at the kth step:
(26)λk=argmaxi=1,2,…,Nrk−1,ϕiϕi.Remark 3.
Minimizing the least squares cost functionJΘΛk=Y−ΦΛkΘΛk2 is equal to taking the derivation of JΘΛk with respect to ΘΛk,
(27)∂JΘΛk∂ΘΛk=ΦΛkTY−ΦΛkΘΛk=ΦΛkTrk=0.
It means that the residualrk is orthogonal to the subinformation matrix ΦΛk. Therefore, the index corresponding to the largest inner product between the residual rk−1 and the normalized column vector ϕi is computed by
(28)λk=argmaxi∈S\Λk−1rk−1,ϕiϕi,S=1,2,…,N.Update the support setΛk and the subinformation matrix ΦΛk by
(29)Λk=Λk−1,λk,ΦΛk=ΦΛk−1,ϕλk.With the obtained index setΛk and the corresponding subinformation matrix ΦΛk, the least squares estimation for the nonzero parameters at the kth iteration can be achieved. Define a cost function:
(30)JΘΛk=Y−ΦΛkΘΛk2.MinimizingJΘΛk should get the least squares estimates of the nonzero parameters. But there exist the unknown instrumental variables zt−i and wt−i and the unmeasurable noise term vt−i in the information vector φt of Φ; the parameter vector ΘΛk cannot be estimated by using the standard OMP methods in the CS theory. The method is to use their estimates z^t−i, w^t−i, and v^t−i to replace them; the computation of the estimates z^t−i, w^t−i, and v^t−i is as follows:
(i)
Expand the ordersnb, nc, and nd in (1), (2), (8), and (9).
(ii)
z^t−i is computed by replacing ai and bi with their estimates a^i,kandb^i,k in (8).
(iii)
v^t−i is computed by replacing ai, bi, ci, and di with their estimates a^i,k, b^i,k, c^i,k, and d^i,k in (2) with (1) being inserted.
(iv)
w^t−i is computed by replacing bi and di with their estimates b^i,k and d^i,k in (9). The results are
(31)z^kt=−∑i=1Lb^i,kz^kt−i+∑i=1naa^i,kyt−1+yt,v^kt=∑i=onaa^i,kyt−i−∑i=0Lb^i,k∑i=1Lc^i,kfkut−i−∑i=1Ld^i,kv^t−i,b^i,0=1w^kt=−∑i=1Lb^i,kw^kt−i+∑i=1Ld^i,kv^t−i+vt,and the estimates of φvt, φt, and Φ at iteration k can be written as
(32)φ^v,kt≔z^t−1−w^t−1,z^t−2−w^t−2,…,z^t−L−w^t−L,v^t−1,v^t−2,…,v^t−LT∈ℝ2L,φ^kt≔φutφ^v,kt≕−yt−1,−yt−2,…,−yt−na,f1ut,f2ut,…,fLut,z^kt−1−w^kt−1,z^kt−2−w^kt−2,…,z^kt−L−w^kt−L,vt−1,vt−2,…,vt−LT,Φ^k≔φ^kT1φ^kT2⋮φ^kTm≕ϕ^1,k,ϕ^2,k,…,ϕ^N,k.With replacing the above unknown variables with their estimates, the least squares estimateΘ^Λk for the parameter vector ΘΛk at the k step is given:
(33)Θ^Λk=Φ^ΛkTΦ^Λk−1Φ^ΛkTY.Remark 4.
At the beginning, the IV-OMP algorithm with the estimatesz^t−1, w^t−i, and v^t−i in the subinformation matrix Φ^Λk causes an inaccurate support atom selecting. With the iteration k increasing, these estimated variables become more accurate, the misselected support atoms certainly unmeet the threshold requirement, and the corresponding element in the estimated parameter support set will be a small nonzero value. Thus, we set an appropriate small threshold ε to filter the parameter estimate Θ^Λk. If θ^h,k<ε (where θ^h,k is the hth element of Θ^ΛK), eliminate θ^h,k from Θ^ΛK and the corresponding ϕ^h,k from Φ^Λk,
(34)Θ^Λkε=Θ^Λkεθ^h,k∉Θ^Λkε,Φ^Λkε=Φ^Λkϕ^h,k∉Φ^Λk.
Becauseθ^h,k is the hth element in Θ^Λk and ϕ^h,k is the hth column in Φ^Λk, the expressions of eliminating θ^h,k from Θ^Λk and ϕ^h,k from Φ^Λk in MATLAB are
(35)Θ^Λkε:Θ^Λkh=,Φ^Λkε:Φ^Λkh=.
Then we have(36)rk=Y−Φ^ΛkεΘ^Λkε.Briefly, the proposed IV-OMP algorithm can be implemented as in Algorithm1.Algorithm 1: (IV-OMP algorithm).
IV-OMP for CS Recovery
Measurement matrix:Φ=ϕ1,ϕ2,…,ϕN,.
output vector:Y, sparsity level: K,
parameter vectorΘ=θ1,θ2,…,θNT.
Initialization:k=0,r0=Y,Λ0=∅,Θ^Λ0=0,Θ^k=0. repeat
k=k+1;
Formφ^kt by Equation (32);
Φ^k=φ^kT1φ^kT2⋮φ^kTm=ϕ^1,k,ϕ^2,k,…,ϕ^N,k;
λk = index of the highest amplitude component of ϕ^i,kTrk−1;
Λk=Λk−1,λk;
Θ^Λk=Φ^ΛkTΦ^Λk−1Φ^ΛkTY;
Ifθ^h,k<ε, let Θ^Λkh= to get Θ^Λkε,.
and letΦ^Λkh= to get Φ^Λkε;
rk=Y−Φ^ΛkεΘ^Λkε;
letΘ^kΛkε=Θ^Λkε to get Θ^k;
a^:,k=Θ^k1:na,c^:,k=Θ^kna+1:na+L,
b^:,k=Θ^kna+L+1:na+2L,
d^:,k=Θ^kna+2L+1:na+3L;
z^kt=−∑i=1Lb^i,kz^kt−i+∑i=1naa^i,kyt−i+yt,
v^kt=∑i=0naa^i,kyt−i−∑i=0Lb^i,k∑i=1Lc^i,kfkut−i
−∑i=1Ld^i,kv^t−i,
w^kt=−∑i=1Lb^i,kw^kt−i+∑i=1Ld^i,kv^t−i+vt;
until k = a fixed positive integer k>K, s,.t. rk<ε;.
Estimated parameter vector:Θ^Λkε and Θ^k with Θ^Λ¯kε=0.


## 5. Examples

Consider the following Hammerstein CARMA system(37)Azyt=Bzxt+Dzvt,xt=−1.70ut+0.80u2t,Az=1+1.20z−1+0.60z−2,Bz=1+0.12z−1−0.06z−2,Dz=1+0.11z−1−0.10z−2,Θ=1.20,0.60⏟na=2,−1.70,0.80,0,0,0⏟L=5,0.12,−0.06,0,0,0⏟L,0.11,−0.10,0,0,0⏟LT.In the simulation, the inputut is taken as an uncorrelated persistently excited signal vector sequence with zero mean and unit variance, and vt is taken as a white noise sequence with zero mean and variance σ2=0.102. Applying the proposed IV-OMP algorithm to estimate the parameters and the orders of this system, the parameter estimates and their errors are shown in Table 1 and Figures 2 and 3. The relative error of the parameters is
(38)δ≔Θ^k−ΘΘ×100%.Table 1
the recursive parameter estimates and errors by using different algorithms.
Algorithm
k
a^1,k
a^2,k
c^1,k
c^2,k
c^3,k
c^4,k
c^5,k
b^1,k
b^2,k
b^3,k
b^4,k
b^5,k
d^1,k
d^2,k
d^3,k
d^4,k
d^5,k
δ (%)
IV-OMP
1
0.7717
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
87.52
2
0.7100
0
0
0
0
0
−0.2365
0
0
0
0
0
0
0
0
0
0
88.71
3
0.824
0
0
0
0
0
−0.2875
0
0
0
0
0
0
0
0.5002
0
0
90.55
4
1.304
0.7615
0
0
0
0
−0.2439
0
0
0
0
0
0
0
0.0488
0
0
82.65
5
1.331
0.8617
0
0
0
0.2315
−0.2468
0
0
0
0
0
0
0
−0.0846
0
0
83.87
6
1.399
0.8219
0
0
0
0.1992
−0.2422
0
0.3025
0
0
0
0
0
−0.2355
0
0
85.66
7
1.086
0.5736
−1.857
0
0
0.3747
0.0527
0
−0.0982
0
0
0
0
0
0.0080
0
0
40.00
8
1.126
0.5786
−1.668
0.8567
0
−0.0220
0
0
−0.0314
0
0
0
0
0
0.0172
0
0
9.46
9
1.124
0.5833
−1.693
0.8531
0
0
0
0
−0.0362
0
0
0
0
0
0.0126
0
0
9.27
10
1.124
0.5855
−1.686
0.8181
0
0
0
0
−0.0370
0
0
0
0
0
0.0104
0
0
9.02
11
1.125
0.5732
−1.697
0.8275
0
0
0
0
−0.0531
0
0
0
0
−0.0238
0.0126
0
0
8.58
12
1.122
0.5754
−1.690
0.8295
0
0
0
0.0303
−0.0550
0
0
0
0
−0.0477
0.0038
0
0
7.54
13
1.163
0.5926
−1.688
0.8159
0
0
0
0.0864
−0.0610
0
0
0
0.0593
−0.0686
−0.0090
0
0
3.48
14
1.187
0.5913
−1.682
0.8096
0
0
0
0.0976
−0.0680
0
0
0
0.0858
−0.0756
−0.0016
0
0
2.18
15
1.196
0.5973
−1.659
0.8035
0
0
0
0.1014
−0.0626
0
0
0
0.0929
−0.0726
0
0
0
2.46
True values
1.20
0.60
−1.70
0.80
0
0
0
0.12
−0.06
0
0
0
0.11
−0.10
0
0
0Figure 2
The parameter estimation errors versusδ.Figure 3
The parameter estimates.From Table1 and Figures 2–3, we can get the following:
(1)
The parameter estimation errors become (generally) smaller and smaller with the iterationk increasing.
(2)
There exist 3 wrongly selected atoms corresponding to parametersc^4,k, c^5,k, and d^3,k (black-colored lines in Figure 3). With the iteration k increasing, the wrongly chosen atoms are deleted, and the parameter estimation errors become smaller.


## 6. Conclusions

This paper is aimed at exploring the identification of both the parameters and orders of the Hammerstein CARMA system with limited sampling data. To solve the parameterizing difficulty caused by parameter coupling between the nonlinear part and the linear part in a Hammerstein system, firstly, by filtering the equation of the linear block with the coefficient function of the controlled term, we separate the parameter coupling between the linear block and the nonlinear block. Moreover, by using two instrumental variables, the Hammerstein system is parameterized into an autoregressive form. To achieve in simultaneously identifying parameters and orders and to promote the computational efficiency of the identification algorithm, an instrumental variable-based orthogonal matching pursuit (IV-OMP) optimization method of compressive sensing is extended to identify parameters and orders of the Hammerstein system. Simulation results illustrate that the investigated method is effective and has advantages of simplicity and efficiency. The proposed IV-OMP optimization method can be extended to the colored noise systems, the networked dynamic systems [42–46], and so on.



---

*Source: 8420426-2018-07-22.xml*