Fully parametric identification for continuous time fractional order Hammerstein systems

Fully parametric identification for continuous time fractional order Hammerstein systems

Journal Pre-proof

Fully parametric identification for continuous time fractional order Hammerstein systems Jiachang Wang, Yiheng Wei, Tianyu Liu, Ang Li, Yong Wang PII: DOI: Reference:

S0016-0032(19)30717-3 https://doi.org/10.1016/j.jfranklin.2019.10.001 FI 4198

To appear in:

Journal of the Franklin Institute

Received date: Revised date: Accepted date:

20 September 2018 9 June 2019 4 October 2019

Please cite this article as: Jiachang Wang, Yiheng Wei, Tianyu Liu, Ang Li, Yong Wang, Fully parametric identification for continuous time fractional order Hammerstein systems, Journal of the Franklin Institute (2019), doi: https://doi.org/10.1016/j.jfranklin.2019.10.001

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd on behalf of The Franklin Institute.

Fully parametric identification for continuous time fractional order Hammerstein systems Jiachang Wang, Yiheng Wei, Tianyu Liu, Ang Li, Yong Wang∗ Department of Automation, University of Science and Technology of China, Hefei 230026, China

Abstract This paper mainly investigates the issue of parameters identification for continuous time fractional order Hammerstein systems. When the commensurate order of linear part in system is prior information, two estimation algorithms are proposed to identify the parameters. Both methods use filters to obtain the derivative of the continuous time signal: one based on the Poisson filter and the other using the auxiliary model to construct the filter. Afterwards, the commensurate order update law is given, which contains the relationship between parameters and orders. In addition, due to the existence of fractional order terms, an irrational item will be introduced, and this paper gives a detailed mathematical derivation to discuss this issue. The effectiveness of proposed methods is demonstrated by the numerical simulations. Keywords: Hammerstein model, fractional order system, continuous time, system identification, least squares.

1. Introduction Mathematical models have always been an important topic in the areas of control, adaptive filtering, and system identification, etc. As one of the ✩ The work described in this paper was fully supported by the National Natural Science Foundation of China (No. 61601431, No. 61573332), the Anhui Provincial Natural Science Foundation (No. 1708085QF141) and the funds of China Scholarship Council (No. 201706340089, No. 201806345002). ∗ Corresponding author Email address: [email protected] (Yong Wang)

Preprint submitted to Journal of The Franklin Institute

October 14, 2019

mathematical models, the Hammerstein model shown in Fig. 1, consisting of a 5

static nonlinear block and a dynamic linear block, is able to accurately describe a wide variety of nonlinear objectives [1], such as fuel cells [2], chemical processes [3, 4], magneto-rheological dampers [5], and etc. Hence it has received extensive attention from engineers and researchers. 1RQOLQHDUEORFN

u t

1ā

/LQHDUEORFN

u t

/ā

H t y t

ͺн

y t

Figure 1: The block diagram of Hammerstein system.

As shown in Fig. 1, the dynamic linear block in conventional Hammer10

stein system is mostly an integer order system [6, 7, 8]. However, fractional order system (FOS) is more common than the traditional integer order system, especially for the application with long memory characteristics, such as supercapacitors [9], thermal diffusion processes [10, 11], and viscoelastic structures [12], etc. More and more researches have begun to focus on fractional order model

15

and found that it has a more accurate description for many physical processes [13, 14]. Therefore, the identification of Hammerstein system with fractional order model is an unavoidable problem. Currently, the research on identification is almost independent of the Hammerstein systems’ research. In terms of Hammerstein system identification,

20

some effective ideas have been proposed, such as the principle of key-term separation [15, 16] and the idea of multi-innovation [17, 18, 19], and used in such system identification. Meanwhile, there are many algorithms have achieved good results in this field, such as subspace method [20], gradient-based methods [21, 22] and so on[23, 24]. As far as FOS identification is concerned, there

25

are also many interesting and meaningful results [13, 14, 25]. Besides, in order to deal with the problem of fractional order derivation, modulation functions [26, 27], Poisson moment functions (PMF) [28, 29], and block pulse functions [30] are introduced into the parameters identification of FOSs. Furthermore,

2

there are also some innovation identification algorithms used for the parame30

ter estimation of FOS. For instances, the fractional order gradient method was proposed in [29, 31], and used for identification of FOSs contaminated by different noises. It should be pointed that the identification of fractional orders in the above researches is not involved or based on the optimization algorithm for order search.

35

Although many of the above results have been fruitful, there are still few identification results of Hammerstein system with fractional order models. In [32], the identification problem for discrete time fractional order Hammerstein system (FOHS) is discussed, but the continuous time model is a direct description for the real systems. Besides, the subspace identification method and a

40

hybrid frequency/time-domain approach are also proposed in [33, 34]. Though these methods are effective, the fractional order identification requires an order search by means of some swarm optimization algorithms, which obviously does not take into account the relationship between parameters and fractional orders. Motivated by the above discussions, this paper aims to design a fully para-

