# Subspace Identification of Hammerstein Model with Unified Discontinuous Nonlinearity

Authors: Borhen Aissaoui; Moêz Soltani; Abdelkader Chaari
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The main aim of this study is to handle the case where the structures of nonlinear systems are unknown. In the many works, the parametric identification of nonlinear systems represented by Hammerstein model, with discontinuous and asymmetric nonlinearity, considers the structures of the nonlinear and linear blocks are known, especially the nonlinear bloc. To solve this problem, a unified form of nonlinearity representing eight cases of nonlinearities can be used. The parameters of both blocks, linear and nonlinear, are estimated using an iterative subspace approach. More importantly, in an attempt to show the extent to which this method is efficient, we apply it to experimental data obtained from the electropneumatic system. As a result, the numerical and experimental examples confirm a good conditioning and computational efficiency.


---

## Body

## 1. Introduction

Many physical and biological systems have been modelled as a Hammerstein structure, such as the Stretch Reflex at the human ankle [1], the mechanical behaviour of lung tissue [2] and the electropneumatic system. In the Hammerstein model, there are two cascaded linear and nonlinear blocks which are an interconnection between the static structure nonlinearity and the Linear Time-Invariant (LTI) dynamic subsystem. Indeed, a variety of methods have been developed to identify Hammerstein systems. Narendra and Gallman [3] developed an iterative identification algorithm. They separate the parameters into two sets: one corresponding to the nonlinear component and the second from the linear element. In addition, the cross correlation-based methods are applied in order to estimate the linear dynamic subsystem [4] where the algorithm updates the nonlinearity’s output estimation at each step using the inverse dynamic of the linear component. Finally, they estimate and predict the inverse dynamics and identify the output nonlinearity. In other studies, the structures of the nonlinear and linear blocks are known in the case of continuous nonlinearity [5], discontinuous nonlinearities (DN) such as saturation with hysteresis (SH) [6], hysteresis-backlash [7], and hysteresis-relay nonlinearities [8].In the real world, the identification of static nonlinearity subsystem remains the major problem in the Hammerstein model. Indeed, the structures of the nonlinear block are unknown and their estimate is important. Rejeb et al. [9] used the least squares method to estimate the parameters of unknown nonlinear block. Later, the same authors [10] proposed a general model combining discontinuous and asymmetric nonlinearities which generates all the possible combinations of the nonlinearities discontinuous elements. In some context, Vörös [11] developed a system with the general backlash which is based on appropriate switching and internal functions.The identification using parametric approaches requires a minimal information priority about the system structure, that is, system order, noise model, and so forth [12]. Then, the LTI state model cannot be estimated through input-output measurements. Hence, the subspace method is introduced in order to overcome this problem. The subspace methods estimate the state space models for linear system with no a priori knowledge about the system [13]. In other words, those approaches are characterized by the simplicity of implementation including its effective results. In many cases, these methods provide a good alternative to the classical nonlinear optimization-based prediction-error methods [14].Subspace identification is quite a well-accepted method for the identification of multivariable LTI systems. This method estimates an LTI state space model directly from input-output measurements. Moreover, Verhaegen and Varga [15] as well as Verhaegen and Westwick [16] proposed the Multivariable Output Error State Space identification model (MOESP) algorithm to identify Hammerstein models. In [1, 17, 18], the authors have treated the case of Single Input Single Output (SISO) Hammerstein model. Other authors proposed some subspace methods adapted to the Hammerstein-Wienner [19, 20]. Moreover, many researchers have developed a new identification subspace-based algorithm for the linear model [21–23].Nevertheless, the structure of nonlinear block is used as known in those methods. Indeed, we developed an identification algorithm for a SISO Hammerstein model with discontinuous nonlinearities such as hysteresis, saturation, preload, and dead-zone. The proposed method uses an iterative subspace algorithm. The determination of the structure of nonlinearity is obtained based on Unified Discontinuous Nonlinearity (UDN) method. The parameters of both blocks are estimated simultaneous using the proposed iterative subspace approach.This paper is organized as follows. Section2 presents the problem formulation. Section 3 describes the proposed method. A numerical simulation validates the proposed algorithm and then applying to estimate the linear and nonlinear blocks model of electropneumatic system are shown in Section 4. Finally, this paper will be concluded in Section 5.


## 2. Problem Formulation

