Soft Comput DOI 10.1007/s00500-017-2749-6

METHODOLOGIES AND APPLICATION

Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain approach Mohammad-Reza Rahmani1 · Mohammad Farrokhi1

© Springer-Verlag GmbH Germany 2017

Abstract In this paper, modeling and identification of nonlinear dynamic systems using neuro-fractional Hammerstein model are considered. The proposed model consists of the neural networks (NNs) as the nonlinear subsystem and the fractional-order state space (FSS) as the linear subsystem. The identification procedure consists of a hybrid frequency/time-domain approach based on the input–output data acquired from the system. First in the frequency domain, the fractional order and fractional degree of the FSS subsystem are determined offline using an iterative linear optimization algorithm. Then, in the time domain, the state-space matrices of the FSS as well as parameters of the NN are estimated using Lyapunov stability theory. Moreover, in order to use only the input–output data from the system, a fractional-order linear observer based on auxiliary model idea is utilized to estimate the system states. The convergence and stability analysis of the proposed method are provided. Simulating and experimental examples show superior performance of the proposed method as compared with the Hammerstein models reported in the literature. Keywords Hammerstein model · Fractional-order differential equations · Neural networks · Frequency-domain identification

Communicated by V. Loia.

B

Mohammad Farrokhi [email protected] Mohammad-Reza Rahmani [email protected]

1

School of Electrical Engineering, Iran University of Science and Technology, Tehran 16846-13114, Iran

1 Introduction In the real world, all physical systems are inherently nonlinear, which can degrade linear system identification methods (Pelt and Bernstein 2010). In the past few decades, identification of nonlinear block-oriented systems such as Hammerstein and Wiener systems as well as their combinations has been a subject of growing research interest (Yu et al. 2013). Among these block-oriented systems, the Hammerstein system has attracted the most attention because of its simplicity and acceptable modeling (Zhao et al. 2014). Generally, a Hammerstein system consists of a nonlinear memoryless subsystem (NMS) followed by a linear dynamic subsystem (LDS). It has been shown that the Hammerstein system can describe a wide variety of nonlinear systems such as chemical processes, electrical drives, and thermal systems (Guo et al 2012). State-space systems have an important role in modeling the dynamics of nonlinear systems. Although many approaches have been presented for linear state-space systems, the identification for nonlinear state-space systems is still developing (Wang and Ding 2015). The major difficulty of the state-space modeling is the simultaneous estimation of the states as well as the parameters. Wang and Ding (2015) and Ding et al. (2016) have adopted a state observer based on the auxiliary model identification idea to overcome this problem. Neural networks (NNs) have been widely employed for the identification of nonlinear systems (Ibnkahla 2012). Two different approaches have been proposed in the literature for identification of Hammerstein systems using the NNs. In the first approach, the NNs are used to represent the NMS (Chen et al. 2014; Cui et al. 2014; Mkadem and Boumaiza 2009). Hong et al. (2014) have reviewed the advantages of using the NNs over the polynomials. On the other hand, in the second

123

M. Rahmani, M. Farrokhi

approach, the NNs are employed to identify the systems that are in Hammerstein form (Jia et al. 2005; Ren and Lv 2011). Recently, there has been increasing interest among researchers on the fractional-order differential equations due to their capacity to accurately mimic the dynamic behavior of the real-world systems; among them, electrical devices (Caponetto et al. 2013; Galvão et al. 2013) and thermal systems (Benchellal et al. 2006; Gabano and Poinot 2011a, b) can be mentioned. The fractional-order Hammerstein system is a kind of Hammerstein-type systems with the LDS being of the fractional-order type. To the best of the authors’ knowledge, little work has been reported on identification of fractionalorder Hammerstein systems (Aoun et al. 2006; Liao et al. 2012; Zhao et al. 2014b), while there is no identification study in the literature on neuro-fractional Hammerstein systems (NFHSs). In Aoun et al. (2006), it is assumed that the fractional order is known a priori. In Liao et al. (2012) and in Zhao et al. (2014b), iterative time-domain algorithms are proposed to estimate the fractional order. However, the number of the LDS states, which is called here the fractional degree, should be known a priori. In Liao et al. (2012), a nonlinear optimization problem in which the NMS and LDS are identified in the time domain via subspace identification technique must be solved to estimate the fractional order. This method is prone to local minima issues due to the use of nonlinear optimization. Moreover, since in every step of the nonlinear optimization, the NMS and LDS are estimated via subspace identification, this algorithm may suffer from long computational time. In Zhao et al. (2014b), two methods have been proposed to estimate the fractional order: (1) a one-dimensional grid search is suggested. In this method, at every point in the grid of the fractional order, the parameters of the NMS and the LDS are estimated iteratively in the time domain based on the least-square algorithms. (2) The fractional order, the NMS, and the LDS are identified iteratively, wherein the fractional-order estimation is based on the so-called P-type order learning law. However, the major disadvantage of this identification method is that the derivatives of the input and output signals appear in the regression model. Moreover, in this reference it is claimed that the convergence of the fractional-order identification is guaranteed for small learning rate, which may make the algorithm very time-consuming. Besides, the overall stability of this iterative algorithm is not shown. In this manuscript, the NNs and the fractional-order linear systems are employed in the Hammerstein model for the identification of nonlinear systems. In the proposed method, the NN (i.e., the NMS) is of the radial-basis-function (RBF) type, while the LDS is of the commensurate fractional-order kind. It incorporates the advantages of the RBF in static continuous function approximation and the benefits of FSS on system identification. It is noteworthy that any other NN