45

metric identification method for continuous time Hammerstein system with fractional order model. Firstly, the fractional order signal of FOHS is derived by PMF, and the parameters are estimated based on the least squares method. Then, the auxiliary model is introduced to identify FOHS, and the corresponding identification method is proposed. Afterwards, an order updating method

50

is given, which points out the relationship between system parameters and fractional order. Moreover, due to its particularity an irrational item will be introduced into the above relationship, which makes the corresponding irrational system unable to be solved directly, and this paper also gives a numerical simulation strategy for this problem.

55

The remainder of this paper is organized as follows. Section 2 provides some preliminaries of fractional calculus and problem description. Two parameter identification algorithms are presented in Section 3. Besides, the order identification method and corresponding simulation strategy are also discussed in this section. In Section 4, two numerical examples are provided to illustrate 3

60

the effectiveness of proposed method, and some conclusions are finally given in Section 5. 2. Preliminaries This section will give some essential fundamental knowledge used afterwards, and then the problem studied in this paper will be described.

65

2.1. Fractional calculus As a generalization of the classical calculus, there are three commonly used definitions for the fractional calculus: Gr¨ unwald–Letnikov definition, Riemann– Liouville definition and Caputo definition [35, 36]. The Riemann–Liouville fractional integral and Caputo fractional derivative are given as follows, which will

70

be used in this study. The α-order Riemann–Liouville fractional integral of a function f (t) ∈ R is defined as α t0 It f (t)

where Γ(α) = α t0 I t 75

R∞ 0

=

1 Γ(α)

Rt

t0

(t − τ )α−1 f (τ )dτ,

(1)

xα−1 e−x dx is the Gamma function. To simplify the notation,

is abbreviated as I α when t0 = 0.

With (1), the Caputo fractional derivative of f (t) can be expressed as α t0 Dt f (t)

m

d = t0 Itm−α dt m f (t) R t 1 = Γ(m−α) (t − τ )(m−α−1) f (m) (τ )dτ , t0

where m − 1 < α < m, m ∈ N+ , and similar to I α , D α is represented to

(2) α t0 Dt

when t0 = 0. 2.2. Problem formulation Reconsidering the Hammerstein system shown in Fig. 1, the linear block in 80

FOHS is a fractional order model. Thus, the system considered in this paper will be described by the following mathematical equation   A(sαi )y0 (t) = B(sαj )¯ u(t), 

y(t) = y0 (t) + ξ(t), 4

(3)

where i = 0, 1, . . . , Na , j = 0, 1, . . . , Nb , y(t) is the measured output signal, and ξ(t) is the additive Gaussian white noise. In (3), sα is defined as fractional order

85

differential operator, and then fractional order polynomials A(sαi ) = sαNa + PNa −1 PNb αi αj and B(sαj ) = i=0 ai s j=0 bj s , respectively. When it comes to the

commensurate order systems researched in this study, polynomials A(sαi ) and B(sαj ) are simplified as  P  A(sα ) = sNa α + Na −1 ai siα , i=0 PNb  α jα B(s ) = j=0 bj s ,

(4)

where the values of constants ai , bj and commensurate order α are unknown and will be identified.

90

Besides, the nonlinear block N(·) in Fig. 1 is a nonlinear function f (u), and it is assumed to be the sum of a series of prior known basis functions fk (·), k = 1, 2, . . . , Nc , which means u ¯(t) in (3) satisfies
PNc
u ¯(t) = cT f u(t) = k=1 ck fk u(t) ,

(5)

where ck are also the coefficients to be determined.

From the previous introduction, the main objective of this paper is to esti95

mate {a, b, c} = {ai , bj , ck } and α for a commensurate FOHS using the measured input and output signal {u(t), y(t)}. 3. Identification methods of FOHSs Firstly, this section will give two parameter identification methods when the fractional order is a prior information. Subsequently, a detailed derivation and

100

identification strategy will be given for the case of unknown fractional order. 3.1. Least squares identification based on PMF As aforementioned, the relationship between input u(t) and output y(t) of the entire system can be easily exported from (3), (4) and (5) sNa α y(t) +

PNa −1 i=0

ai siα y(t) =

PNb PNc j=0

5

jα k=1 bj ck s fk (u(t))

+ ν(t),

(6)

where ν(t) = sNa α ξ(t) + 105

sured noise.

PNa −1 i=0

ai siα ξ(t) is an additional item caused by mea-

In the identification of continuous time system, the significant difficulty lies in the incapacity to measure the fractional order differentiation of signals directly. Meanwhile, the fractional derivative obtained by definition or approximation methods will amplify the effect of noise. In this subsection, the PMF 110

approach for fractional derivatives is chosen to deal with this issue. The transfer function of PMF can be expressed as   β m Gm (s) = L gm (t) = s+λ ,

(7)

where β ∈ R+ , λ ∈ R+ , s is the Laplace variable and gm (t) is the m-order poisson moment function. Thus, (7) can be viewed as a cascade of m low pass filters with transfer function of 115