Figure1 presents a discrete-time SISO Hammerstein system, which is described by the following forms: (1)S:xk+1=Axk+Bvk+ωk,yk=Cxk+Dvk+εk,where y(k)∈R, x(k)∈Rn, v(k)∈R, ω(k)∈Rn, and ε(k)∈R are the LTI system output, state, output of the nonlinearity, process noise, and output measurement noise vectors at discrete time k. A∈Rn×n, B∈Rn×1, C∈R1×n, and D∈R are the unknown system matrices of the state space.Figure 1
A discrete-time Hammerstein model.The discontinuous and asymmetric nonlinearityN[·] can be characterized by a linear function as (2)vk=Nuk=gTφuk,with(3)g=g1,…,grT,φuk=φ1uk,…,φrukT,where g∈Rr×1 is the unknown nonlinearity parameter vector and φu(k)∈Rr×1 is the observation vector.In our work, the process noiseω(k) is considered zero. Then, replacing (2) in (1), we obtain(4)S:xk+1=Axk+BgTφuk,yk=Cxk+DgTφuk+εk.The identification problem is how to estimate the nonlinear parameter vector g, the dynamic matrices A,B,C,D, and the model order n from input-output data set.


## 3. Methods

In this section, we present a subspace identification algorithm for a SISO Hammerstein model with Unified Discontinuous Nonlinearities.


### 3.1. Unified Discontinuous Nonlinearity

Figure2 illustrates the v(k) evolution with respect to the input u(k) of the nonlinear block using UDN [9]. This latter contains an asymmetric piecewise function with a hysteresis, a saturation, a preload, a dead-zone, and so forth.Figure 2
Unified Discontinuous Nonlinearity.The UDN is characterized by the unknown parametersp, a, b, Z1, Z2, e1, and e2. Z1 and Z2 are scalar positive parameters that characterize the saturation, a≥0 and b≥0 define the width and the centre of the hysteresis as well as the dead-zone width, p is the hysteresis slope. e1 and e2 define the preload thresholds. The different choices of these parameters lead to the definition of eight nonlinearities which are listed in Table 1 and shown in Figure 3.Table 1
Different generated nonlinearities.
DN Choices of parameters UDN p ≠ 0, a≠0, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN1 p ≠ 0, a=b, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN2 p ≠ 0, a=b, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN3 p ≠ 0, a≠0, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN4 p ≠ 0, a≠0, b=0, Z1≠0, Z2≠0, e1=0, e2=0 DN5 p ≠ 0, a=0, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN6 p ≠ 0, a=0, b=0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN7 p ≠ 0, a=0, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN8 p ≠ 0, a=0, b=0, Z1≠0, Z2≠0, e1=0, e2=0Figure 3
Different cases of nonlinearities.Using Table1, the general form of v(k) with respect to u(k) and Δu(k) for the different cases of nonlinearities can be defined as(5)vk=pZ1+e1ifuk>Z1+b+a,Δuk<0orZ1+b-a<uk<Z1+b+a,Δuk<0puk-a-b+e1ifb+a≤uk≤Z1+b+a,Δuk>0puk-a+b+e2ifZ2-b+a≤uk≤-b+a,Δuk>0puk+a-b+e1ifb-a≤uk≤Z1+b-a,Δuk<0puk+a+b+e2ifZ2-b-a≤uk≤-b-a,Δuk<0pZ2+e2ifuk<Z2-b-a,Δuk<0orZ2-b<uk<Z2-b+a,Δuk>0e1ifb-a<uk<b+a,Δuk>0e1if-b-a<uk<-b+a,Δuk<00if-b+a≤uk≤b-a∀Δuk,with α=h[Δu(k)] and Δu(k)=u(k)-u(k-1).Using (2), (5) can be rewritten as(6)vk=pukf1k-paf2k+pbf3k+pZ1h1k+pZ2h2k+e1h1prk+e2h2prk=gUDNTφUDNk,with(7)f1k=α1-huk-Z1+b+a+hb+a-uk+1-hZ2-b+a-uk+huk--b+a+1-α1-huk-Z1+b-a+hb-a-uk+1-hZ2-b-a-uk+huk--b-a,f2k=α1-huk-Z1+b+a+hb+a-uk+1-hZ2-b+a-uk+huk--b+a-1-α1-huk-Z1+b-a+hb-a-uk+1-hZ2-b-a-uk+huk--b-a,f3k=α1-huk-Z1+b+a+hb+a-uk-1-hZ2-b+a-uk+huk--b+a+1-α1-huk-Z1+b-a+hb-a-uk-1-hZ2-b-a-uk+huk--b-a,h1k=αhuk-Z1+b+a+1-αhuk-Z1+b-a,h2k=αhZ2-b+a-uk+1-αhZ2-b-a-uk,h1prk=huk-b-a,h2prk=h-b+a-uk,hβ=1ifβ>00ifβ≤0,where hβ is a switching sequence.Using (6), (3) are rewritten as(8)gUDNT=p,pa,pb,pZ1,pZ2,e1,e2,φUDNk=ukf1k,-f2k,f3k,h1k,h2k,h1prk,h2prkT.


### 3.2. SH Nonlinearity