123

can be used in replacement of the RBF, provided that the NN is linear in unknown parameters. The advantages of the proposed Hammerstein model are as follows: (1) There is no need for the nonlinear system to be in the Hammerstein form. (2) The identification algorithm is based on the input– output data acquired from the system. (3) There is no need for persistent excitation of the nonlinear system, while in many reported methods (e.g., Ding et al. 2013, 2014) persistent excitation is a requirement. In addition, in almost all papers (e.g., Chen et al. 2014; Cui et al. 2014; Hong et al. 2014) the LDS is designed as a discrete-time filter, while in this manuscript the LDS is a continuous system, which can provide more realistic description of the real-life continuous systems (Garnier and Young 2014). In this paper, in comparison with the aforementioned papers (e.g., Liao et al. 2012; Zhao et al. 2014b), the fractional order and the degree of the LDS are identified in the frequency domain using genetic algorithm (GA). It is noteworthy that there exist many alternative evolutionary and biologically inspired algorithms, i.e., ant colony optimization (ACO), bee colony algorithm (BCA), artificial fish swarm optimization (AFSO) among which particle swarm optimization (PSO) has been used widely (Lorenzo et al. 2017). One of the major advantages of the proposed method is that it drastically reduces the computational time. In order to avoid solving a nonlinear optimization problem and not getting trapped in local minima, a two-dimensional stochastic search (i.e., GA) for the fractional order and the degree of LDS is employed. Then, for each point on the grid, a linear optimization is performed to estimate the NMS and LDS. Moreover, a frequency-domain identification scheme is proposed that does not need the derivatives of the input and output signals. Afterward, the parameters of the NMS and the coefficients of LDS are determined online using the Lyapunov stability theory. This paper is organized as follows: In Sect. 2, the model structure is introduced and the problem statement is presented. In Sect. 3, the proposed frequency-/time-domain identification algorithm is given. Some simulation results are presented in Sect. 4. Finally, Sect. 5 concludes the paper.

2 Problem statement Consider the following NFHS: ⎧ N    ⎪ ⎪ wi exp −σi2 (u(t) − μi )2 = w T ϕ (u(t)) , ⎨ v(t) = i=1

⎪ D α x(t) = Ax(t) + bv(t), ⎪

⎩ y(t) = 0 · · · 0 1 x(t),

Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain…

⎡ ⎤ 0 x1 (t) ⎢1 ⎢ ⎢ x2 (t) ⎥ ⎢ ⎢ ⎥ x(t) = ⎢ . ⎥ , A = ⎢ 0 ⎢. ⎣ .. ⎦ ⎣ .. xn (t) 0 ⎡ ⎤ b0 ⎢ b1 ⎥ ⎢ ⎥ b=⎢ . ⎥ ⎣ .. ⎦ ⎡

bn−1 D x(t) = D α x1 (t) α

···

0 ··· 0 ··· 1 0 .. .. . . 0 0

0 0 ··· .. . 1

D α xn (t) ,

−a0 −a1 −a2 .. .

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦

−an−1

3.1 Estimation of LDS structural parameters

(1)

where u(t) is the system input, ϕ = [ φ1 . . . φ N ]T are

T known functions of the input u(t), w = w1 . . . w N are weights of NN, x(t) is state vector, y(t) is the output of the system, v(t) is unmeasurable intermediate variable, A and b are state-space matrix/vector, respectively, and the symbol D α denotes the Caputo fractional-order derivative defined as (Li et al. 2009; Podlubny 1999) D α f (t) =

 t 1 ([α] − α + 1) 0  [α]+1  d 1 f (τ ) dτ , (t − τ )α−[α] dτ [α]+1

(2)

where[α] denotes the integer part of α and (·) represents the Gamma function. Using the Laplace transform of an αorder derivative of signal f (t) with zero initial conditions defined as L{D α f (t)} = s α F(s), the strictly proper transfer function of the LDS can be represented as (Djamah et al. 2009) G(s) =

n−1 iα Y (s) i=0 bi s = n−1 iα . V (s) s nα + i=0 ai s

online. The first step is performed via a linear optimization algorithm in the frequency domain, while the second step is carried out using the Lyapunov stability theory in the time domain.

(3)

The system identification problem considered in this paper is as follows: Using the input–output measurements {u(t), y(t)} of system (1), the goal is to find the fractional order α, the degree n, a realization of the state-space matrices (A, b) for the LDS, and the weights of the NMS (w). Then, the proposed NFHS modeling and identification algorithm are developed to identify general nonlinear dynamic systems.

3 Parameters estimation The proposed identification method of the NFHS in (1) is accomplished in two steps: (1) estimation of the LDS structural parameters (i.e., the fractional order α and the degree n) in offline and (2) identifying the coefficients of the transfer function of the LDS and the weights of the RBF in NSM in