β s+λ .

The purpose of introducing PMF is to pre-filter the measurement signals. This not only suppresses the effect of noise, but also transfers the fractional derivative of measured signal y(t) and u(t) to gm (t). Taking the output y(t) as an example, one can acquire gm (t) ∗ D α y(t) = D α gm (t) ∗ y(t),

(8)

where ∗ stands for the convolution product. 120

In (8), the α-th order derivative of y(t) is transmitted to gm (t), which will avoid the issue of directly calculating fractional derivative of the measured signal. Since PMF is essentially equivalent to a low pass filter, its parameters should ensure that the bandwidth of the filter is not less than the one of the actual system. Some properties and more detailed knowledge about PMF can

125

be seen in [37]. By performing the pre-filter on (6), the following equation can be obtained ygNa α (t) +

PNa −1 i=0

ai ygiα (t) =

PNb PNc j=0

6

jα k=1 bj ck fkg (u(t))

+ νg (t),

(9)

where    ygiα (t) = D iα gm (t) ∗ y(t),    jα fkg (u(t)) = D jα gm (t) ∗ fk (u(t)),      ν (t) = D Na α g (t)∗ξ(t)+PNa −1 a D iα g (t)∗ξ(t). g m i m i=0

It is obvious that the estimated output can be easily written in a linear regression form from (9) ygNa α (t) = ΦT (t)θ + νg (t), 130

(10)

where

and

   T T T T  Φ = φT ∈ RNc (Nb +1)+Na , 1 , φ2 , . . . , φNc , φy  T  θ = bT ⊗ cT , aT ∈ RNc (Nb +1)+Na , 
T 0  φk=f Nb α u(t)
, f (Nb −1)α u(t)
, . . . , fkg u(t) ∈RNb +1, kg kg  φ =  − y (Na −1)α (t), −y (Na −2)α (t), . . . , −y 0 (t)T ∈RNa . y g g g

Then, if there are T sampling points, the target parameters will be obtained from least squares principle, and the corresponding estimator can be defined as θb =

135

h P T 1 T

t=1

Φ(t)ΦT (t)

i−1h P T 1 T

t=1

i Φ(t)ygNa α (t) .

(11)

b, b To acquire a b and cb directly, c1 is always set to 1, which is also discussed

b then a b, b in [22]. Denote θbq as the q-th element in θ, b and cb will be obtained b that is from θ,   T  b = θbNc (Nb +1)+1 , θbNc (Nb +1)+2 , . . . , θbNc (Nb +1)+Na , a      T b b = θb1 , θb2 , . . . , θbNb +1 ,   h i    cb = 1 PNb +1 1, θbNb +1+m , θb2Nb +2+m , . . . , θb(Nc −1)Nb +(Nc −1)+m T . m=1 b b b Nb +1 θm

θm

(12)

θm

It should be noted that the above estimator shown in (11) will give an unbiased estimate, which can be easily obtained from the following Lemma 1.

7

Lemma 1. [28] If a system whose impulse response is given as g(t) is excited 140

by Gaussian white noise ξ(t) with a mean 0 and variance σ 2 , the system will still output a signal with a mean of 0. That is E{g(t) ∗ ξ(t)} = 0,

(13)

where E{·} represents the mathematical expectation. With this lemma, it is easy to get hold of E{ν(t)} = 0, and ulteriorly, E{νg (t)} = 0. Meanwhile, one can also acquire the following equation from (11) n P  o −1  1 PT T Na α b E θ(t) =E T1 t=1 Φ(t)ΦT (t) (t) Φ(t)y g t=1 T n P   o P −1 T T T T =E t=1 Φ(t)Φ (t) t=1 Φ(t)Φ (t)θ n P −1  o T T Φ(t)νg (t) , +E t=1 Φ(t)Φ (t)

145

(14)

 b which means E θ(t) = θ. Thus, (11) is an unbiased estimator.

According to the previous derivation, the identification algorithm for FOHS

based on PMF can be summarized in Algorithm 1.

Algorithm 1: Least squares identification method for FHOS based on 150

PMF.

• Step 1: Collect the sampled data {u(t), y(t)}, and set parameters for PMF in (7). • Step 2: Pre-filter the input and output signal using PMF as shown in (8). 155

• Step 3: Construct the data matrices according to (9) and (10). • Step 4: Perform the least squares estimation according to (11) to obtain b the parameter estimation vector θ.

8

3.2. Identification for FOHS based on auxiliary model 160

Firstly, based on the prediction error method, an error relationship between the measured output and the estimated output can be easily represented as ε(t) = y(t) −

B(sα ) ¯(t) A(sα ) u

(15)

= A(sα )yf (t) − B(sα )¯ uf (t), where yf (t) = F (sα )y(t), u ¯f (t) = F (sα )¯ u(t), and F (sα ) =

1 A(sα )

is constructed