The nonlinearity outputv(k), represented by Figure 3 (DN4), can be developed in piecewise as follows: (9)vk=pZ1ifuk>Z1+aorZ1-a<uk<Z1+a,Δuk<0puk-aifZ2+a≤uk≤Z1+a,Δuk>0puk+aifZ2-a≤uk≤Z1-a,Δuk<0pZ2ifuk<Z2-aorZ2-a<uk<Z2+a,Δuk>0.The set of the above piecewise defining v(k) can be summarized in one equation as (10)vk=pukf1k-paf2k+pZ1f3k+pZ2f4k,with(11)f1k=α1-huk-Z1+a+hZ2+a-uk+1-α1-huk-Z1-a+hZ2-a-uk,f2k=α1-huk-Z1+a+hZ2+a-uk-1-α1-huk-Z1-a+hZ2-a-uk,f3k=αhuk-Z1+a+1-αhuk-Z1-a,f4k=αhZ2+a-uk+1-αhZ2-a-uk.The parameter and observation vectors can be defined, respectively, as follows:(12)gSHT=p,pa,pZ1,pZ2,φSHk=ukf1k,f2k,f3k,f4kT.


### 3.3. The Proposed Algorithm

In this section, we propose an iterative subspace identification algorithm in order to find the unknown parameters of UDN, order of the system, and the system matricesA,B,C,D.Substituting (8) in (4), we obtain (13)xk+1=Axk+B~φUDNk,yk=Cxk+D~φUDNk+εk,where B~=Bg1UDN,Bg2UDN,…,BgrUDN and D~=Dg1UDN,Dg2UDN,…,DgrUDN can be perceived as the LTI system state space realization. We can use MOESP algorithm [24] to estimate the order of the linear part and the component of state space matrices A^, B~^, C^, and D~^ from input-output data.The proposed algorithm minimizes the following cost function:(14)B^,D^,g^UDN=argminB,D,gUDNΘ~^BD-Θ~BD22,with (15)Θ~BD=Bg1UDN,…,BgrUDNDg1UDN,…,DgrUDN.The solution to this optimization problem (14) can be resolved using the Singular Value Decomposition (SVD) [25] of the matrix Θ^BD as(16)Θ~^BD=U1U2Σ100Σ2V1TV2T,where Σ1=diagσ1,…,σr, U1∈Rn×1, and V1∈Rr×1.The parameters of the different blocks are estimated iteratively using the proposed algorithm. Then, (8) are rewritten as follows: (17a)gsUDNT=p^s,p^sa^s,p^sb^s,p^sZ^1s,p^sZ^2s,e^1s,e^2s,(17b)φ^sUDNk=ukf^s1k,-f^s2k,f^s3k,h^s1k,h^s2k,h^s1prk,h^s2prkT,with(18)f^s1k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk+1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s+1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk+1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,f^s2k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk+1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s-1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk+1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,f^s3k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk-1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s+1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk-1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,h^s1k=αhuk-Z^s1+b^s+a^s+1-αhuk-Z^s1+b^s-a^s,h^s2k=αhZ^s2-b^s+a^s-uk+1-αhZ^s2-b^s-a^s-uk,h^s1prk=huk-b^s-a^s,h^s2prk=h-b^s+a^s-uk.The proposed iterative algorithm is summarized as follows.Step 1 (initialization).
g ^   0 U D N T and φ^(0)UDN(k) and the stopping criteria λ.Step 2.
s = s + 1; construct the information vector φ^(s)UDN(k).Step 3.
Estimate the system matricesA^(s),B~^(s),C^(s),D~^(s) and the model ordern using MOESP algorithm.Step 4.
Compute the SVD ofΘ~^(s)BD using (16).Step 5.
Computes the matrices parameterB^(s), D^(s), and g^(s)UDN as(19)B^sD^s=Us1Σs1,g^sUDN=Vs1,Θ^s=A^sB^sC^sD^s.Step 6.
Calculate Mean Square Error (MSE) using(20)MSEs=1N∑k=1Nyk-y^sk2,with(21)x^sk+1=A^sx^sk+B~^sφ^snewk,y^sk=C^sx^sk+D~^sφ^snewk,whereN is the number of data ykand y^(s)k denotes the real and the estimated outputs.
IfMSE(s)-MSE(s-1)<λ, then Θ^=Θ^(s) and g^UDN=g^(s)UDN. Otherwise, return to Step 2.


## 3.1. Unified Discontinuous Nonlinearity