The parameter estimation of the LDS structure is performed in two steps: nonparametric estimation and iterative parametric estimation. First, by mapping the sampled input–output data to the frequency domain, the nonparametric frequency response of the LDS is obtained. Then, an iterative linear optimization algorithm is proposed to define an appropriate commensurate transfer function in the mean-square sense. The nonparametric frequency response function of the NFHS is obtained as (Chen and Wang 2009) ∗ ˆ jωk ) = Y (k)U (k) , G( U (k)U ∗ (k)

(4)

√ L−1 where ωk = 2π k/Lts , j = −1, {U (k), Y (k)}k=0 stands L , for the discrete Fourier transform of {u(its ), y(its )}i=1 ∗ respectively, and U (k) denotes the conjugate of U (k). The sampling time at which the input–output data are sampled is represented by ts . Moreover, L is the total number of the input–output data. The basic idea behind this formulation is the best linear approximation (BLA), which is discussed in detail in Pintelon and Schoukens (2012). Remark 1 In order to obtain a satisfactory performance, it is preferred to apply an input signal with the following properties: 1. The input signal should have a frequency spectrum, which covers the frequency interval of interest such that the frequency components of the LDS are captured. In other words, the input–output date should be rich in the frequency response. 2. The amplitude of the input signal must be small enough such that it does not excite the nonlinear modes of the system. Next, the LDS structural parameters are determined via a constrained stochastic search, namely GA, by minimizing the cost function given in Theorem 1. The stability of the estimated transfer function is based on the constraint. In other words, it must be ensured that the denominator of (3) in terms of s α has no zeroes in the right half of the s-plane. This constraint is a consequence of Assumption A1 that is stated in the next section. The entire search region defined as  is a two-dimensional space    of the GA α, ˆ nˆ ∈ α, ˆ nˆ  0 < αˆ ≤ 1, 1 ≤ nˆ ≤ n max , where n max is

123

M. Rahmani, M. Farrokhi

an arbitrary integer. However in this paper, n max = 4 is ˆ jωk ) with a relatively selected in order to approximate G( low-order transfer function. In this way, for a moderate value of the fractional order, i.e., α = 0.5, the maximum order that appears in the estimated parametric transfer function is n max α = 2. Hence, the resulted transfer function will be of the well-known second-order type. The linear optimization algorithm can be obtained by the following re-parameterization of G(s) in (3): G( jω) = w= 

b = a = z=

n−1  iα M( jω) oT b i=0 b ( jω) = = , n i  N ( jω) 1 + z T a 1 + i=1 ai ( jω)iα

T ∈ Cn , 1 · · · ( jω)(n−1)α

 T  b0 · · · bn−1 ∈ Rn , 

T a1 · · · an ∈ Rn and

T ∈ Cn . (5) ( jω)α · · · ( jω)nα

The following theorem provides the linear optimization algorithm that evaluates the cost function for various values of (α, ˆ n). ˆ Theorem 1 Consider the fractional-order system defined in (5) with a priori known values of α and n. Let θˆ be an estimate T  that is given by of the parameter vector θ = b T a T −1  T FT (g I + g R ) θˆ = F F

(6)

Proof Equation (5) is a complex-valued equation, which can be decomposed into the real and imaginary parts as  G R ( jωl ) = oTR ( jωl )b − G R ( jωl )zTR ( jωl )  −G I ( jωl )zTI ( jωl ) a ,  G I ( jωl ) = oTI ( jωl )b − G R ( jωl )zTI ( jωl )  +G I ( jωl )zTR ( jωl ) a .

(8)

(9)

The addition of the last two equations can be written in matrix form as g = F θ,

(10)

where g and F are the same as in (6). To find the estimated value of θ in (10), the least-squares optimization method can be used here: T −1 T θˆ = (F F) F g.

(11)

This concludes the proof.



The proposed frequency-domain identification can be summarized as follows: ˆ jω) Step 1. Obtain a nonparametric frequency response G( defined in (4). Step 2. Perform the constrained  GA on a two-dimensional    space α, ˆ nˆ ∈ α, ˆ nˆ  0 < αˆ ≤ 1, 1 ≤ nˆ ≤ n max , with the cost function given in (7).

where

T g = G( jω1 ) G( jω2 ) · · · G( jω L ) ,    ⎤ ⎡ T o R ( jω1 ) + oTI ( jω1 ) G I ( jω1 ) zTI ( jω1 ) − zTR ( jω1 ) − G R ( jω1 ) zTI ( jω1 ) + zTR ( jω1 ) ⎢ T    ⎥ ⎢ o ( jω2 ) + oT ( jω2 ) G I ( jω2 ) zT ( jω2 ) − zT ( jω2 ) − G R ( jω2 ) zT ( jω2 ) + zT ( jω2 ) ⎥ I I R I R ⎢ R ⎥ ⎥ F=⎢ ⎢ ⎥ .. .. ⎢ ⎥ . . ⎣ ⎦     T T T T T T o R ( jω L ) + o I ( jω L ) G I ( jω L ) z I ( jω L ) − z R ( jω L ) − G R ( jω L ) z I ( jω L ) + z R ( jω L ) in which subscripts R and I denote the real and imaginary parts, respectively, L ≥ 2n is the number of data and ωl (l = 1, . . . , L) are selected such that FT F is nonsingular. Then, the cost function   T  g − Fθˆ ε = g − Fθˆ

(7)

3.2 Lyapunov-based parameter estimation of LDS and NMS After obtaining αopt and n opt , the updating laws for the statespace matrices (A, b) of the LDS and the weights of the NMS (w) are derived based on the Lyapunov stability theory. To do so, the following assumptions must be made first:

is minimized. A1. Matrix A, defined in (1), is stable.

123

Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain…

A2. αopt = α and n opt = n, i.e., there are no estimation errors for α and n. A3. In order to get a consistent and unique estimation, without loss of generality, the coefficient b0 in (1) is set to one. This assumption is common in the related literature (e.g., see Chen et al. 2014; Chen 2014; Wang et al. 2013). A4. A(t), b(t), and w (t) in (1) are time invariant, i.e., D α A(t), b(t), w(t) = 0. A5. The system states (x) are measurable.

learning laws guarantee that the state prediction errors x˜ (t) and the parameter estimation errors tend to zero as t → ∞, i.e., ˜ ˜ ˜ b(t), w(t) = 0, lim x˜ (t), A(t),

t→∞

(17)

˜ ˆ ˜ ˜ ˆ where A(t) = A(t) − A(t), w(t) = w(t) − w(t) and b(t) = ˆ b(t) − b(t).

Consequently, the following model is considered:

Proof From (1) and (12), the error dynamic equation can be obtained as

⎧   N  ⎪ ⎪ ˆ T (t)ϕ (u(t)) , ˆ = wˆ i (t) exp −σi2 (u(t) − μi )2 = w ⎨ v(t)

D α x˜ = D α x − D α xˆ T T ˜ x + Ax˜ + b˜ w ˆ ϕ + bˆ w ˜ T ϕ. ˜ ϕ + b˜ w = Aˆ

i=1

α ˆ ˆ v(t), ⎪ ⎪

ˆ ⎩ D xˆ (t) = A(t)ˆx(t) + b(t) yˆ (t) = 0 · · · 0 1 xˆ (t), ⎤ ⎡ xˆ 1 (t) ⎢ xˆ2 (t) ⎥ ⎥ ⎢ xˆ (t) = ⎢ . ⎥ , ⎣ .. ⎦ xˆn (t) ⎡ 0 0 ··· 0 −aˆ 0 (t) ⎢1 0 ··· 0 −aˆ 1 (t) ⎢ ⎢ 0 1 0 · · · − aˆ 2 (t) ˆ A(t) =⎢ ⎢. . . .. . .. .. ⎣ .. .. .

⎡

0

0

0 ⎤

1

Consider the following Lyapunov function:   1  ˜ T −1 ˜  1 ˜ w ˜ = x˜ T P1 x˜ + tr A P2 A V x˜ , A, 2 2 1 T −1 1 ˜ P4 w, ˜ + b˜ T P3−1 b˜ + w 2 2

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦

Dα V ≤ (12)

ˆ b, ˆ and yˆ are the estimates of their correˆ xˆ , A, wherev, ˆ w, sponding functions. Lemma 1 Consider x(t) ∈ R to be a derivable and continuous function. Then 1 α 2 α t D x (t) ≤ x(t)t0 Dt x(t), ∀α ∈ (0, 1] , t ≥ t0 . 20 t Proof See Aguila-Camacho et al. (2014).

(13) 

Theorem 2 Consider the NFHS in (1) and the identification ˆ bˆ and w ˆ model in (12). The following updating laws for A, under Assumptions A1–A5 guarantee stability and convergence of the proposed identification method: ˆ D A(t) = P2 P1 x˜ (t)ˆx (t), αˆ ˆ D b(t) = P3 P1 x˜ (t)ϕT (t)w(t),

(14)

ˆ = P4 ϕ(t)bT (t)P1 x˜ (t). D α w(t)

(16)

T

(19)

where tr(A) denotes the trace of A. Using Lemma 1, the α-order derivative of (19) can be written as

−aˆ n−1 (t)

1 ⎢ bˆ1 (t) ⎥ ⎢ ⎥ ˆ b(t) =⎢ ⎥, .. ⎣ ⎦ . ˆbn−1 (t)

α

(18)

  1 1  α T D x˜ P1 x˜ + x˜ T P1 D α x˜ 2  2     ˜ T P−1 A ˜ + D α b˜ T P−1 b˜ +tr D α A 2 3   ˜ T P4−1 w. ˜ + Dα w

(20)

According to the Lyapunov stability theory, condition D α V < 0 must be fulfilled to guarantee asymptotic convergence. This can be insured when the right-hand side of (20) is nega˜ = −D α A, ˆ tive. Based on Assumption A4, we have D α A ˆ and D α w ˜ = −D α w. ˆ Consequently, by subD α b˜ = −D α b, stituting learning laws (14)–(16) and error dynamic equation (18) into the right side of (20) and using some trace properties and disregarding the higher-order error terms of the form T b˜ w ˜ in (18), the following inequality holds:      ˜ T P−1 D α A ˆ ˜ T P1 x˜ xˆ T − tr A D α V ≤ −˜x T Q˜x + tr A 2   −1 T T T α ˆ − b˜ P +b˜ P1 x˜ ϕ w D bˆ 3

  ˜ T P4−1 D α w ˆ < 0, ˜ T ϕbˆ P1 x˜ − w +w T

(21)

(15)

where Pi (i = 1, . . . , 4) are arbitrary positive definite weight matrices and x˜ (t) = x(t) − xˆ (t). In other words, the above

where Q = (−AT P1 − P1 A)/2. It is noteworthy that for the stable matrix A and any positive definite matrix Q, there exists a positive definite matrix P1 that satisfies (AT P1 + P1 A)/2 = −Q (Fu et al. 2013). Using the fol-

123

M. Rahmani, M. Farrokhi

⎡

lowing equations: 

 ˜ T P−1 D α A ˆ =A ˜ T P1 x˜ xˆ T , A 2   ˆ b˜ T P3−1 D α bˆ = b˜ T P1 x˜ ϕT w,

  T ˜ T P4−1 D α w ˆ =w ˜ T ϕbˆ P1 x, w ˜

(22)

(24)

123

⎡ ⎤ ⎤ ··· 0 −aˆ 0 (t) ⎢ −aˆ 1 (t) ⎥ ··· 0 ⎥ ⎢ ⎥ ⎥ , aˆ ba (t) = ⎢ ⎥, .. 0 ···⎥ ⎣ ⎦ ⎦ . .. .. . . −aˆ n−2 (t)

· · · 1 , aˆ aa (t) = −aˆ n−1 (t),

(29)

where ε is the desired accuracy.

(26)

where κ(t) is the suitable time-varying gain vector of the observer and 0 0 1 .. .

Moreover, the following stopping criterion for the parameter updates is defined to improve the identification efficiency:

i=1

xo (t) = η(t) + κ(t)y(t), D α η(t) = (Abb − κ(t)aab ) η(t)   + (Abb − κ(t)aab ) κ(t) + aˆ ba (t) − κ(t)aˆ aa (t) y(t)   ˆ T (t)ϕ, + bb (t) − κ(t)bˆa (t) w (27)

0 ⎢1 ⎢ Abb = ⎢ 0 ⎣ .. . aab = 0 0

(28)

(25)

The structure of the reduced-order state observer is as follows

⎡

⎥ ⎥ ˆ ⎥ , ba (t) = bˆn−1 (t). ⎦

N    wˆ i (t) − wˆ i (t − 1) < ε,

Hence, if (22)–(24) hold, then the convergence properties ˆ α b, ˆ and D α w ˆ are guaranteed. Solving (22)–(24) for D α A,D leads to the updating laws (14), (15), and (16), respectively.  Usually, the system states (x) are unmeasurable, and hence, the aforementioned updating laws in (14)–(16) are not applicable in this situation. To make the proposed method more practical (i.e., using only the input–output data from the system for identification purposes), a fractional-order linear observer based on the idea of the auxiliary model is utilized here. The auxiliary model, where the unmeasurable variables are replaced with their estimates, is widely used in various Hammerstein system identification strategies (Hu and Ding 2014; Hu et al. 2013). It can be seen from (1) that xn = y. Thus, x(t) is replaced

T by xo (t) y (t) where xo (t) is achieved by using the reduced-order state observer for the following state-space model:   ˆ ˆ ˆ T (t)ϕ . + b(t) w D α x(t) = A(t)x(t)

⎤

bˆn−2 (t)

(23)

equation (21) reduces to D α V ≤ −˜x T Q˜x < 0.

⎢ ⎢ bˆ b (t) = ⎢ ⎣

1 bˆ1 (t) .. .

ˆ Remark 2 To initialize the algorithm, w(0) is taken as arbitrary numbers, and according to Assumption A3 and due to the fact that observable canonical state space of LDS is used, ˆ ˆ b(0)and A(0) are set to be ⎡

0 ⎢1 ⎢ b ˆ ⎢ ˆ b(0) =  , A(0) = ⎢0 ⎢. b0 ⎣ ..

0 0 1 .. .

0

0

··· 0 ··· 0 0 ··· .. . . . . 0 1

−1/an −a1 /an −a2 /an .. .

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦

 /an −an−1

(30) where b and ai are the frequency-domain estimation of θ =  T   for αopt , n opt . b T a T The whole procedure for the proposed NFHS identification can be summarized as follows:   1. Apply the frequency-domain algorithm to estimate α, ˆ nˆ   for the optimal LDS structural parameters αopt , n opt using Sect. 3.1.   2. Using the LDS structural parameters αopt , n opt , apply the time-domain algorithm described in Sect. 3.2 to estimate A, b and w. 3. Return to step (2) until the stopping criterion (29) is reached.

4 Simulating examples The effectiveness of the proposed NFHS identification algorithm is demonstrated in two examples. The stability and convergence properties of the presented algorithm are illustrated on a simulated NFHS. Then, a comparative study for the proposed NFH modeling of an experimental heating system, whose data set is taken from (De moor 2015), is given. Example 1 (identification of the proposed NFHS). Consider the following NFHS:

Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain… 



2 (u(t))2



+ 3e



−0.42 (u(t)−5)2



(31)

1.2 1 0.8 0.6 0.4 0.2

3 Real

1

⎤  0.9898 0 ˆ ⎣ ⎦ ˆ f = 1.9504 , A f = w 1 2.9248

a 20

30

bˆ f = 0.0508,   w − w ˆ f = 0.0242, δb = 0.0153, δA = 0.0112 δw =

w (33) which show that the time-domain algorithm identifies the parameters of the RBF accurately and improves the accuracy of the frequency-domain estimation of the LDS parameters. Furthermore, to verify the validation of this model, a random input with uniform distribution on the interval [−10, 10] sampled at 1 sample per second is applied to the system and to the model with estimated parameters given by (33). Figure 2 demonstrates the effectiveness of the proposed identification algorithm. As this figure shows, the identified model mimics the input–output behavior of the system in (31) with very good accuracy. The performance of the proposed algorithm

10

20

30

40

40

0

10

20

30

40

0

10

20

30

40

18

20

Time (seconds) 0.055 1

0.05

b

a

1

-60 -62 0

10

20

30

0.045

40

Fig. 1 Parameter estimation results of Example 1 x 10

-3

6 4 2

2

0

4

6

8

10

12

14

16

System output Model output

-5

10 5 0 -5

0

2

4

6

8

10

12

14

16

18

20

Time (seconds)



−494.3996 , −59.4770

0 -495 -500 -505 -510 -515

0

3

w

10

x 10

Figure 1 shows the estimation results. According to this figˆ f , bˆ f ˆ f,A ure, the final values of the estimated parameters w and their corresponding errors δw,A,b are as follows: ⎡

2

0

0

(32)

40

1

Outputs

ˆ = 0.1 1 1 P1 = 10 I2 , P3 = 0.1, P4 = 10, w(0)     1 0 0 −516.3424 ˆ , A(0) = , P2 = 104 0 10 1 −63.2463   1 ˆ b(0) = , ε = 10−4 . 0.0486

T

30

3

Prediction error



2

20

10

2 1

Estimated

0

Following the frequency-domain step, a random input with uniform distribution on the interval [3.6, 4.4]sampled at 1 sample per second is applied and the input–output pair L {u(its ), y(its )}i=1 with ts = 1e − 4 and L = 38 is recorded. The frequency-domain  optimization is carried  out    n ˆ ≤ 4 , and over the space α, ˆ nˆ ∈ α, ˆ nˆ  0 <αˆ ≤ 1, 1 ≤  the optimal values are found as αopt , n opt = (0.7, 2). It is noteworthy that the calculation time for the proposed GA-based optimization is about 35 seconds, while the computational time for the equivalent grid search is 142 seconds. The computations are performed in MATLAB environment using a 2.50 GHz processor with 8 GB RAM. These results show that the GA drastically reduces the run time when compared to the deterministic grid search method. In the next step, the NFHS is excited by the sinusoidal input u(t) = 10 sin(10t), and the online identification algorithm is applied to estimateA, b, and w with the following initialization:

2

2

w

2

1



v(t) = e −0.2 (u(t)+5) + 2e −0.3 0.05s 0.7 +1 G(s) = VY (s) (s) = s 1.4 +60s 0.7 +500

w



Fig. 2 Validation results of Example 1

is measured in terms of the root-mean-square error (RMSE) and the relative error (RE) as follows: 1 RMSE = N

N  

yˆ − y i

i

2

!1/2 ,

i=1

 !1/2 i i 2 i=1 yˆ − y ,  N  i 2 i=1 y

N  RE =

(34)

where N , yˆ i , and y i are the number of data pair, the output of the model, and the system at the ith sample, respectively. The RMSE and RE are equal to 4.0401e − 005 and 0.0116, respectively. Example 2 (identification of the experimental heating system). In this example, the proposed NFHS identification

123

10

200

8

150

6 4 2

0

200

400

600

800

1000

1200

1400

200 150 100 50

0

200

400

600

800

1000

1200

1400

System output Model output

100 50 0 0

1600

Prediction error

output (degree C)

Outputs

Input (Volts)

M. Rahmani, M. Farrokhi

400

600

800

1000

1200

1400

1600

200

400

600

800

1000

1200

1400

1600

150 100 50 0

1600

200

0

Time (seconds)

Time (seconds)

Fig. 3 Input–output data set of Example 2

Fig. 4 validation results of Example 2 without NMS 200

Following the proposed frequency-domain identification, the aforementioned input–output data are transformed to the frequency domain with ts = 2 s and L = 801. The frequencydomain optimization is carried out, and the following linear fractional model as the best linear approximation of the heating system is estimated: ˆ G(s) =

105.5152s 1.2 − 56.8586s 0.6 + 16.6237 , (35) 16.0854s 1.8 + 9.3146s 1.2 + 4.5809s 0.6 + 1

ˆ As Fig. 4 shows, G(s) can mimic the dynamic behavior of the heating system with RMSE and RE equal to 49.6351 and 0.3590, respectively. In the time domain, the NMS is represented by

v(t) ˆ =

3 

  wˆ i (t) exp −0.52 (u(t) − 3 − 3(i − 1))2 , (36)

i=1

and the following initialization is used to estimate A, b and w:

123

Outputs

150 100 50 System output Model output

0 -50

Prediction error

algorithm is employed for the identification of a nonlinear heating system. The input is the voltage that drives a 300-Watt Halogen lamp suspended several inches above a thin metal plate and measured in volts. The output is the temperature in degrees Celsius measured by a thermocouple that is mounted on the back of the plate. This experiment is described in detail in Dullerud and Smith (1996), and its input–output data set is taken from Daisy: database for identification of systems (De moor 2015). The data set contains 801 samples of the input–output data that are measured with a sampling time of 2 seconds (Fig. 3).

0

200

400

600

800

1000

1200

1400

1600

200

400

600

800

1000

1200

1400

1600

10 5 0 -5 -10

Time (seconds)

Fig. 5 validation results of Example 2 with NMS

ˆ P1 = I3 , P3 = 0.01, P4 = 0.01, w(0) = 0, ⎡ ⎤ 0 0 −0.0622 ˆ = ⎣ 1 0 −0.2848 ⎦ , P2 = 10−4 I3 , A(0) 0 1 −0.5791 ⎡ ⎤ 1 ˆ b(0) = ⎣ 6.3473 ⎦ . −3.4203

(37)

Figure 5 shows the effectiveness of the proposed timedomain identification. The output of the identified neurofractional Hammerstein model follows the real system output with RMSE and RE equal to 6.1230 and 0.0442, respectively. To perform a fair comparison with respect to the widely used Hammerstein model, whose NMS is represented by the sum of known basis functions and its LDS is described by an infinite impulse response (IIR) filter, it is assumed that the NMS of the two models is the same. Moreover, the identification algorithm auxiliary model-based recursive least-squares

Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain… 200

integer-order Hammerstein models. The future work includes robustification of the proposed identification method and using a more sophisticated ANN.

Outputs

150 100 50

-50

Prediction error

System output Proposed method AM-RLS method

0 0

200

400

600

800

1000

1200

1400

1600

Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest.

Proposed method AM-RLS method

10 0

References

-10 0

500

1000

1500

Time (seconds)

Fig. 6 Comparison results between the proposed method and AMRLS method Table 1 RMSE and RE for identification methods of Example 2 Identification method

RMSE

RE

Proposed method

6.1230

0.0442

AM-RLS

8.3441

0.0603

algorithm (AM-RLS) described by Hu and Ding (2014) is utilized to estimate the coefficients of the IIR filter and weights of the NMS of the commonly used Hammerstein model. The identification quality of the proposed method in comparison with the AM-RLS method is illustrated in Fig. 6. As this figure shows, the proposed identification method performs much better than the AM-RLS method in terms of the RMSE and RE (Table 1). Moreover, only eight parameters in the proposed method must be updated, while 19 parameters must be updated using the IIR filter with order 8. Hence, it can be concluded that the proposed method can identify the heating system more accurately with less number of parameters.

5 Conclusion This paper investigated the neuro-fractional Hammerstein modeling and identification of nonlinear systems from input– output data. The proposed method comprised of two steps: frequency domain and time domain. In the frequency domain, the fractional order and the degree of the LDS were estimated. Then in the time domain and based on the Lyapunov method, other parameters of Hammerstein model were determined. In order to make the proposed method more practical, an observer was designed using an auxiliary model to estimate the system states. The simulating example showed that the proposed modeling method is superior to the sole

Aguila-Camacho N, Duarte-Mermoud MA, Gallegos JA (2014) Lyapunov functions for fractional order systems. Commun Nonlinear Sci Numer Simul 19:2951–2957. doi:10.1016/j.cnsns.2014.01. 022 Aoun M, Maltl R, Cois O (2002) Oustaloup A System identification using fractional Hammerstein models. In: Proceedings of the 15th IFAC world congress Benchellal A, Poinot T, Trigeassou JC (2006) Approximation and identification of diffusive interfaces by fractional models. Sig Process 86:2712–2727. doi:10.1016/j.sigpro.2006.02.025 Caponetto R, Graziani S, Sapuppo F, Tomasello V (2013) An enhanced fractional order model of ionic polymer-metal composites actuator. Adv Math Phys 2013:1–6. doi:10.1155/2013/717659 Chen H, Xiao Y, Ding F (2014) Hierarchical gradient parameter estimation algorithm for Hammerstein nonlinear systems using the key term separation principle. Appl Math Comput 247:1202–1210. doi:10.1016/j.amc.2014.09.070 Chen J (2014) Several gradient parameter estimation algorithms for dual-rate sampled systems. J Franklin Inst 351:543–554. doi:10. 1016/j.jfranklin.2013.08.016 Chen S, Hong X, Gao J, Harris CJ (2014b) Complex-valued B-spline neural networks for modeling and inverting Hammerstein systems. IEEE Trans Neural Netw Learn Syst 25:1673–1685 Chen Y-C, Wang J-S (2009) A Hammerstein–Wiener recurrent neural network with frequency-domain eigensystem realization algorithm for unknown system identification. J UCS 15:2547–2565 Cui M, Liu H, Li Z, Tang Y, Guan X (2014) Identification of Hammerstein model using functional link artificial neural network. Neurocomputing 142:419–428. doi:10.1016/j.neucom.2014.03.051 De moor BLR (2015) Daisy: database for the identification of systems, department of electrical engineering (ed) ESAT/SISTA, KU Leuven, Belgium. url: http://homes.Esat.Kuleuven.Be/smc/daisy Ding F, Liu X, Chen H, Yao G (2014) Hierarchical gradient based and hierarchical least squares based iterative parameter identification for CARARMA systems. Sig Process 97:31–39. doi:10.1016/j. sigpro.2013.10.018 Ding F, Liu X, Chu J (2013) Gradient-based and least-squares-based iterative algorithms for Hammerstein systems using the hierarchical identification principle. IET Control Theory Appl 7:176–184. doi:10.1049/iet-cta.2012.0313 Ding F, Liu X, Ma X (2016) Kalman state filtering based least squares iterative parameter estimation for observer canonical state space systems using decomposition. J Comput Appl Math 301:135–143 Djamah T, Mansouri R, Bettayeb M (2009) Djennoune S State space realization of fractional order systems. In: 2nd Mediterranean conference on intelligent systems and automation (CISA’09) vol 1. AIP Publishing, pp 37–42 Dullerud G, Smith R (1996) Sampled-data model validation: an algorithm and experimental application. Int J Robust Nonlinear Control 6:1065–1078

123

M. Rahmani, M. Farrokhi Fu Z-J, Xie W-F, Luo W-D (2013) Robust on-line nonlinear systems identification using multilayer dynamic neural networks with twotime scales. Neurocomputing 113:16–26 Gabano JD, Poinot T (2011a) Estimation of thermal parameters using fractional modelling. Sig Process 91:938–948. doi:10.1016/j. sigpro.2010.09.013 Gabano JD, Poinot T (2011b) Fractional modelling and identification of thermal systems. Sig Process 91:531–541. doi:10.1016/j.sigpro. 2010.02.005 Galvão RKH, Hadjiloucas S, Kienitz KH, Paiva HM, Afonso RJM (2013) Fractional order modeling of large three-dimensional RC networks. IEEE Trans Circuits Syst I Regul Pap 60:624–637 Garnier H, Young PC (2014) The advantages of directly identifying continuous-time transfer function models in practical applications. Int J Control 87:1319–1338 Guo Y, Wang Y, Sun G, Mao J (2012) Modeling of rate-dependent hysteresis for piezoelectric actuator with MPI model-based Hammerstein system. In: Intelligent robotics and applications. Springer, Berlin, pp 281-290 Hong X, Chen S, Harris CJ, Khalaf EF (2014) Single-carrier frequencydomain equalization for Hammerstein communication systems using complex-valued neural networks. IEEE Trans Signal Process 62:4467–4478. doi:10.1109/tsp.2014.2333555 Hu H, Ding R (2014) Least squares based iterative identification algorithms for input nonlinear controlled autoregressive systems based on the auxiliary model. Nonlinear Dyn 76:777–784 Hu P, Ding F, Sheng J (2013) Auxiliary model based least squares parameter estimation algorithm for feedback nonlinear systems using the hierarchical identification principle. J Franklin Inst 350:3248–3259. doi:10.1016/j.jfranklin.2013.07.013 Ibnkahla M (2012) Stochastic analysis of neural network modeling and identification of nonlinear memoryless MIMO systems. EURASIP J Adv Signal Process 2012:1–22 Jia L, Chiu M-S, Ge SS (2005) A noniterative neuro-fuzzy based identification method for Hammerstein processes. J Process Control 15:749–761. doi:10.1016/j.jprocont.2005.03.006 Li Y, Chen Y, Podlubny I (2009) Mittag-Leffler stability of fractional order nonlinear dynamic systems. Automatica 45:1965–1969 Liao Z, Zhu Z, Liang S, Peng C, Wang Y (2012) Subspace identification for fractional order Hammerstein systems based on instrumental variables. Int J Control Autom Syst 10:947–953. doi:10.1007/ s12555-012-0511-5

123

Lorenzo PR, Nalepa J, Kawulok M, Ramos LS, Pastor JR (2017) Particle swarm optimization for hyper-parameter selection in deep neural networks. In: Proceedings of the genetic and evolutionary computation conference. ACM, pp 481–488 Mkadem F, Boumaiza S (2009) Extended Hammerstein behavioral model using artificial neural networks. IEEE Trans Microw Theory Tech 57:745–751 Pintelon R, Schoukens J (2012) System identification: a frequencydomain approach. Wiley, New York Podlubny I (1999) Fractional differential equations : an introduction to fractional derivatives, to methods of their solution and some of their applications. Mathematics in Science and Engineering. Academic Press, New York Ren X, Lv X (2011) Identification of extended Hammerstein systems using dynamic self-optimizing neural networks. IEEE Trans Neural Netw 22:1169–1179 Van Pelt TH, Bernstein DS (2010) Non-linear system identification using Hammerstein and non-linear feedback models with piecewise linear static maps. Int J Control 74:1807–1823. doi:10.1080/ 00207170110089798 Wang D, Ding F, Ximei L (2013) Least squares algorithm for an input nonlinear system with a dynamic subspace state space model. Nonlinear Dyn 75:49–61. doi:10.1007/s11071-013-1048-8 Wang X, Ding F (2015) Recursive parameter and state estimation for an input nonlinear state space system using the hierarchical identification principle. Sig Process 117:208–218 Yu F, Mao Z, Jia M (2013) Recursive identification for Hammerstein– Wiener systems with dead-zone input nonlinearity. J Process Control 23:1108–1115. doi:10.1016/j.jprocont.2013.06.014 Zhao J, Ma X, Zhao S, Fei J (2014a) Hammerstein identification of supercharged boiler superheated steam pressure using LaguerreFuzzy model. Int J Heat Mass Transf 70:33–39. doi:10.1016/j. ijheatmasstransfer.2013.10.056 Zhao Y, Li Y, Chen Y (2014b)Complete parametric identification of fractional order Hammerstein systems. In: 2014 International conference on fractional differentiation and its applications (ICFDA). IEEE. pp 1–6