from the estimated parameters and equivalent to a filter essentially. Therefore, the identification purpose is equivalent to minimize ε(t). 165

However, the real pivotal problem lies in how we can get F (sα ), since both A(sα ) and B(sα ) are unknown, and it is the estimates of them that available. Therefore, an auxiliary model can be used in the identification process. The core idea of auxiliary model is to use the output, generated by the auxiliary model, to construct the information matrix, and the filter F (sα ) with the auxiliary

170

model’s parameters. Actually, it is an iterative process, and the auxiliary model is to generate the output as the instrumental variable in each iteration, and prefilter the measured signals. When the end of iteration is reached, the parameters of auxiliary model will approximate the real ones. During the iteration process, the transfer function of auxiliary model in the

175

b q (sα ) = q-th iteration is denoted as G

output is

bq (sα ) B bq (sα ) , A

then the corresponding noise-free

PNb jα bq (sα ) B j=0 bj,q s ˆ ˆ ybq (t) = u ¯(t) = u ¯(t), P N −1 a bq (sα ) A sNa α + i=0 ai,q siα

and the filter generated by auxiliary model is Fq (sα ) =

1 bq (sα ) . A

(16) According to

(15) and (16), the following filtered system equation can be expressed as Na α ybf,q (t)+

180

PNa −1 PNb PNc jα iα bf,q (t)= j=0 i=0 ai,q y k=1 bj,q ck,q fkf,q (u(t)),

(17)

jα iα in which, the variable ybf,q (t) = Fq (sα )[D iα yb(t)], and fkf,q (u(t)) = Fq (sα )[D jα fk (u(t))].

Thus, one can get the linear regression form of (17) Na α bT ybf,q (t) = Φ q (t)θq ,

9

(18)

  b q = φbT , φbT , . . . , φbT , φbT T∈RNc (Nb +1)+Na is the associated instruwhere Φ y,q 1,q 2,q Nc ,q

mental variable vector, and 
T 0  φbk,q=f Nb α u(t)
, f (Nb −1)α u(t)
, . . . , fkf,q u(t) , kf,q kf,q  φb

y,q =



(N −1)α

− ybf,qa

(N −2)α

(t), −b yf,qa

T 0 (t), . . . , −b yf,q (t) .

(19)

Taking T sample points into account, then, the (q + 1)-th estimated parameter value θbq+1 can be obtained by solving the following optimization problem

h P i h P i

T T Na 1 b q (t)ΦT b arg min T1 t=1 Φ q (t) θ− T t=1 Φq (t)yf,q (t) , θ

185

(20)

  T T T T where Φq = φT ∈ RNc (Nb +1)+Na . Note that φk,q = φbk,q 1,q , φ2,q , . . . , φNc ,q , φy,q

for the reason that fractional order α of linear part is fixed in this section.  T (N −1)α (N −2)α 0 Besides, φy,q = − yf,qa (t), −yf,qa (t), . . . , −yf,q (t) . Therefore, (20) results in the following solution

i−1 h P i hP T T Na α T b b θbq+1 = t=1 Φq (t)yf,q (t) . t=1 Φq (t)Φq (t)

(21)

The above equation (21) furnishes with the consistent estimate in each iteration 190

under the following conditions  P   T  b q (t)ΦT (t)  T1 t=1 E Φ is non singular, q T →∞     1 PT E Φ b q (t)ξ(t) = 0. t=1 T

(22)

T →∞

Though the issue handled here includes fractional order part and nonlinear

part, when the coupling between b and c is considered as a wholeness, the analysis will be similar to [38]. Therefore, the statistical property of the algorithm is consistent with the refined instrumental variable method discussed in [38]. 195

Now, the entire identification process can be summarized by Algorithm 2.

Algorithm 2: Identification for FHOS based on auxiliary model

• Step 1: Let q = 0, initialize the prefilter F0 (sα ), and set iteration number 200

Q, terminate parameter κ.

10

• Step 2: Construct the data matrices according to (18) and (20). • Step 3: Let q = q + 1, calculate θbq utilizing (21), and acquire Fq (sα ) by θbq .

• Step 4: Determine if q > Q or 205

to Step 2.

kθbq −θbq−1 k kθbq k

< κ is held. If neither met, return

3.3. Fractional order estimation As for the identification of commensurate order α, it might be not an arduous task when some swarm optimization algorithms are performed. However, 210

these methods ignore the relationship between the system parameters and the fractional order directly. In this section, an effective order identification method for the linear part of FOHS will be given, which will discuss the relationship between order update and parameter estimation. When both the order α and parameters {ai , bj , ck } need to be identified,

215

the corresponding algorithm can be divided into two identification stages: one for order update, the other for parameters estimation. Neither the previous Algorithm 1 nor Algorithm 2 involves fractional order update. One can define an error function as follows e(t) = y(t) − yb(t)

(23)