Figure2 illustrates the v(k) evolution with respect to the input u(k) of the nonlinear block using UDN [9]. This latter contains an asymmetric piecewise function with a hysteresis, a saturation, a preload, a dead-zone, and so forth.Figure 2
Unified Discontinuous Nonlinearity.The UDN is characterized by the unknown parametersp, a, b, Z1, Z2, e1, and e2. Z1 and Z2 are scalar positive parameters that characterize the saturation, a≥0 and b≥0 define the width and the centre of the hysteresis as well as the dead-zone width, p is the hysteresis slope. e1 and e2 define the preload thresholds. The different choices of these parameters lead to the definition of eight nonlinearities which are listed in Table 1 and shown in Figure 3.Table 1
Different generated nonlinearities.
DN Choices of parameters UDN p ≠ 0, a≠0, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN1 p ≠ 0, a=b, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN2 p ≠ 0, a=b, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN3 p ≠ 0, a≠0, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN4 p ≠ 0, a≠0, b=0, Z1≠0, Z2≠0, e1=0, e2=0 DN5 p ≠ 0, a=0, b≠0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN6 p ≠ 0, a=0, b=0, Z1≠0, Z2≠0, e1≠0, e2≠0 DN7 p ≠ 0, a=0, b≠0, Z1≠0, Z2≠0, e1=0, e2=0 DN8 p ≠ 0, a=0, b=0, Z1≠0, Z2≠0, e1=0, e2=0Figure 3
Different cases of nonlinearities.Using Table1, the general form of v(k) with respect to u(k) and Δu(k) for the different cases of nonlinearities can be defined as(5)vk=pZ1+e1ifuk>Z1+b+a,Δuk<0orZ1+b-a<uk<Z1+b+a,Δuk<0puk-a-b+e1ifb+a≤uk≤Z1+b+a,Δuk>0puk-a+b+e2ifZ2-b+a≤uk≤-b+a,Δuk>0puk+a-b+e1ifb-a≤uk≤Z1+b-a,Δuk<0puk+a+b+e2ifZ2-b-a≤uk≤-b-a,Δuk<0pZ2+e2ifuk<Z2-b-a,Δuk<0orZ2-b<uk<Z2-b+a,Δuk>0e1ifb-a<uk<b+a,Δuk>0e1if-b-a<uk<-b+a,Δuk<00if-b+a≤uk≤b-a∀Δuk,with α=h[Δu(k)] and Δu(k)=u(k)-u(k-1).Using (2), (5) can be rewritten as(6)vk=pukf1k-paf2k+pbf3k+pZ1h1k+pZ2h2k+e1h1prk+e2h2prk=gUDNTφUDNk,with(7)f1k=α1-huk-Z1+b+a+hb+a-uk+1-hZ2-b+a-uk+huk--b+a+1-α1-huk-Z1+b-a+hb-a-uk+1-hZ2-b-a-uk+huk--b-a,f2k=α1-huk-Z1+b+a+hb+a-uk+1-hZ2-b+a-uk+huk--b+a-1-α1-huk-Z1+b-a+hb-a-uk+1-hZ2-b-a-uk+huk--b-a,f3k=α1-huk-Z1+b+a+hb+a-uk-1-hZ2-b+a-uk+huk--b+a+1-α1-huk-Z1+b-a+hb-a-uk-1-hZ2-b-a-uk+huk--b-a,h1k=αhuk-Z1+b+a+1-αhuk-Z1+b-a,h2k=αhZ2-b+a-uk+1-αhZ2-b-a-uk,h1prk=huk-b-a,h2prk=h-b+a-uk,hβ=1ifβ>00ifβ≤0,where hβ is a switching sequence.Using (6), (3) are rewritten as(8)gUDNT=p,pa,pb,pZ1,pZ2,e1,e2,φUDNk=ukf1k,-f2k,f3k,h1k,h2k,h1prk,h2prkT.


## 3.2. SH Nonlinearity

The nonlinearity outputv(k), represented by Figure 3 (DN4), can be developed in piecewise as follows: (9)vk=pZ1ifuk>Z1+aorZ1-a<uk<Z1+a,Δuk<0puk-aifZ2+a≤uk≤Z1+a,Δuk>0puk+aifZ2-a≤uk≤Z1-a,Δuk<0pZ2ifuk<Z2-aorZ2-a<uk<Z2+a,Δuk>0.The set of the above piecewise defining v(k) can be summarized in one equation as (10)vk=pukf1k-paf2k+pZ1f3k+pZ2f4k,with(11)f1k=α1-huk-Z1+a+hZ2+a-uk+1-α1-huk-Z1-a+hZ2-a-uk,f2k=α1-huk-Z1+a+hZ2+a-uk-1-α1-huk-Z1-a+hZ2-a-uk,f3k=αhuk-Z1+a+1-αhuk-Z1-a,f4k=αhZ2+a-uk+1-αhZ2-a-uk.The parameter and observation vectors can be defined, respectively, as follows:(12)gSHT=p,pa,pZ1,pZ2,φSHk=ukf1k,f2k,f3k,f4kT.


## 3.3. The Proposed Algorithm

In this section, we propose an iterative subspace identification algorithm in order to find the unknown parameters of UDN, order of the system, and the system matricesA,B,C,D.Substituting (8) in (4), we obtain (13)xk+1=Axk+B~φUDNk,yk=Cxk+D~φUDNk+εk,where B~=Bg1UDN,Bg2UDN,…,BgrUDN and D~=Dg1UDN,Dg2UDN,…,DgrUDN can be perceived as the LTI system state space realization. We can use MOESP algorithm [24] to estimate the order of the linear part and the component of state space matrices A^, B~^, C^, and D~^ from input-output data.The proposed algorithm minimizes the following cost function:(14)B^,D^,g^UDN=argminB,D,gUDNΘ~^BD-Θ~BD22,with (15)Θ~BD=Bg1UDN,…,BgrUDNDg1UDN,…,DgrUDN.The solution to this optimization problem (14) can be resolved using the Singular Value Decomposition (SVD) [25] of the matrix Θ^BD as(16)Θ~^BD=U1U2Σ100Σ2V1TV2T,where Σ1=diagσ1,…,σr, U1∈Rn×1, and V1∈Rr×1.The parameters of the different blocks are estimated iteratively using the proposed algorithm. Then, (8) are rewritten as follows: (17a)gsUDNT=p^s,p^sa^s,p^sb^s,p^sZ^1s,p^sZ^2s,e^1s,e^2s,(17b)φ^sUDNk=ukf^s1k,-f^s2k,f^s3k,h^s1k,h^s2k,h^s1prk,h^s2prkT,with(18)f^s1k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk+1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s+1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk+1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,f^s2k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk+1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s-1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk+1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,f^s3k=α1-huk-Z^s1+b^s+a^s+hb^s+a^s-uk-1-hZ^s2-b^s+a^s-uk+huk--b^s+a^s+1-α1-huk-Z^s1+b^s-a^s+hb^s-a^s-uk-1-hZ^s2-b^s-a^s-uk+huk--b^s-a^s,h^s1k=αhuk-Z^s1+b^s+a^s+1-αhuk-Z^s1+b^s-a^s,h^s2k=αhZ^s2-b^s+a^s-uk+1-αhZ^s2-b^s-a^s-uk,h^s1prk=huk-b^s-a^s,h^s2prk=h-b^s+a^s-uk.The proposed iterative algorithm is summarized as follows.Step 1 (initialization).
g ^   0 U D N T and φ^(0)UDN(k) and the stopping criteria λ.Step 2.
s = s + 1; construct the information vector φ^(s)UDN(k).Step 3.
Estimate the system matricesA^(s),B~^(s),C^(s),D~^(s) and the model ordern using MOESP algorithm.Step 4.
Compute the SVD ofΘ~^(s)BD using (16).Step 5.
Computes the matrices parameterB^(s), D^(s), and g^(s)UDN as(19)B^sD^s=Us1Σs1,g^sUDN=Vs1,Θ^s=A^sB^sC^sD^s.Step 6.
Calculate Mean Square Error (MSE) using(20)MSEs=1N∑k=1Nyk-y^sk2,with(21)x^sk+1=A^sx^sk+B~^sφ^snewk,y^sk=C^sx^sk+D~^sφ^snewk,whereN is the number of data ykand y^(s)k denotes the real and the estimated outputs.
IfMSE(s)-MSE(s-1)<λ, then Θ^=Θ^(s) and g^UDN=g^(s)UDN. Otherwise, return to Step 2.


## 4. Result and Discussions

In this section, we are going to examine the performance of the proposed method developed above.In the first step, we consider the unified form of the nonlinearity to determine the structure of the nonlinearity of the system. Once the nonlinearity is recognized, in the second step, we use it to identify the model with the exact form of the nonlinearity.


### 4.1. Example  1