b αb )¯ = y(t) − G(s u(t).

Then the identification problem will be converted into minimizing the following 220

criterion function J=

1 2

PT

t=1

e2 (t).

(24)

With (23) and (24), the (m + 1)-th order α bm+1 can be easily obtained by the following update law

∂J α bm+1 = α bm − λH −1 ∂α , α=b α m

11

(25)

where λ is a positive real constant, and 
 ∂J = ∂e(t) T e(t) is the gradient, ∂α ∂α
T ∂e(t)  H = ∂e(t) ∂α ∂α is Hessian matrix.

(26)

The previous update law is obviously derived from the well known Gauss-

225

Newton algorithm, and other optimization algorithm, such as gradient descent method, may be equally effective. Thus, the complete identification procedure can be summarized as

Algorithm 3: Fully parametric identification algorithm for FHOS. 230

• Step 1: Let m = 0, set the maximum of order iteration M , and step size λ; initialize the order α bm = α0 ; set terminate parameter κ1 , κ2 .

• Step 2: Calculate the parameter vector θm using Algorithm 1 or 2, and compute Jm . 235

• Step 3: Calculate α bm+1 with (25) and (26), and perform Step 2 obtaining Jm+1 , and let λ = λ2 .

• Step 4: If Jm+1 > Jm and m + 1.

kθbm+1 −θbm k kθbm k

kb αm −b αm−1 k α bm

• Step 5: Determine if m > M or 240

> κ1 , return to Step 3; else m =

return to Step 2.

< κ2 is held. If neither met,

However, the real difficulty is not in the choice of the order update law when it comes to the fractional order estimation. Reconsidering the equation (26), one will find that 245

∂e(t) ∂α

plays a pivotal role in the fractional order update, and

the following equation can be acquired  ∂e(t) ∂ ∂α = ∂α y(t) =

h

 b α )¯ − G(s u(t)

b α ) ∂ A(s b α) B(s b2 (sα ) ∂α A

−

12

b α ) ∂ B(s b α) A(s b2 (sα ) ∂α A

i u ¯(t).

(27)

With (4), one can get  b α
P  ∂ A(s ) = Na sNa α + Na −1 iai siα ln(s), i=0 ∂α
PNb b α)  ∂ B(s jα = ln(s), j=0 jbj s ∂α

(28)

where ln(·) is natural logarithm function. Then, (27) can be denoted as b α ) ln(s)¯ u(t), = G(s

∂e(t) ∂α

b α) = where G(s

Na −1 b α )(Na sNa α + P iai siα ) B(s i=0

b2 (sα ) A

∂e(t) ∂α

The existence of ln(s) causes

250

−

Nb P

(29)

jbj sjα

j=0

b α) A(s

.

to become the output of an irrational sys-

tem, which is difficult to be handled with. A scheme of numerical approximation to this problem will be derived by the following lemma. Lemma 2. [39] For an s-domain function

ln(s) sα ,

where α > 0, its inverse

Laplace transform can be expressed as L −1 where ψ(α) = 255

d 1 Γ(α) dα Γ(α)

 ln(s) sα

=

tα−1 Γ(α) [ψ(α)

− ln(t)],

(30)

is digamma function.

In order to solve the above problem, the crucial task is to obtain output r(t) of a irrational system ln(s) with input u ¯(t), that is r(t) = ln(s)¯ u(t) = sα ln(s) ¯(t), sα u One can denote w(t) = w(t) = =

ln(s) ¯(t), sα u

 tα−1

and then

Γ(α) [ψ(α)

Rt

(31)

− ln(t)] ∗ u ¯(t)

(t−τ )α−1 Γ(α) [ψ(α) 0

(32)

− ln(t − τ )]¯ u(τ )dτ .

The following equation can be acquired from (1) and (32) w(t) = ψ(α)I α u ¯(t) −

1 Γ(α)

Rt 0

(t − τ )α−1 ln(t − τ )¯ u(τ )dτ.

(33)

Since the second part of the right side in (33) cannot be calculated directly, 260

it can be approximated by discretization, which means t = nh where n is ordinal 13

of the sampling point, and h is the sampling interval equaling to the sampling time Ts . Thus, the above equation can be approximated by w(nh) = ψ(α)I α u ¯(nh) −

Pn−1 α−1 ln(nh − ih)¯ u(ih)h. i=0 (nh − ih)

1 Γ(α)

(34)

Therefore, (27) can finally be calculated using the following formula ∂e(t) ∂α

b α )r(t) = G(s b α )[I α w(t)]. = G(s

(35)

To sum up, the entire identification process can be summarized as a block diagram shown in Fig. 2. f u

y t

G  sD

u t

н ͺ

T

D

3DUDPHWHUV ,GHQWLILFDWLRQ

yÖt

l  s D G

lf u

T  Dl

l D

et

J T Dl T Dl

2UGHU8SGDWH/DZ

Figure 2: The entire identification process of a FHOS. 265

Remark 1. In this paper, a comprehensive research on the identification of continuous time FOHSs is carried out. The main work can be summarized as follows: • The Poisson filter is used to obtain the derivative of the measured signal, 270

and then the estimated parameters are acquired by least squares method. • The concept of auxiliary model is adopted, and the corresponding estimation algorithm is proposed for the identification of FOHSs. • A commensurate order update method is given for the linear part of the system, which is completely different from the conventional method of glob-

275

ally searching by some swarm optimization algorithms. What is more, a numerical approximation method is given for the problem that it is difficult to solve irrational item output when the order updating. 14

4. Numerical examples and analysis In this section, the validity of the proposed method will be demonstrated 280

by two numerical examples. The commensurate order is assumed as prior information in the first example, and we will compare the performance of two identification algorithms. Then, the fully parametric identification will be performed using the order update law mentioned in previous subsection.

285

Consider a commensurate FOHS, whose linear part is described as follows   D 2α y0 (t) + a1 D α y0 (t) + a0 y0 (t) = b0 u ¯(t), (36)  y(t) = y0 (t) + ξ(t).

where ξ(t) is the Gaussian white noise, the commensurate order α = 1.5, and

a1 = 3, a0 = 2, b0 = 1.2 are the actual parameters of linear model. Considering that many nonlinearities can be approximated by Taylor expansion using polynomials, the nonlinear part of this example is described by a third order polynomial

290


u ¯(t) = f u(t) = c1 u(t) + c2 u2 (t) + c3 u3 (t),

(37)

where c1 = 1, c2 = 0.5, c3 = 0.25.

Example 1: Parameters estimation with the known commensurate order. In simulation, the simulation time is set to 100s, the sampling time Ts = 0.1s, and the input of system is a Gaussian random signal. The parameters of PMF are set to β = λ = 1.5, m = 4, and the parameters of auxiliary model are set to 295

1. Both Algorithm 1 and Algorithm 2 are performed to identify parameters, and the following relative error of parameters estimation is used to quantify the identification performance δ=

b kθ − θk × 100%. kθk

(38)

Monte Carlo experiment is carried out 200 times under different noise level, and 300

the corresponding identified results are as shown in Table. 1. In addition, to express the identification effect more intuitively, the simulation curves under SNR = 10dB are also recorded in Fig. 3 and Fig. 4. 15

Table 1: The results of 200 times Monte Carlo trials under different noise level

SNR

Algorithm

a0 = 2 ——————– mean std

a1 = 3 ——————– mean std

b0 = 1.2 ——————– mean std

c2 = 0.5 ——————– mean std

c3 = 0.25 ——————– mean std

1

1.9836 0.0231

2.9803 0.0296

1.1887 0.0331

0.5009 0.0132

0.2505 0.0128 0.73%

2

2.0003 0.0187

3.0000 0.0240

1.2018 0.0253

0.4996 0.0097

0.2495 0.0088 0.05%

1

1.9982 0.0069

2.9977 0.0082

1.1974 0.0124

0.5007 0.0051

0.2505 0.0042 0.10%

2

2.0004 0.0062

3.0005 0.0073

1.2008 0.0082

0.4998 0.0034

0.2499 0.0028 0.03%

δ

10dB

20dB

2

1.5

1

0.5

0

-0.5

-1

0

10

20

30

40

50

60

70

Figure 3: The comparison of output.

16

80

90

100

10-3

6

4

2

0

-2

-4

-6

0

10

20

30

40

50

60

70

80

90

100

Figure 4: Identified error with different algorithms.

It can be seen from Table. 1 that the two algorithms can accurately identify the system parameters at different noise levels, and the relative error δ reflects 305

the higher accuracy of Algorithm 2. As shown in Fig. 3 and Fig. 4, the output of identified models completely overcomes the effects of noise, whose magnitude of estimated error is close to 10−3 . At the same time, the error of identification model with Algorithm 1 is significantly smaller than that with Algorithm 2, which is consistent with the conclusion that the relative error of the parameter

310

in Table. 1 reflects. From the above results, the effectiveness of the two algorithms can be easily verified, but Algorithm 2 has more advantages. Obviously, the introduction of the auxiliary model in Algorithm 2 significantly improves the performance of parameter identification, which is confirmed both in the identification error

315

curve and in the statistical characteristics of estimated parameters. This result is consistent with the common sense, because the filter constructed based on the auxiliary model in Algorithm 2 has a bandwidth closer to the original system, which is helpful to retain more system information. Example 2: Fully parametric estimation for FOHS.

17

320

The system (36) is still researched, but the commensurate order α have no longer been known here. In this simulation, Algorithm 2 and the order update shown in (25) are used to complete the fully parametric identification of FOHS, and some initial parameters are set to: the fractional order α0 = 1, the step size λ = 10 and the terminate parameter κ1 = κ2 = 10−5 .

325

Remark 2. In general, the closer the initial value is selected to the actual order , the better the optimization, because it means fewer iterations. However, the problem is that we dont have any prior knowledge on actual order, so choosing a value close to the real order as the initial value is just an ideal idea. For simplicity, the initial value for α is set to α0 = 1 in our numerical simulation.

330

Since the fractional order is not a prior information, its update will introduce an irrational item ln(s). Firstly, the order update method and implementation method for irrational system mentioned in (25) and (34) will be verified under SNR = 20dB. Fig. 5 shows the order update curve in each iteration. It is obvious that the forementioned methods are very effective, because the order

335

update eventually converges, and at the end of the iteration, the final updated order α bend = 1.4941, which is very close to the real order α = 1.5.

Besides, Monte Carlo experiment with SNR = 20dB is also performed 200

times, and the identified parameters are counted in Fig. 6. As can be seen from Fig. 6, all of the system parameters are recovered despite some deviations. The

340

reason for existence of bias is not difficult to understand, because the fractional order system in the simulation is implemented by frequency domain approximation, while the identification algorithm is performed in the time domain. Thus, the approximation error is bound to affect the identification accuracy. Fortunately, the relative error is only δ = 0.58%, which is quite acceptable.

345

5. Conclusions In this paper, the identification issue of Hammerstein systems with fractional order model has been discussed, and two identification methods are proposed for

18

1.6

1.5

1.4

1.5

1.3

1.49

1.2

30

31

32

1.1

1

5

10

15

20

25

Figure 5: The estimated fractional order in each iteration.

2.10

1.52 3.10

1.50

1.48

3.00

2.00

2.90 1.90

1.47

1.26

1.20

2.80

0.27

0.52

0.25

0.50

0.24 1.14

0.48

0.23

Figure 6: The results of fully parametric identification.

19

30

the case where the commensurate order is given. Furthermore, the corresponding order identification method is given for cases where the order is unknown. 350

Under this scheme, a numerical simulation strategy is proposed to handle with the irrational item when estimating the commensurate order. The simulation results show that the proposed fully parametric identification algorithm and simulation strategy are completely effective. In the further work, the identification issues of system with parameter correlation and the improvement of

355

estimation accuracy are worth researching.

References References [1] F. Giri, E. W. Bai, Block-oriented nonlinear system identification, Springer, London, 2010. 360

[2] F. Jurado, A method for the identification of solid oxide fuel cells using a Hammerstein model, Journal of Power Sources 154 (1) (2006) 145–152. [3] J. G. Smith, S. Kamat, K. Madhavan, Modeling of pH process using wavenet based Hammerstein model, Journal of Process Control 17 (6) (2007) 551–561.

365

[4] Z. Y. Zou, G. P. Liu, N. Guo, Predictive control of nonlinear Hammerstein systems and application to pH processes, in: 2003 European Control Conference (ECC), IEEE, Cambridge, UK, 2003, pp. 3064–3069. [5] J. D. Wang, A. Sano, T. W. Chen, B. Huang, Identification of Hammerstein systems without explicit parameterisation of non-linearity, International

370

Journal of Control 82 (5) (2009) 937–952. [6] B. Q. Mu, H. F. Chen, L. Y. Wang, G. Yin, W. X. Zheng, Recursive identification of Hammerstein systems: convergence rate and asymptotic normality, IEEE Transactions on Automatic Control 62 (7) (2017) 3277– 3292. 20

375

[7] F. Ding, X. P. P. Liu, G. J. Liu, Identification methods for Hammerstein nonlinear systems, Digital Signal Processing 21 (2) (2011) 215–238. [8] Y. Liu, E. W. Bai, Iterative identification of Hammerstein systems, Automatica 43 (2) (2007) 346–354. [9] C. F. Zou, L. Zhang, X. S. Hu, Z. P. Wang, T. Wik, M. Pecht, A review

380

of fractional-order techniques applied to lithium-ion batteries, lead-acid batteries, and supercapacitors, Journal of Power Sources 390 (2018) 286– 296. [10] S. Z. Chen, F. W. Liu, I. Turner, X. L. Hu, Numerical inversion of the fractional derivative index and surface thermal flux for an anomalous heat

385

conduction model in a multi-layer medium, Applied Mathematical Modelling 59 (2018) 514–526. [11] S. Victor, R. Malti, H. Garnier, A. Oustaloup, Parameter and differentiation order estimation in fractional models, Automatica 49 (4) (2013) 926–935.

390

[12] J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute 354 (17) (2017) 7917–7945. [13] R. Shalaby, E. H. Mohammad, A. Z. Belal, Fractional order modeling and control for under-actuated inverted pendulum, Communications in Non-

395

linear Science and Numerical Simulationdoi:10.1016/j.cnsns.2019.02. 023. [14] T. Abuaisha, J. Kertzscher, Fractional-order modelling and parameter identification of electrical coils, Fractional Calculus and Applied Analysis 22 (1) (2019) 193–216.

400

[15] H. B. Chen, Y. S. Xiao, F. Ding, Hierarchical gradient parameter estimation algorithm for Hammerstein nonlinear systems using the key term separation principle, Applied Mathematics and Computation 247 (2014) 1202–1210. 21

[16] F. Ding, H. B. Chen, L. Xu, J. Y. Dai, Q. S. Li, T. Hayat, A hierarchical least squares identification algorithm for Hammerstein nonlinear systems 405

using the key term separation, Journal of the Franklin Institute 355 (8) (2018) 3737–3752. [17] F. Ding, T. W. Chen, Performance analysis of multi-innovation gradient type identification methods, Automatica 43 (1) (2007) 1–14. [18] F. Ding, Several multi-innovation identification methods, Digital Signal

410

Processing 20 (4) (2010) 1027–1039. [19] P. Ma, F. Ding, New gradient based identification methods for multivariate pseudo-linear systems using the multi-innovation and the data filtering, Journal of the Franklin Institute 354 (3) (2017) 1568–1583. [20] I. Goethals, K. Pelckmans, J. A. Suykens, B. De Moor, Subspace identifica-

415

tion of Hammerstein systems using least squares support vector machines, IEEE Transactions on Automatic Control 50 (10) (2005) 1509–1519. [21] F. Ding, Y. Shi, T. W. Chen, Gradient-based identification methods for Hammerstein nonlinear ARMAX models, Nonlinear Dynamics 45 (1-2) (2006) 31–43.

420

[22] S. S. Cheng, Y. H. Wei, D. Sheng, Y. Q. Chen, Y. Wang, Identification for Hammerstein nonlinear ARMAX systems based on multi-innovation fractional order stochastic gradient, Signal Processing 142 (2018) 1–10. [23] L. Vanbeylen, R. Pintelon, J. Schoukens, Blind maximum likelihood identification of Hammerstein systems, Automatica 44 (12) (2008) 3139–3146.

425

[24] S. Boubaker, Identification of nonlinear Hammerstein system using mixed integer-real coded particle swarm optimization: application to the electric daily peak-load forecasting, Nonlinear Dynamics 90 (2) (2017) 797–814. [25] S. M. Fahim, S. Ahmed, S. A. Imtiaz, Fractional order model identification using the sinusoidal input, ISA Transactions 83 (2018) 35–41. 22

430

[26] D. Y. Liu, T.-M. Laleg-Kirati, O. Gibaru, W. Perruquetti, Identification of fractional order systems using modulating functions method, in: American Control Conference, IEEE, Washington, USA, 2013, pp. 1679–1684. [27] Y. Dai, Y. H. Wei, Y. S. Hu, Y. Wang, Modulating function-based identification for fractional order systems, Neurocomputing 173 (2016) 1959–1966.

435

[28] Y. S. Hu, Y. Fan, Y. H. Wei, Y. Wang, Q. Liang, Subspace-based continuous-time identification of fractional order systems from nonuniformly sampled data, International Journal of Systems Science 47 (1) (2016) 122–134. [29] R. Z. Cui, Y. H. Wei, Y. Q. Chen, S. S. Cheng, Y. Wang, An innovative pa-

440

rameter estimation for fractional-order systems in the presence of outliers, Nonlinear Dynamics 89 (1) (2017) 453–463. [30] Y. G. Tang, H. F. Liu, W. W. Wang, Q. S. Lian, X. P. Guan, Parameter identification of fractional order systems using block pulse functions, Signal Processing 107 (2015) 272–281.

445

[31] R. Z. Cui, Y. H. Wei, S. S. Cheng, Y. Wang, An innovative parameter estimation for fractional order systems with impulse noise, ISA Transactions 82 (2018) 120–129. [32] D. Ivanov, Identification discrete fractional order Hammerstein systems, in: International Siberian Conference on Control and Communications, IEEE,

450

Omsk, Russia, 2015, pp. 1–4. [33] Z. Liao, Z. T. Zhu, S. Liang, C. Peng, Y. Wang, Subspace identification for fractional order Hammerstein systems based on instrumental variables, International Journal of Control, Automation and Systems 10 (5) (2012) 947–953.

455

[34] M.-R. Rahmani, M. Farrokhi, Identification of neuro-fractional Hammerstein systems: a hybrid frequency-/time-domain approach, Soft Computing (2017) 1–10. 23

[35] Y. H. Wei, Y. Q. Chen, S. S. Cheng, Y. Wang, A note on short memory principle of fractional calculus, Fractional Calculus and Applied Analysis 460

20 (6) (2017) 1382–1404. [36] I. Podlubny, Fractional Differential Equations: an Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, Academic Press, San Diego, 1999. [37] T. Bastogne, H. Garnier, P. Sibille, A PMF-based subspace method for

465

continuous-time model identification. application to a multivariable winding process, International Journal of Control 74 (2) (2001) 118–132. [38] P. C. Young, H. Garnier, M. Gilson, Identification of Continuous-time Models from Sampled Data, Springer, London, 2008. [39] A. D. Polyanin, A. V. Manzhirov, Handbook of Integral Equations: Second

470

Edition, CRC Press, New York, 2008.

24