We considered Hammerstein model where the parameters characterizing the DN4 nonlinearity are p=2, a=0.5, b=0, Z1=1.5, Z2=-1, e1=0, and e2=0.The LTI state space modelS(·) can be determined by looking at the system transfer function in factored form(22)Gz=z-1b1+b2z-11+a1z-1+a2z-2=z-11+0.5z-11-0.2z-1+0.35z-2.The iterative identification procedure was carried out for 2000 samples of uniformly distributed random inputs with u(k)≤3. ε(k) is a white noise with zero mean and variance σ2, which is added to the output system at different SNR levels (SNR = 10 and 8 dB), in order to demonstrate the effectiveness of the proposed algorithm in a noisy environment.The Monte Carlo experiments are carried out with 50 runs. The average estimated value (li-) and the standard deviation value (σli) obtained in the noisy case are(23)li-=1M∑m=1Mli,m,σli=1M∑m=1Mli,m-li-2,whereM is the number of replication and l^i=li-±σli.The structure of the nonlinearity is assumed unknown. For this reason, we begin with the first step to find the type of the nonlinearity using the UDN.The true parameters and the estimated parameters of nonlinear and linear blocks for the Hammerstein model are listed in Table2. Figure 4 represents the nonlinearity estimated by model UDN for different values of SNR.Table 2
Parameters estimation for the Hammerstein model with SNR = 10 and 8 dB.
SNR  Parameters for the nonlinear model Denominator Numerator MSE p a b Z 1 Z 2 e 1 e 2 a 1 a 2 b 1 b 2 Without noise True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.0004 l i - 2.0000 0.5001 0.0000 1.5001 −1.0000 0.0006 0.0007 −0.2001 0.3500 1.0013 0.4998 σ l i 0.0001 0.0003 0.0001 0.0005 0.0000 0.0004 0.0007 0.0001 0.0000 0.0002 0.0002  10 dB True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.0437 l i - 2.0027 0.4760 0.0326 1.4874 −1.0546 0.0538 −0.0368 −0.2006 0.3496 1.0029 0.4898 σ l i 0.0709 0.0618 0.0208 0.0264 0.0188 −0.0573 0.0189 0.0017 0.0010 0.0091 0.0073  8 dB True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.1596 l i - 2.0741 0.5341 0.1018 1.4210 −0.8515 −0.0715 −0.0290 −0.1980 0.3488 1.0385 0.4887 σ l i 0.0388 0.0106 0.0287 0.0941 0.0963 0.1068 0.0327 0.0016 0.0022 0.0247 0.0783Figure 4
DN4 estimated with UDN model.From Figure4, we can observe that the nonlinearity of this system corresponds to a saturation with hysteresis (DN4) in the absence of noise. As shown in Table 2 and Figure 4, the proposed method present good performance for the lower level of noise. With higher level of noise (SNR = 8 dB), the proposed algorithms exhibit poor performance. Despite that, we were able to identify the model structure of nonlinear block. Once the structure of nonlinear block is determined, we will resume again the identification steps using the exact model.We considered the parameters characterizing the SH nonlinearity:a=0.5, p=2, Z1=1.5, and Z2=-1. The true and the estimated parameters of nonlinear and linear blocks for the Hammerstein model with SH nonlinearity are summarized in Table 3.Table 3
Parameters estimation for the Hammerstein model with SNR = 10, 8, and 5 dB.
SNR  Parameters for the nonlinear model Denominator Numerator MSE p a Z 1 Z 2 a 1 a 2 b 0 b 1 Without noise True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 4.10 −6 l i - 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 0.9998 0.5001 σ l i 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0001 0.0001  10 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0129 l i - 2.0125 0.4955 1.4884 −0.9971 −0.2039 0.3468 0.9844 0.4939 σ l i 0.0149 0.0169 0.0154 0.0218 0.0067 0.0041 0.0248 0.0149  8 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0263 l i - 1.9881 0.4925 1.4877 −0.9964 −0.1964 0.3455 1.0198 0.5081 σ l i 0.0127 0.0304 0.0275 0.0278 0.0246 0.0213 0.0354 0.0358  5 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0396 l i - 2.0776 0.4824 1.4724 −0.9270 −0.2232 0.3328 0.9686 0.5198 σ l i 0.0144 0.0731 0.0308 0.0205 0.0455 0.0458 0.0153 0.0215The comparison simulation results with SH nonlinearity are given in Table3 with different SNR levels.As seen in Table3, our method with the parameters of the exact structure of nonlinear block gives the best MSE performance without and with noise. A similar analysis can be seen also in Figure 5.Figure 5
Nonlinearity SH estimated without SNR and with SNR = 10, 8 et 5 dB.


### 4.2. Example  2: Electropneumatic System

We use the proposed method to identify the parameters of an industrial process. The process is an electropneumatic system and shown in Figure6. The actual system consists of a carriage, a proportional valve 5/3, and a linear position sensor (linear potentiometer).Figure 6
Electropneumatic system.The proportional distributor 5/3, shown in Figure7, converts the analog electrical input signal into the openings of the cross sections adequate to the outputs. In midnominal voltage (5 V), the distributor takes its pneumatic position in which all the control edges are closed. In both directions (ascending and descending), the distributor witnesses a voltage-variation and this leads to the emergence of a hysteresis phenomenon.Figure 7
Distributor connection scheme.The linear potentiometer is used to measure the movement of the carriage. The latter converts the linear position into an output voltage which is proportional to that of the supply and the cursor position. This is performed by using the voltage divider principle.Figure8 shows the static characteristic of the proportional distributor given by the manufacturer corresponding to a single orifice 4. A distributor test is achieved by varying the voltage in both directions (ascending and descending) with opening orifice 4.Figure 8
Manufacturer test of static characteristic of the proportional distributor 5/3.We try to structurally identify the discontinuous static nonlinearity in the first step. To do this, we use the model of Hammerstein based on the UDN model. Two sets of inputs/outputs datasets (Figures9 and 10) are used to the system identification and validation. For each set, the first one is used as training data and the second one is used as test data.Figure 9
First set of measurement input after treatment.Figure 10
First set of measurement output after treatment.The identified parameters obtained by our method are listed in Table4. This table shows the average estimated value (li-) for the UDN and SH nonlinearities.Table 4
Estimates of the parameters for the nonlinear model.
Nonlinearity Parameters nonlinear block MSE p a b Z 1 Z 2 e 1 e 2 UDN 6.2850 0.6536 0.0701 0.5136 −0.5058 0.2785 −0.2512 0.0475 SH 8.9668 0.5740 — 0.4016 −0.3973 — — 0.0366Through Figure11, we can conclude that the UDN corresponds to the structure of DN4 regarding Figure 3. As shown in Figure 11, the DN4 gives a good characterization of the nonlinearity (the hysteresis phenomenon of the distributor 5/3) for the electropneumatic system. The MSEs of the estimates of the model using UDN is equal to 0.0475 and using SH is 0.0366 when the Hammerstein is taken with an exact nonlinearity.Figure 11
Nonlinearity with model UDN and SH structures.Using the second datasets (Figures9 and 10) for system validation, the performance of the proposed method is shown in Figure 12. We notice that the identification quality improves, especially at the level of nonlinearity. From Figure 12, it can be seen that our method provides a good approximation modelling accuracy.Figure 12
The real and the estimated outputs for the training and testing data sets.


## 4.1. Example  1

We considered Hammerstein model where the parameters characterizing the DN4 nonlinearity are p=2, a=0.5, b=0, Z1=1.5, Z2=-1, e1=0, and e2=0.The LTI state space modelS(·) can be determined by looking at the system transfer function in factored form(22)Gz=z-1b1+b2z-11+a1z-1+a2z-2=z-11+0.5z-11-0.2z-1+0.35z-2.The iterative identification procedure was carried out for 2000 samples of uniformly distributed random inputs with u(k)≤3. ε(k) is a white noise with zero mean and variance σ2, which is added to the output system at different SNR levels (SNR = 10 and 8 dB), in order to demonstrate the effectiveness of the proposed algorithm in a noisy environment.The Monte Carlo experiments are carried out with 50 runs. The average estimated value (li-) and the standard deviation value (σli) obtained in the noisy case are(23)li-=1M∑m=1Mli,m,σli=1M∑m=1Mli,m-li-2,whereM is the number of replication and l^i=li-±σli.The structure of the nonlinearity is assumed unknown. For this reason, we begin with the first step to find the type of the nonlinearity using the UDN.The true parameters and the estimated parameters of nonlinear and linear blocks for the Hammerstein model are listed in Table2. Figure 4 represents the nonlinearity estimated by model UDN for different values of SNR.Table 2
Parameters estimation for the Hammerstein model with SNR = 10 and 8 dB.
SNR  Parameters for the nonlinear model Denominator Numerator MSE p a b Z 1 Z 2 e 1 e 2 a 1 a 2 b 1 b 2 Without noise True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.0004 l i - 2.0000 0.5001 0.0000 1.5001 −1.0000 0.0006 0.0007 −0.2001 0.3500 1.0013 0.4998 σ l i 0.0001 0.0003 0.0001 0.0005 0.0000 0.0004 0.0007 0.0001 0.0000 0.0002 0.0002  10 dB True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.0437 l i - 2.0027 0.4760 0.0326 1.4874 −1.0546 0.0538 −0.0368 −0.2006 0.3496 1.0029 0.4898 σ l i 0.0709 0.0618 0.0208 0.0264 0.0188 −0.0573 0.0189 0.0017 0.0010 0.0091 0.0073  8 dB True 2.0000 0.5000 0.0000 1.5000 −1.0000 0.0000 0.0000 −0.2000 0.3500 1.0000 0.5000 0.1596 l i - 2.0741 0.5341 0.1018 1.4210 −0.8515 −0.0715 −0.0290 −0.1980 0.3488 1.0385 0.4887 σ l i 0.0388 0.0106 0.0287 0.0941 0.0963 0.1068 0.0327 0.0016 0.0022 0.0247 0.0783Figure 4
DN4 estimated with UDN model.From Figure4, we can observe that the nonlinearity of this system corresponds to a saturation with hysteresis (DN4) in the absence of noise. As shown in Table 2 and Figure 4, the proposed method present good performance for the lower level of noise. With higher level of noise (SNR = 8 dB), the proposed algorithms exhibit poor performance. Despite that, we were able to identify the model structure of nonlinear block. Once the structure of nonlinear block is determined, we will resume again the identification steps using the exact model.We considered the parameters characterizing the SH nonlinearity:a=0.5, p=2, Z1=1.5, and Z2=-1. The true and the estimated parameters of nonlinear and linear blocks for the Hammerstein model with SH nonlinearity are summarized in Table 3.Table 3
Parameters estimation for the Hammerstein model with SNR = 10, 8, and 5 dB.
SNR  Parameters for the nonlinear model Denominator Numerator MSE p a Z 1 Z 2 a 1 a 2 b 0 b 1 Without noise True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 4.10 −6 l i - 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 0.9998 0.5001 σ l i 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0001 0.0001  10 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0129 l i - 2.0125 0.4955 1.4884 −0.9971 −0.2039 0.3468 0.9844 0.4939 σ l i 0.0149 0.0169 0.0154 0.0218 0.0067 0.0041 0.0248 0.0149  8 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0263 l i - 1.9881 0.4925 1.4877 −0.9964 −0.1964 0.3455 1.0198 0.5081 σ l i 0.0127 0.0304 0.0275 0.0278 0.0246 0.0213 0.0354 0.0358  5 dB True 2.0000 0.5000 1.5000 −1.0000 −0.2000 0.3500 1.0000 0.5000 0.0396 l i - 2.0776 0.4824 1.4724 −0.9270 −0.2232 0.3328 0.9686 0.5198 σ l i 0.0144 0.0731 0.0308 0.0205 0.0455 0.0458 0.0153 0.0215The comparison simulation results with SH nonlinearity are given in Table3 with different SNR levels.As seen in Table3, our method with the parameters of the exact structure of nonlinear block gives the best MSE performance without and with noise. A similar analysis can be seen also in Figure 5.Figure 5
Nonlinearity SH estimated without SNR and with SNR = 10, 8 et 5 dB.


## 4.2. Example  2: Electropneumatic System

We use the proposed method to identify the parameters of an industrial process. The process is an electropneumatic system and shown in Figure6. The actual system consists of a carriage, a proportional valve 5/3, and a linear position sensor (linear potentiometer).Figure 6
Electropneumatic system.The proportional distributor 5/3, shown in Figure7, converts the analog electrical input signal into the openings of the cross sections adequate to the outputs. In midnominal voltage (5 V), the distributor takes its pneumatic position in which all the control edges are closed. In both directions (ascending and descending), the distributor witnesses a voltage-variation and this leads to the emergence of a hysteresis phenomenon.Figure 7
Distributor connection scheme.The linear potentiometer is used to measure the movement of the carriage. The latter converts the linear position into an output voltage which is proportional to that of the supply and the cursor position. This is performed by using the voltage divider principle.Figure8 shows the static characteristic of the proportional distributor given by the manufacturer corresponding to a single orifice 4. A distributor test is achieved by varying the voltage in both directions (ascending and descending) with opening orifice 4.Figure 8
Manufacturer test of static characteristic of the proportional distributor 5/3.We try to structurally identify the discontinuous static nonlinearity in the first step. To do this, we use the model of Hammerstein based on the UDN model. Two sets of inputs/outputs datasets (Figures9 and 10) are used to the system identification and validation. For each set, the first one is used as training data and the second one is used as test data.Figure 9
First set of measurement input after treatment.Figure 10
First set of measurement output after treatment.The identified parameters obtained by our method are listed in Table4. This table shows the average estimated value (li-) for the UDN and SH nonlinearities.Table 4
Estimates of the parameters for the nonlinear model.
Nonlinearity Parameters nonlinear block MSE p a b Z 1 Z 2 e 1 e 2 UDN 6.2850 0.6536 0.0701 0.5136 −0.5058 0.2785 −0.2512 0.0475 SH 8.9668 0.5740 — 0.4016 −0.3973 — — 0.0366Through Figure11, we can conclude that the UDN corresponds to the structure of DN4 regarding Figure 3. As shown in Figure 11, the DN4 gives a good characterization of the nonlinearity (the hysteresis phenomenon of the distributor 5/3) for the electropneumatic system. The MSEs of the estimates of the model using UDN is equal to 0.0475 and using SH is 0.0366 when the Hammerstein is taken with an exact nonlinearity.Figure 11
Nonlinearity with model UDN and SH structures.Using the second datasets (Figures9 and 10) for system validation, the performance of the proposed method is shown in Figure 12. We notice that the identification quality improves, especially at the level of nonlinearity. From Figure 12, it can be seen that our method provides a good approximation modelling accuracy.Figure 12
The real and the estimated outputs for the training and testing data sets.


## 5. Conclusion

In this paper, a new iterative subspace algorithm for the Hammerstein model identification is presented. The application of the proposed method improves a major advantage which is the simultaneous identification of the parameters of the linear and nonlinear blocks. Indeed, these two blocks are assumed unknown. The use of UDN model with Hammerstein model can be identified the adapted structure for a given system with and without noise. Furthermore, an iterative subspace algorithm is used to finally tune the parameters of model.Two examples are presented to verify the effectiveness of the proposed approach. As a result, the proposed method is reliable to describe approximately the real nonlinearity of such as systems and can be applied to other various industrial systems to a level of well-defined noise.



---

*Source: 1794921-2016-11-30.xml*