# Sequential Parameter Identification of Fractional-Order Duffing System Based on Differential Evolution Algorithm

Authors: Li Lai; Yuan-Dong Ji; Su-Chuan Zhong; Lu Zhang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Using the dynamic properties of fractional-order Duffing system, a sequential parameter identification method based on differential evolution optimization algorithm is proposed for the fractional-order Duffing system. Due to the step by step parameter identification method, the dimension of parameter identification of each step is greatly reduced and the search capability of the differential evolution algorithm has been greatly improved. The simulation results show that the proposed method has higher convergence reliability and accuracy of identification and also has high robustness in the presence of measurement noise.


---

## Body

## 1. Introduction

Since the gradual development of chaos control theory in the 1990s, nonlinear system identification has gradually developed into an important direction of modern control theory [1]. Compared with traditional integer-order calculus, fractional-order calculus can better depict various materials and processes with memory; hence, theory of fractional-order calculus has become a popular topic in the development of nonlinear science [2]. With the continuous development of fractional-order calculus theory, fractional-order nonlinear system identification has gradually attracted the attention of scholars [3–5]. The general system identification problem can be converted into an optimization problem, which can be solved; traditional identification methods, such as step response, impulse response, frequency response, least squares, and maximum likelihood, exhibit good recognition performance for linear system identification problems; however, they cannot frequently obtain satisfactory identification results for nonlinear system identification [6].In recent years, intelligent optimization algorithms, such as the genetic algorithm, the ant colony algorithm, and particle swarm optimization (PSO), have been applied to system identification, thereby opening up a new train of thought for nonlinear system identification problems [7–9]. However, the introduction of fractional orders and the complexity of nonlinear systems have made fractional-order nonlinear system identification more difficult than integer-order system identification; consequently, achieving accurate identification is more difficult when using traditional intelligent optimization algorithms for nonlinear system parameters.PSO is a swarm intelligence optimization algorithm based on the foraging behavior of birds; it has easy-to-understand features, less parameters, and simple programming, and easy implementation [10]. However, PSO also has defects, such as low optimization efficiency and the tendency to be easily trapped in the local optimum when the problem dimension is too high [11]. Yuan and Yang combined PSO and active control theory to realize fractional-order chaotic system parameter identification and synchronous control [12]. Huang et al. used PSO based on quantum parallel characteristics for the parameter identification of fractional-order Lorenz system and Chen system [13].The differential evolution (DE) algorithm, proposed by Storn et al. in 1997, is a new method based on “the survival of the fittest” rule of the intelligent evolutionary algorithm; it has received considerable attention and has extensive applications because of its simple operation, less control parameters, high stability, and powerful global optimization capability [14, 15]. Compared with other evolution algorithms, DE applies a global search strategy based on population. It uses real number coding based on the simple mutation operation of difference and one-on-one competition survival strategy to reduce the complexity of genetic operations. Simultaneously, the unique memory capability of DE enables it to dynamically track current search conditions to adjust its search strategy, thereby achieving powerful global convergence capability and robustness. In 2012, Tang et al. realized parameter identification on a fractional-order Liu system in the same order using DE [16]. Zhu et al. completed parameter identification on a fractional-order chaotic system with measurement noise based on an improved DE algorithm [17].As a matter of fact, most of the existing optimization algorithms based parameter identification methods for fractional-order chaotic systems use only one objective function. The major challenge is that an intelligent optimization algorithm may be trapped in the local optima of the objective function, which leads to a slow convergence speed and low identification accuracy. This issue is particularly challenging when the dimension of the problem is high and there are numerous local optima [18]. Although the DE algorithm has the advantages of simple structure, ease of use, and high speed of convergence, the DE algorithm may also easily be trapped in the local optimum because of the high identification dimension [19].The Duffing system is a mathematical model that is widely used in mechanical vibration, physics, biology, neuroscience, and other fields; its rich nonlinear dynamic behavior has attracted the interest of numerous scholars. Aguilar-Ibáñez et al. achieved parameter identification on an integer-order Duffing system using an improved gradient algorithm [20]. The traditional integer-order damping model can be described accurately using fractional-order calculus. Scholars have recently studied stability, bifurcation, chaos, and other dynamic characteristics of the fractional-order Duffing equation [21–24]; however, research on the parameter identification problem of a fractional-order Duffing system has not yet been conducted.For the fractional-order Duffing system, [25] studied its steady attractor and found out that the state transition from “vibration stopped (VS) state” to “small-scale periodic (SSP) state” of steady attractor is sensitive to regular signals but immune to noise. In the adjacent area of the steady attractor, a small perturbation of the parameters of the periodic driving force can lead to a fundamental change of the system dynamic behavior. Based on these properties of the steady attractor, the authors of [25] proposed a detecting method for weak periodic sinusoidal signal. Moreover, our study has found that a reasonable application of the above properties of the steady attractor can also provide an efficient parameter estimation method for amplitude and frequency estimation of a sinusoidal signal. Note that the amplitude and frequency of the sinusoidal signal are also parameters to be identified of the fractional-order Duffing system.Therefore, inspired by the properties of the steady attractor of the fractional-order Duffing system and the relationship between the two output state, this study firstly finds out that the original 6-dimensional identification problem of the fractional-order Duffing system can be transformed to a successive step by step identification problem, which only has low identification dimensions. And then a sequential parameter identification method is established based on the DE optimization algorithm, which we called sequential differential evolution (SDE) method. The simulation results show that when compared with traditional DE based optimization algorithm (we called DE method), the proposed SDE method in this study can obtain better parameter identification results with and without noise. Simultaneously, the comparison of the simulation results of the sequential parameter identification method based on PSO algorithm (we called sequential particle swarm optimization (SPSO) method) also shows that the DE algorithm exhibits performance superior to that of PSO in nonlinear parameter identification problems.


## 2. Fractional-Order Duffing System and Its Numerical Algorithm

In this paper, we consider the following state equations given by the fractional-order Duffing system [23]:(1)Dp0x1=x2,Dp0x2=a0x1-b0x13-k0x2+r0cos⁡ω0t,where Dp0x is the fractional-order derivative based on the definition of Caputo [2]:(2)Dp0x=1Γn-p0∫0tfnτt-τ1+p0-ndτ,n-1<p0<ndnxdtn,p0=nin which, Γ(·) is the gamma function and Γ(z)=∫0+∞tz-1e-1dt. When p0=1, (1) is the classical Duffing equation.The observation equation of (1) with additive noise is(3)Y=X+W,where X=[x1(t),x2(t)]T is state vector, Y=[y1(t),y2(t)]T is the output observation vector of system, and W=[w1(t),w2(t)]T is the Gaussian white noise vector, in which, w1(t) and w2(t) are not related to each other.Here,a0, b0, k0, amplitude r0, angular frequency ω0, and order p0 are unknown parameters to be identified. In this paper, we suppose these parameters are all positive real numbers and p0∈(0,1].The numerical algorithm for solving the fractional-order Duffing equation (1) is based on the numerical algorithm given in [23]:(4)x1tk=x2tk-1hp0-∑j=vkcjp0x1tk-j,x2tk=a0x1tk-b0x1tk3-k0x2tk-1+r0cosω0tkhp0-∑j=vkcjp0x2tk-j,where (5)c0p0=1,cjp0=1-1+p0jcj-1p0,j=1,2,….


## 3. Differential Evolution Algorithm

Consider a general optimization problem:(6)Θ^=argΘminfΘ,where f is the objective function, Θ=[θ1,θ2,…,θJ] is the J-dimensional parameter vector, and θj∈[Minθj,Maxθj], j=1,2,…,J.The DE algorithm is a kind of evolutionary algorithms, which has the advantages of simple structure, ease of use, and high speed of convergence [14]. The main steps of DE algorithm to search the extreme value of the objective function are given as follows.(1) Initialization Operation. The DE algorithm randomly generates an initial population with size of NP, and then the initial value of the ith individual is(7)Θi,0=θ1i,0,θ2i,0,…,θJi,0,i=1,2,…,NP,in which, θj(i,0)=Minθj+[Maxθj-Minθj]×rand(1), j=1,2,…,J, and rand(1) denotes random numbers generated by a uniform distribution. At the same time, set up the mutation factor F∈[0,2], crossover probability CR, and maximum evolution times L.(2) Mutation Operation. The DE algorithm uses two difference individual vectors randomly chosen from the population as the random change in the source of third individual; the weighted difference vector of the two individual vectors is added to the third individual according to certain rules. For each individual target Θ(i,l), the mutation vector is produced by the following difference strategy:(8)Vi,l+1=v1i,l+1,v2i,l+1,…,vJi,l+1=Θr1,l+F×Θr2,l-Θr3,l,where the random selected numbers r1, r2, and r3 are not the same, r1, r2, r3, and the sequence number i of the target vector are also not the same. In the process of evolution, in order to ensure the effectiveness of the solution, if the mutation vector does not meet the boundary conditions, a new mutation vector is regenerated using the same random method in the step of initialization.(3) Crossover Operation. In order to increase the diversity of the parameter vector, for the lth generation and its mutation population, lth generation test vector was introduced by the crossover operation (9)Ui,l+1=u1i,l+1,u2i,l+1,…,uJi,l+1,in which,(10)uji,l+1=vji,l+1,if  rand0,1≤CRθji,l,otherwise,j=1,2,…,J.(4) Selection Operation. The DE algorithm selects the next generation of the population using the greedy algorithm; if the objective function values of the test individuals are less than those of the target individual, then the test individuals replace the target individual in the next generation; otherwise, target individuals still preserved:(11)Θi,l+1=Ui,l+1,iffUi,l+1≤fΘi,lΘi,l,otherwise,i=1,2,…,NP.


## 4. Parameter Identification of Fractional-Order Duffing System

Based on the dynamic characteristics of the fractional-order Duffing system, the parameter identification method proposed in this paper will be divided into three steps:Step 1.
Identification of orderp0;Step 2.
Identification of amplituder0 and angular frequency ω0;Step 3.
Identification ofa0, b0, and k0.Based on the above three steps, the original 6-dimensional identification problem is transformed into a successive step by step identification problem, including a 1-dimension identification problem, a 2-dimension identification problem, and a 3-dimensional identification problem, thereby greatly reducing the search dimension and improving the search efficiency.


### 4.1. Identification of Orderp0

From the first equation of the fractional-order Duffing system (1) and the measurement equation (3) we can see, in the absence of noise, the second component y2(t) of the final output vector Y(t) is the fractional derivative of the first component y1(t) with order p0, that is, Dp0y1=y2.SupposeDpy1=yp, then if and only if p=p0, the energy difference of yp(t) and y2(t) reaches its minimum value. The estimated value of order p0 can be obtained by adjusting the order of p and searching for the minimum value of the energy difference:(12)p^0=argminp∈0,1yp-y22,where x(t)2 denotes the energy of x(t).Therefore, in the first step of our identification method, the objective functionf1 of DE algorithm is chosen as the energy difference of the fractional-order derivative Dpy1 and y2(t): (13)fSDE1p=Dpy1-y22;the order p0 is then estimated by using the DE algorithm.


### 4.2. Identification of Amplituder0 and Angular Frequency ω0

#### 4.2.1. Oscillation Variance of the Output of the Fractional-Order Duffing System

For eachp0∈(0,1], there are three equilibrium points of the fractional-order Duffing system (1), (0,0), and (±a0/b0,0).For the equilibrium point(±a0/b0,0), its Jacobian matrix is given by(14)H=01-2a0-k0.Then there are two cases:(1) Ifk02≥8a0, the characteristic values of H are two negative numbers:(15)λ1,2=-k0±k02-8a02,whose arguments are both π.(2) Ifk02<8a0, the characteristic values of H are a pair of conjugate complex numbers:(16)λ1,2=-k0±i8a0-k022,whose arguments are(17)arg⁡λ1=π2+atan⁡k08a0-k02,arg⁡λ2=π+atan⁡8a0-k02k0.Since0<p≤1, then both the two cases satisfy argλi>p/2π, i=1,2. From the conclusions in [26], the equilibrium points (±a0/b0,0) are stable focus points; that is, the phase trajectory from any initial value will eventually be attracted and converge to these two points.For the equilibrium point(0,0), its Jacobian matrix is given by(18)H=01a0-k0;the characteristic values of H are one negative and one positive number λ1,2=-k0±k02+4a0/2. Thus, (0,0) is a saddle point; it is neither an attractor nor a repeller.Similar to the integer-order Duffing system, the outputX(t) of the fractional-order Duffing equation (1) varies as the amplitude r0 of the periodic driving force varies. In the following section, we set the parameters as a0=1, b0=1, k0=1, and ω0=1, and the fractional-order p0 is 0.96. We set the amplitude r0 as 0, 0.1, 0.7, 0.922, 0.933, and 1.8, respectively, and then compare and analyze the first component x1(t) of the output X(t) of the system.(1) Whenr0=0, for any given initial value, the output response x1(t) oscillates around one of the stable focus points and the oscillation amplitude gradually reduces to zero; the system will converge to the phase trajectory of the focus, which reached a steady state. The corresponding phase trajectory and output time response are shown in Figures 1(a) and 1(b), respectively.Figure 1
The phase trajectory and output time response of the fractional-order system withr0=0 and r0=0.1.
(a)
(b)
(c)
(d)Whenr0>0, since the Duffing system is a bistable system, which has two stable focus points and a saddle point, in the presence of the periodic driven force r0cos⁡(ω0t), there is a critical threshold Ac in the bistable system: if 0<r0<Ac, then the output response x1(t) will produce local periodic motion near one steady attractor when the time is long enough. Then x1(t) has the same frequency of the periodic driving force. At this time, the phase plane trajectory shows the attractor under Poincaré mapping. The corresponding phase trajectory and output time response are shown in Figures 1(c) and 1(d), respectively.(2) Whenr0>Ac, x1(t) will create a wide-range transition around the two steady attractors; as the value of r0 increases, the oscillation amplitude of x1(t) also increases, and the phase trajectory will go through four states in turns, namely, “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state,” as shown in Figure 2.Figure 2
The phase trajectory and output time response of the fractional-order system of “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state.”
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It can be observed that the first componentx1(t) of system output around the equilibrium point has increasing oscillation amplitude as the amplitude r0 increases. Note that the horizontal coordinate of the focus point can be regarded as the mean value of x1(t) at time domain, and the variance of x1(t) at time domain exactly characterizes the oscillation degree of the signal around its mean value.Suppose the first component of the output of (1) with initial value X0=(x0,y0) is x1(t); we may assume that its steady state response starts at time T0; then theoscillation variance of x1(t) on [T0,T] is defined by(19)σx=1T-T0∫T0Tx1t-mx2dt,in which, mx=1/T-T0∫T0Tx1(t)dt denotes the mean value of x1(t) in [T0,T].With a fixed angular frequencyω0, the oscillation variance σx is a function of the input amplitude r0. The smaller the oscillation variance of output time response is, the closer the phase state is to the steady state, and the closer the corresponding input signal is to a zero value signal and vice versa. In this paper, we use theoscillation varianceσx of x1(t) to estimate the amplitude of the input periodic driving force.Figure3 shows the oscillation variance σx of fractional-order Duffing system with different orders varies as a function of r0, in which, sampling interval h=0.01 s and sampling time T=20π s.Figure 3
Curves of the oscillation variance of the output of the fractional-order Duffing system as functions of the amplituder0 of the input periodic driving force.From Figure3 we can see, for any fractional-order Duffing system, the oscillation variance σx is an increasing function of absolute value of the amplitude r0 and reaches its minimum value at r0=0. In Section 4.2.2, this dynamic characteristic of fractional-order Duffing system is used to estimate the unknown parameters of its periodic driving force.


#### 4.2.2. Identification ofr0 and ω0 Based on the Oscillation Variance

Inputu=-rcos(ωt) into the fractional-order Duffing system (1); we obtained the following controlled system:(20)Dp0x1=x2,Dp0x2=a0x1-b0x13-k0x2+Rcos⁡ω0t+ϕ.Supposer0+Δr=r,ω0+Δω=ω, then(21)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2ϕ=atan⁡r0+Δrsin⁡Δωtr0-r0+Δrcos⁡Δωt.Since1≥cos(Δωt)≥-1 then we have(22)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2≥Δr≥0.Thus, R=0 if and only if (23)Δr=0,Δω=0⟷r=r0,ω=ω0.For anyr, ω, the oscillation variance of the first output component x1(t) of (20) is a function of r, ω:(24)σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt.Then from the above discussion in Section 4.2.1, we know that R=0 if and only if r=r0, ω=ω0, while R=0 if and only if the output is in the steady state, and in this case, σ(r,ω) reaches at its minimum. Thus, the amplitude r0 and angular frequency ω0 can be estimated by searching the minimum value of the output oscillation variance of x1(t):(25)r^0,ω^0=argminr,ωσr,ω,r∈0,rmax,ω∈0,ωmax,in which, rmax and ωmax are the search upper bounds of r and ω, respectively.Therefore, in the second step of our identification method, the objective function of DE algorithm is chosen as the oscillation variance of (20): (26)fSDE2r,ω=σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt;the amplitude r0 and angular frequency ω0 are then estimated by using the DE algorithm.


### 4.3. Identification ofa0, b0, and k0

After identification of the orderp0, the amplitude r0, and angular frequency ω0, we obtained the estimated values p^0, r^0, and ω^0, and then the remaining parameters to be identified are a0, b0, and k0. The identification of these three parameters can still be considered as a multidimensional continuous optimization problem. Thus, the approximate response system corresponding to system (3) with noise is(27)Dp^0z1=z2,Dp^0z2=az1-bz13-kz2+r^0cos⁡ω^0t.Therefore, in the last step of our identification method, the objective function of DE algorithm is chosen as the energy difference between the first component of the outputy1(t) in driving system (3) and the first component z1(t) of the output in approximate response system (27):(28)fSDE3a,b,k=y1-z12.The remaining parameters a0, b0, and k0 are then estimated by using the DE algorithm.


### 4.4. Design Steps of the Proposed SDE Method

The whole design steps for parameter identification of fractional-order Duffing system using SDE method can be summarized as follows.Step 1.
Identification of orderp0.
Step 1.1. Choose control parameters of DE, that is, population size NP, mutation factor F, crossover probability CR, and the maximum evolution times L. Randomly initialize NP individuals in the search space as Θ(i,l)=p(i,l), i=1,2,…,NP. Let iteration number l=0.
Step 1.2. Solve the systems (1) and (3) under initial state X0 and obtain output states Y=[y1,y2]T using the numerical method given by (4)-(5). For the ith individual, obtain output states of Dp(i,l)y1=yp from initial state X0, and then calculate its objective function according to (13).
Step 1.3. Perform mutation operation according to (8). Perform crossover operation according to (9) to obtain crossover trial vectors. Calculate the objective function of crossover trial vectors by (13). Perform selection operation according to (11) to generate individuals for next generation.
Step 1.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution p^0; otherwise, let l=l+1 and go to Step 1.3.Step 2.
Identification of amplituder0 and angular frequency ω0.
Step 2.1. Randomly initialize NP individuals as Θ(i,l)=[r(i,l),ω(i,l)], i=1,2,…,NP. Let iteration number l=0.
Step 2.2. For the ith individual, input u=-r(i,l)cos(ω(i,l)t) into the fractional-order Duffing system (1) and get the estimated system (20), obtain output states of (20) from initial state X0, and then calculate its objective function according to (26).
Step 2.3. The same as Step  1.3.
Step 2.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution r^0, ω^0; otherwise, let l=l+1 and go to Step 2.3.Step 3.
Identification of parametersa0, b0, and k0.
Step 3.1. Randomly initialize NP individuals Θ(i,k)=[a(i,l),b(i,l),k(i,l)], i=1,2,…,NP. Let iteration number l=0.
Step 3.2. For the ith individual, solve the approximate response system (27) with estimated parameters p^0, r^0, and ω^0 under initial state X0 and obtain output states, and then calculate its objective function according to (28).
Step 3.3. It is the same as Step 1.3.
Step 3.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution a^0,b^0,k^0; otherwise, let l=l+1 and go to Step 3.3.


## 4.1. Identification of Orderp0

From the first equation of the fractional-order Duffing system (1) and the measurement equation (3) we can see, in the absence of noise, the second component y2(t) of the final output vector Y(t) is the fractional derivative of the first component y1(t) with order p0, that is, Dp0y1=y2.SupposeDpy1=yp, then if and only if p=p0, the energy difference of yp(t) and y2(t) reaches its minimum value. The estimated value of order p0 can be obtained by adjusting the order of p and searching for the minimum value of the energy difference:(12)p^0=argminp∈0,1yp-y22,where x(t)2 denotes the energy of x(t).Therefore, in the first step of our identification method, the objective functionf1 of DE algorithm is chosen as the energy difference of the fractional-order derivative Dpy1 and y2(t): (13)fSDE1p=Dpy1-y22;the order p0 is then estimated by using the DE algorithm.


## 4.2. Identification of Amplituder0 and Angular Frequency ω0

### 4.2.1. Oscillation Variance of the Output of the Fractional-Order Duffing System

For eachp0∈(0,1], there are three equilibrium points of the fractional-order Duffing system (1), (0,0), and (±a0/b0,0).For the equilibrium point(±a0/b0,0), its Jacobian matrix is given by(14)H=01-2a0-k0.Then there are two cases:(1) Ifk02≥8a0, the characteristic values of H are two negative numbers:(15)λ1,2=-k0±k02-8a02,whose arguments are both π.(2) Ifk02<8a0, the characteristic values of H are a pair of conjugate complex numbers:(16)λ1,2=-k0±i8a0-k022,whose arguments are(17)arg⁡λ1=π2+atan⁡k08a0-k02,arg⁡λ2=π+atan⁡8a0-k02k0.Since0<p≤1, then both the two cases satisfy argλi>p/2π, i=1,2. From the conclusions in [26], the equilibrium points (±a0/b0,0) are stable focus points; that is, the phase trajectory from any initial value will eventually be attracted and converge to these two points.For the equilibrium point(0,0), its Jacobian matrix is given by(18)H=01a0-k0;the characteristic values of H are one negative and one positive number λ1,2=-k0±k02+4a0/2. Thus, (0,0) is a saddle point; it is neither an attractor nor a repeller.Similar to the integer-order Duffing system, the outputX(t) of the fractional-order Duffing equation (1) varies as the amplitude r0 of the periodic driving force varies. In the following section, we set the parameters as a0=1, b0=1, k0=1, and ω0=1, and the fractional-order p0 is 0.96. We set the amplitude r0 as 0, 0.1, 0.7, 0.922, 0.933, and 1.8, respectively, and then compare and analyze the first component x1(t) of the output X(t) of the system.(1) Whenr0=0, for any given initial value, the output response x1(t) oscillates around one of the stable focus points and the oscillation amplitude gradually reduces to zero; the system will converge to the phase trajectory of the focus, which reached a steady state. The corresponding phase trajectory and output time response are shown in Figures 1(a) and 1(b), respectively.Figure 1
The phase trajectory and output time response of the fractional-order system withr0=0 and r0=0.1.
(a)
(b)
(c)
(d)Whenr0>0, since the Duffing system is a bistable system, which has two stable focus points and a saddle point, in the presence of the periodic driven force r0cos⁡(ω0t), there is a critical threshold Ac in the bistable system: if 0<r0<Ac, then the output response x1(t) will produce local periodic motion near one steady attractor when the time is long enough. Then x1(t) has the same frequency of the periodic driving force. At this time, the phase plane trajectory shows the attractor under Poincaré mapping. The corresponding phase trajectory and output time response are shown in Figures 1(c) and 1(d), respectively.(2) Whenr0>Ac, x1(t) will create a wide-range transition around the two steady attractors; as the value of r0 increases, the oscillation amplitude of x1(t) also increases, and the phase trajectory will go through four states in turns, namely, “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state,” as shown in Figure 2.Figure 2
The phase trajectory and output time response of the fractional-order system of “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state.”
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It can be observed that the first componentx1(t) of system output around the equilibrium point has increasing oscillation amplitude as the amplitude r0 increases. Note that the horizontal coordinate of the focus point can be regarded as the mean value of x1(t) at time domain, and the variance of x1(t) at time domain exactly characterizes the oscillation degree of the signal around its mean value.Suppose the first component of the output of (1) with initial value X0=(x0,y0) is x1(t); we may assume that its steady state response starts at time T0; then theoscillation variance of x1(t) on [T0,T] is defined by(19)σx=1T-T0∫T0Tx1t-mx2dt,in which, mx=1/T-T0∫T0Tx1(t)dt denotes the mean value of x1(t) in [T0,T].With a fixed angular frequencyω0, the oscillation variance σx is a function of the input amplitude r0. The smaller the oscillation variance of output time response is, the closer the phase state is to the steady state, and the closer the corresponding input signal is to a zero value signal and vice versa. In this paper, we use theoscillation varianceσx of x1(t) to estimate the amplitude of the input periodic driving force.Figure3 shows the oscillation variance σx of fractional-order Duffing system with different orders varies as a function of r0, in which, sampling interval h=0.01 s and sampling time T=20π s.Figure 3
Curves of the oscillation variance of the output of the fractional-order Duffing system as functions of the amplituder0 of the input periodic driving force.From Figure3 we can see, for any fractional-order Duffing system, the oscillation variance σx is an increasing function of absolute value of the amplitude r0 and reaches its minimum value at r0=0. In Section 4.2.2, this dynamic characteristic of fractional-order Duffing system is used to estimate the unknown parameters of its periodic driving force.


### 4.2.2. Identification ofr0 and ω0 Based on the Oscillation Variance

Inputu=-rcos(ωt) into the fractional-order Duffing system (1); we obtained the following controlled system:(20)Dp0x1=x2,Dp0x2=a0x1-b0x13-k0x2+Rcos⁡ω0t+ϕ.Supposer0+Δr=r,ω0+Δω=ω, then(21)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2ϕ=atan⁡r0+Δrsin⁡Δωtr0-r0+Δrcos⁡Δωt.Since1≥cos(Δωt)≥-1 then we have(22)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2≥Δr≥0.Thus, R=0 if and only if (23)Δr=0,Δω=0⟷r=r0,ω=ω0.For anyr, ω, the oscillation variance of the first output component x1(t) of (20) is a function of r, ω:(24)σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt.Then from the above discussion in Section 4.2.1, we know that R=0 if and only if r=r0, ω=ω0, while R=0 if and only if the output is in the steady state, and in this case, σ(r,ω) reaches at its minimum. Thus, the amplitude r0 and angular frequency ω0 can be estimated by searching the minimum value of the output oscillation variance of x1(t):(25)r^0,ω^0=argminr,ωσr,ω,r∈0,rmax,ω∈0,ωmax,in which, rmax and ωmax are the search upper bounds of r and ω, respectively.Therefore, in the second step of our identification method, the objective function of DE algorithm is chosen as the oscillation variance of (20): (26)fSDE2r,ω=σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt;the amplitude r0 and angular frequency ω0 are then estimated by using the DE algorithm.


## 4.2.1. Oscillation Variance of the Output of the Fractional-Order Duffing System

For eachp0∈(0,1], there are three equilibrium points of the fractional-order Duffing system (1), (0,0), and (±a0/b0,0).For the equilibrium point(±a0/b0,0), its Jacobian matrix is given by(14)H=01-2a0-k0.Then there are two cases:(1) Ifk02≥8a0, the characteristic values of H are two negative numbers:(15)λ1,2=-k0±k02-8a02,whose arguments are both π.(2) Ifk02<8a0, the characteristic values of H are a pair of conjugate complex numbers:(16)λ1,2=-k0±i8a0-k022,whose arguments are(17)arg⁡λ1=π2+atan⁡k08a0-k02,arg⁡λ2=π+atan⁡8a0-k02k0.Since0<p≤1, then both the two cases satisfy argλi>p/2π, i=1,2. From the conclusions in [26], the equilibrium points (±a0/b0,0) are stable focus points; that is, the phase trajectory from any initial value will eventually be attracted and converge to these two points.For the equilibrium point(0,0), its Jacobian matrix is given by(18)H=01a0-k0;the characteristic values of H are one negative and one positive number λ1,2=-k0±k02+4a0/2. Thus, (0,0) is a saddle point; it is neither an attractor nor a repeller.Similar to the integer-order Duffing system, the outputX(t) of the fractional-order Duffing equation (1) varies as the amplitude r0 of the periodic driving force varies. In the following section, we set the parameters as a0=1, b0=1, k0=1, and ω0=1, and the fractional-order p0 is 0.96. We set the amplitude r0 as 0, 0.1, 0.7, 0.922, 0.933, and 1.8, respectively, and then compare and analyze the first component x1(t) of the output X(t) of the system.(1) Whenr0=0, for any given initial value, the output response x1(t) oscillates around one of the stable focus points and the oscillation amplitude gradually reduces to zero; the system will converge to the phase trajectory of the focus, which reached a steady state. The corresponding phase trajectory and output time response are shown in Figures 1(a) and 1(b), respectively.Figure 1
The phase trajectory and output time response of the fractional-order system withr0=0 and r0=0.1.
(a)
(b)
(c)
(d)Whenr0>0, since the Duffing system is a bistable system, which has two stable focus points and a saddle point, in the presence of the periodic driven force r0cos⁡(ω0t), there is a critical threshold Ac in the bistable system: if 0<r0<Ac, then the output response x1(t) will produce local periodic motion near one steady attractor when the time is long enough. Then x1(t) has the same frequency of the periodic driving force. At this time, the phase plane trajectory shows the attractor under Poincaré mapping. The corresponding phase trajectory and output time response are shown in Figures 1(c) and 1(d), respectively.(2) Whenr0>Ac, x1(t) will create a wide-range transition around the two steady attractors; as the value of r0 increases, the oscillation amplitude of x1(t) also increases, and the phase trajectory will go through four states in turns, namely, “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state,” as shown in Figure 2.Figure 2
The phase trajectory and output time response of the fractional-order system of “homoclinic orbit state,” “chaotic state,” “chaotic critical state,” and “large-scale periodic state.”
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It can be observed that the first componentx1(t) of system output around the equilibrium point has increasing oscillation amplitude as the amplitude r0 increases. Note that the horizontal coordinate of the focus point can be regarded as the mean value of x1(t) at time domain, and the variance of x1(t) at time domain exactly characterizes the oscillation degree of the signal around its mean value.Suppose the first component of the output of (1) with initial value X0=(x0,y0) is x1(t); we may assume that its steady state response starts at time T0; then theoscillation variance of x1(t) on [T0,T] is defined by(19)σx=1T-T0∫T0Tx1t-mx2dt,in which, mx=1/T-T0∫T0Tx1(t)dt denotes the mean value of x1(t) in [T0,T].With a fixed angular frequencyω0, the oscillation variance σx is a function of the input amplitude r0. The smaller the oscillation variance of output time response is, the closer the phase state is to the steady state, and the closer the corresponding input signal is to a zero value signal and vice versa. In this paper, we use theoscillation varianceσx of x1(t) to estimate the amplitude of the input periodic driving force.Figure3 shows the oscillation variance σx of fractional-order Duffing system with different orders varies as a function of r0, in which, sampling interval h=0.01 s and sampling time T=20π s.Figure 3
Curves of the oscillation variance of the output of the fractional-order Duffing system as functions of the amplituder0 of the input periodic driving force.From Figure3 we can see, for any fractional-order Duffing system, the oscillation variance σx is an increasing function of absolute value of the amplitude r0 and reaches its minimum value at r0=0. In Section 4.2.2, this dynamic characteristic of fractional-order Duffing system is used to estimate the unknown parameters of its periodic driving force.


## 4.2.2. Identification ofr0 and ω0 Based on the Oscillation Variance

Inputu=-rcos(ωt) into the fractional-order Duffing system (1); we obtained the following controlled system:(20)Dp0x1=x2,Dp0x2=a0x1-b0x13-k0x2+Rcos⁡ω0t+ϕ.Supposer0+Δr=r,ω0+Δω=ω, then(21)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2ϕ=atan⁡r0+Δrsin⁡Δωtr0-r0+Δrcos⁡Δωt.Since1≥cos(Δωt)≥-1 then we have(22)R=r02-2r0r0+Δrcos⁡Δωt+r0+Δr2≥Δr≥0.Thus, R=0 if and only if (23)Δr=0,Δω=0⟷r=r0,ω=ω0.For anyr, ω, the oscillation variance of the first output component x1(t) of (20) is a function of r, ω:(24)σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt.Then from the above discussion in Section 4.2.1, we know that R=0 if and only if r=r0, ω=ω0, while R=0 if and only if the output is in the steady state, and in this case, σ(r,ω) reaches at its minimum. Thus, the amplitude r0 and angular frequency ω0 can be estimated by searching the minimum value of the output oscillation variance of x1(t):(25)r^0,ω^0=argminr,ωσr,ω,r∈0,rmax,ω∈0,ωmax,in which, rmax and ωmax are the search upper bounds of r and ω, respectively.Therefore, in the second step of our identification method, the objective function of DE algorithm is chosen as the oscillation variance of (20): (26)fSDE2r,ω=σr,ω=1T-T0∫T0Tx1t-1T-T0∫T0Tx1tdt2dt;the amplitude r0 and angular frequency ω0 are then estimated by using the DE algorithm.


## 4.3. Identification ofa0, b0, and k0

After identification of the orderp0, the amplitude r0, and angular frequency ω0, we obtained the estimated values p^0, r^0, and ω^0, and then the remaining parameters to be identified are a0, b0, and k0. The identification of these three parameters can still be considered as a multidimensional continuous optimization problem. Thus, the approximate response system corresponding to system (3) with noise is(27)Dp^0z1=z2,Dp^0z2=az1-bz13-kz2+r^0cos⁡ω^0t.Therefore, in the last step of our identification method, the objective function of DE algorithm is chosen as the energy difference between the first component of the outputy1(t) in driving system (3) and the first component z1(t) of the output in approximate response system (27):(28)fSDE3a,b,k=y1-z12.The remaining parameters a0, b0, and k0 are then estimated by using the DE algorithm.


## 4.4. Design Steps of the Proposed SDE Method

The whole design steps for parameter identification of fractional-order Duffing system using SDE method can be summarized as follows.Step 1.
Identification of orderp0.
Step 1.1. Choose control parameters of DE, that is, population size NP, mutation factor F, crossover probability CR, and the maximum evolution times L. Randomly initialize NP individuals in the search space as Θ(i,l)=p(i,l), i=1,2,…,NP. Let iteration number l=0.
Step 1.2. Solve the systems (1) and (3) under initial state X0 and obtain output states Y=[y1,y2]T using the numerical method given by (4)-(5). For the ith individual, obtain output states of Dp(i,l)y1=yp from initial state X0, and then calculate its objective function according to (13).
Step 1.3. Perform mutation operation according to (8). Perform crossover operation according to (9) to obtain crossover trial vectors. Calculate the objective function of crossover trial vectors by (13). Perform selection operation according to (11) to generate individuals for next generation.
Step 1.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution p^0; otherwise, let l=l+1 and go to Step 1.3.Step 2.
Identification of amplituder0 and angular frequency ω0.
Step 2.1. Randomly initialize NP individuals as Θ(i,l)=[r(i,l),ω(i,l)], i=1,2,…,NP. Let iteration number l=0.
Step 2.2. For the ith individual, input u=-r(i,l)cos(ω(i,l)t) into the fractional-order Duffing system (1) and get the estimated system (20), obtain output states of (20) from initial state X0, and then calculate its objective function according to (26).
Step 2.3. The same as Step  1.3.
Step 2.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution r^0, ω^0; otherwise, let l=l+1 and go to Step 2.3.Step 3.
Identification of parametersa0, b0, and k0.
Step 3.1. Randomly initialize NP individuals Θ(i,k)=[a(i,l),b(i,l),k(i,l)], i=1,2,…,NP. Let iteration number l=0.
Step 3.2. For the ith individual, solve the approximate response system (27) with estimated parameters p^0, r^0, and ω^0 under initial state X0 and obtain output states, and then calculate its objective function according to (28).
Step 3.3. It is the same as Step 1.3.
Step 3.4. Check if the stopping criterion is met (i.e., l=L). If it is met, then output the optimal solution a^0,b^0,k^0; otherwise, let l=l+1 and go to Step 3.3.


## 5. Numerical Simulations and Analysis

To verify the validity of the proposed sequential parameter identification method based on the DE algorithm (i.e., SDE method), a numerical simulation experiment on parameter identification in a fractional-order Duffing system is performed using the SDE method. The experiment result is compared with the identification results obtained by directly applying the parameter identification method of the DE algorithm (i.e., DE method), and sequential parameter identification method based on PSO (i.e., SPSO method).In the procedure of DE method, suppose the output noisy observation vector of system (1) is Y=[y1(t),y2(t)]T; the estimated system corresponding to system (1) with noise is(29)Dpz1=z2Dpz2=az1-bz13-kz2+rcos⁡ωt.Thus, the state vector of the estimated system is Z=[z1(t),z2(t)]T; the objective function of the DE method with 6 parameters is defined as follows:(30)fDEp,r,ω,a,b,k=y1-z12.The parameters a0, b0, and k0 are then estimated by using the DE algorithm.The SPSO method has the same steps and objective functions as those of the SDE method. The only difference between SPSO and SDE method is that the SPSO method searches optimal parameters of each step by using PSO algorithm, and the SDE method searches optimal parameters by using DE algorithm.The design steps for parameter identification using SDE method, SPSO method, and DE method are shown in Figure4, respectively.Figure 4
Design steps for parameter identification using SDE, SPSO, and DE methods.Following simulation experiments, the sampling steph=0.01. The real value of the order is p0=0.96, and the variation range is (0,1]. The real values of the other parameters are set as a0=1, b0=1, k0=1, r0=1, and ω0=1. The variation range of each parameter is [0,3]. Other main parameters of the SDE, SPSO, and DE methods are listed in Table 1.Table 1
Main parameters of the SDE, SPSO, and DE methods.
Genetic factor
Mutation factor
Maximum iteration
SDE
F
=
0.8
C
R
=
0.9
L
=
300
Learning factor
Inertia weight
Maximum iteration
SPSO
c
1
=
c
2
=
2
w
=
0.6
L
=
300
Genetic factor
Mutation factor
Maximum iteration
DE
F
=
0.8
C
R
=
0.2
L
=
300


### 5.1. Numerical Simulation 1

Firstly, assuming there is no measurement noise, the fractional-order Duffing system is allowed to evolve freely. In order to find the relationship between the identification errors and the population size, we did 9 sets of numerical simulations; the population sizes of three method are all set as 30, 50, and 100, respectively. After the system achieves a transient state, a point is arbitrarily selected as the initial valueX0. From the initial value, the evolution time is set to T = 3 s. We calculate the convergence curves of the logarithmic identification errors (LIE) based on the aforementioned three methods, as shown in Figures 5 and 6. In these two figures, the cutoff of curve denotes that the error has been reduced to zero and the logarithmic results are negative infinity and, thus, cannot be displayed.Figure 5
Evolutionary convergence curves of the logarithmic identification errors (LIE) of the SDE and SPSO methods.
(a)
(b)
(c)Figure 6
Evolutionary convergence curves of the logarithmic identification errors (LIE) of DE method.As can be seen from Figure5, the convergence speeds and the last LIE of the SPSO method and the SDE method with population sizes 30, 50, and 100 are almost the same. It follows that when the population size is greater than 30, the population size has little effect on the final estimation accuracy and convergence speed of the SDE and SPSO methods. Thus, in the following simulation experiments, we only consider the population size of the SDE and SPSO methods as NP = 30.In particular, Figure5(a) shows the convergence curves of the LIE achieved by the SDE and SPSO methods to identify the order p0. The SDE with population size 30 proposed in this study converges to the real value after 48 iterations. Its evolution is highly efficient. By contrast, SPSO with population size 30 converges to the real value after 149 iterations. Figure 5(b) shows the convergence curve of the LIE of the SDE and SPSO methods to identify parameters r0, ω0. SDE with population size 30 converges to the real value after 115 iterations, whereas SPSO with population size 30 converges after 196 iterations. Figure 5(c) shows the convergence curve of the LIE of SDE and SPSO to identify a0, b0, and k0. SDE with population size 30 converges to the real value after 185 iterations, whereas the SPSO with population size 30 converges after 300 iterations to reduce the identification error to a magnitude of 10-31.Figure6 shows the convergence curve of the LIE of DE method to identify all the 6 parameters (p0,r0,ω0,a0,b0,k0); the population sizes of DE method are set as 30, 50, and 100, respectively. The DE method with population size 30 converges after 114 iterations and reduces identification error only to a magnitude of 10-1.11, the DE method with population size 50 converges after 212 iterations and reduces identification error to a magnitude of 10-1.29, and the DE method with population size 100 converges after 178 iterations and reduces identification error to a magnitude of 10-1.32. Thus, the DE method with population size 100 has the best estimation performance and, in the following simulation experiments, we fixed the population size of DE method as NP = 100.In summary, the successive identification methods proposed in this study, that is, SDE and SPSO methods, exhibit higher identification precision and convergence speed than the direct parameter identification method—the DE method. This result can be attributed to the use of the fractional-order Duffing equation dynamic characteristics, thereby achieving successive step by step identification, which significantly reduces identification dimension and, consequently, improves identification speed and accuracy.


### 5.2. Numerical Simulation 2

To obtain more accurate results to compare with the calculation results of the DE and SPSO methods, we also conducted 20 independent numerical experiments using the three methods and use the average value as the final identification result of parameters in Table2. The root mean squared errors (RMSE) of all the parameters of all the three methods are also given in Table 3. Here, the RMSE of parameter θ is given by(31)RMSE  ofθ=∑i=1Nθ^i-θ02N,where θ0 is the true value of θ and θ^i is the ith estimated value of θ, i=1,2,…,N.Table 2
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (without noise).
True values
SDE
SPSO
DE
p
0
0.96
0.960000000000001
0.960000000000001
0.906261312663094
r
0
1
1.000000000000000
1.000000000000000
0.999794594584971
ω
0
1
1.000000000000000
1.000000000000000
1.012465240090026
a
0
1
1.000000000000000
1.000000000000000
1.016663593275141
b
0
1
1.000000000000000
1.000000000000000
1.013125743617082
k
0
1
1.000000000000000
1.000000000000000
1.165684838310568
Cost time (seconds)
109
5015
244Table 3
RMSE comparisons of the calculation results of SDE, SPSO, and DE (without noise).
SDE
SPSO
DE
p
0
0
0
0.229349161461870
r
0
0
0
0.345647703975610
ω
0
0
0
0.285615515196729
a
0
0.25470 × 10−15
0.78505 × 10−15
0.245726407787693
b
0
0
0
0.232148495895464
k
0
0
0
0.664922349600648The main parameters of the SDE, SPSO, and DE methods are the same as numerical simulation 1, the population sizes of the SDE, SPSO methods are both set as 30, and the population size of the DE method is set as 100. It can be observed from Tables2 and 3 that the SDE and SPSO methods can achieve almost entirely accurate identification for all the parameters. The average values and RMSE (except the RMSE of a0) of all the parameters obtained by SDE method are the same as those obtained by the SPSO method, while the RMSE of a0 obtained by the SDE method is a little higher than that of the SPSO method. As seen from the last columns of Tables 2 and 3, simply using the DE algorithm (i.e., the DE method) can also achieve parameter identification. However, the estimation accuracy of the DE method is considerably lower than the aforementioned sequential parameter identification algorithm. In addition, the SDE algorithm has the shortest computation time. The DE consumes twice time as much as the SDE method does, while SPSO consumes the longest time, which is nearly 46 times that of the SDE method.To further test the validity of the algorithm, and considering the effect of noise in practical application results, we add a Gaussian white noise with a standard deviation of 0.05 into (3) and recalculate the average identification results via an independent numerical simulation using the three aforementioned methods after 20 iterations, as shown in Table 4. The RMSE of all the parameters of all the three methods are also given in Table 5. As can be seen from Tables 4 and 5, in the presence of noise, the proposed SDE method and the SPSO method can still achieve parameter identification with high accuracy. The average values of all the parameters obtained by the proposed SDE method are closer to the true values than those obtained by the SPSO method. Compared with the SPSO method, the SDE method can achieve a litter bigger RMSE for the order p0, but lower RMSEs for r0, ω0, and a0,b0,k0, which shows that the stability of PSO is lower than that of the DE algorithm with the increase of identification dimension. However, the SPSO method takes too long a computation time, which is nearly 27 times that of the SDE method. Meanwhile, the accuracy of the identification provided by DE method is considerably lower than those of the SDE and SPSO methods.Table 4
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (with noise,D=0.05).
True values
SDE
SPSO
DE
p
0
0.96
0.960046346271175
0.959492674736478
0.955663495569419
r
0
1
1.002265003070312
1.010175006124644
1.174963453288619
ω
0
1
1.002160216098992
1.004045047242484
0.957180452883223
a
0
1
0.993668488134535
0.983730930793959
1.030124861207625
b
0
1
0.996382591525655
0.984664492992783
1.024923059605036
k
0
1
0.996729257854148
1.029990427269699
1.221970425815988
Cost time (seconds)
182
5086
454Table 5
RMSE comparisons of the calculation results of SDE, SPSO, and DE (with noise,D=0.05).
SDE
SPSO
DE
p
0
0.001584953335274
0.001293305038638
0.023784962594894
r
0
0.017375130670679
0.020862720773417
0.440801581318434
ω
0
0.009230747360439
0.012166424039078
0.120480160254245
a
0
0.037996942936807
0.055568589135587
0.199587022659009
b
0
0.032648336359859
0.048481351301828
0.202914504433050
k
0
0.050175281400924
0.057163873938459
0.748553299332483These results (Tables2–5) also suggest that compared with the SPSO and DE methods, the proposed SDE method has faster convergence speed and better identification accuracy, both in the presence and absence of noise.


## 5.1. Numerical Simulation 1

Firstly, assuming there is no measurement noise, the fractional-order Duffing system is allowed to evolve freely. In order to find the relationship between the identification errors and the population size, we did 9 sets of numerical simulations; the population sizes of three method are all set as 30, 50, and 100, respectively. After the system achieves a transient state, a point is arbitrarily selected as the initial valueX0. From the initial value, the evolution time is set to T = 3 s. We calculate the convergence curves of the logarithmic identification errors (LIE) based on the aforementioned three methods, as shown in Figures 5 and 6. In these two figures, the cutoff of curve denotes that the error has been reduced to zero and the logarithmic results are negative infinity and, thus, cannot be displayed.Figure 5
Evolutionary convergence curves of the logarithmic identification errors (LIE) of the SDE and SPSO methods.
(a)
(b)
(c)Figure 6
Evolutionary convergence curves of the logarithmic identification errors (LIE) of DE method.As can be seen from Figure5, the convergence speeds and the last LIE of the SPSO method and the SDE method with population sizes 30, 50, and 100 are almost the same. It follows that when the population size is greater than 30, the population size has little effect on the final estimation accuracy and convergence speed of the SDE and SPSO methods. Thus, in the following simulation experiments, we only consider the population size of the SDE and SPSO methods as NP = 30.In particular, Figure5(a) shows the convergence curves of the LIE achieved by the SDE and SPSO methods to identify the order p0. The SDE with population size 30 proposed in this study converges to the real value after 48 iterations. Its evolution is highly efficient. By contrast, SPSO with population size 30 converges to the real value after 149 iterations. Figure 5(b) shows the convergence curve of the LIE of the SDE and SPSO methods to identify parameters r0, ω0. SDE with population size 30 converges to the real value after 115 iterations, whereas SPSO with population size 30 converges after 196 iterations. Figure 5(c) shows the convergence curve of the LIE of SDE and SPSO to identify a0, b0, and k0. SDE with population size 30 converges to the real value after 185 iterations, whereas the SPSO with population size 30 converges after 300 iterations to reduce the identification error to a magnitude of 10-31.Figure6 shows the convergence curve of the LIE of DE method to identify all the 6 parameters (p0,r0,ω0,a0,b0,k0); the population sizes of DE method are set as 30, 50, and 100, respectively. The DE method with population size 30 converges after 114 iterations and reduces identification error only to a magnitude of 10-1.11, the DE method with population size 50 converges after 212 iterations and reduces identification error to a magnitude of 10-1.29, and the DE method with population size 100 converges after 178 iterations and reduces identification error to a magnitude of 10-1.32. Thus, the DE method with population size 100 has the best estimation performance and, in the following simulation experiments, we fixed the population size of DE method as NP = 100.In summary, the successive identification methods proposed in this study, that is, SDE and SPSO methods, exhibit higher identification precision and convergence speed than the direct parameter identification method—the DE method. This result can be attributed to the use of the fractional-order Duffing equation dynamic characteristics, thereby achieving successive step by step identification, which significantly reduces identification dimension and, consequently, improves identification speed and accuracy.


## 5.2. Numerical Simulation 2

To obtain more accurate results to compare with the calculation results of the DE and SPSO methods, we also conducted 20 independent numerical experiments using the three methods and use the average value as the final identification result of parameters in Table2. The root mean squared errors (RMSE) of all the parameters of all the three methods are also given in Table 3. Here, the RMSE of parameter θ is given by(31)RMSE  ofθ=∑i=1Nθ^i-θ02N,where θ0 is the true value of θ and θ^i is the ith estimated value of θ, i=1,2,…,N.Table 2
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (without noise).
True values
SDE
SPSO
DE
p
0
0.96
0.960000000000001
0.960000000000001
0.906261312663094
r
0
1
1.000000000000000
1.000000000000000
0.999794594584971
ω
0
1
1.000000000000000
1.000000000000000
1.012465240090026
a
0
1
1.000000000000000
1.000000000000000
1.016663593275141
b
0
1
1.000000000000000
1.000000000000000
1.013125743617082
k
0
1
1.000000000000000
1.000000000000000
1.165684838310568
Cost time (seconds)
109
5015
244Table 3
RMSE comparisons of the calculation results of SDE, SPSO, and DE (without noise).
SDE
SPSO
DE
p
0
0
0
0.229349161461870
r
0
0
0
0.345647703975610
ω
0
0
0
0.285615515196729
a
0
0.25470 × 10−15
0.78505 × 10−15
0.245726407787693
b
0
0
0
0.232148495895464
k
0
0
0
0.664922349600648The main parameters of the SDE, SPSO, and DE methods are the same as numerical simulation 1, the population sizes of the SDE, SPSO methods are both set as 30, and the population size of the DE method is set as 100. It can be observed from Tables2 and 3 that the SDE and SPSO methods can achieve almost entirely accurate identification for all the parameters. The average values and RMSE (except the RMSE of a0) of all the parameters obtained by SDE method are the same as those obtained by the SPSO method, while the RMSE of a0 obtained by the SDE method is a little higher than that of the SPSO method. As seen from the last columns of Tables 2 and 3, simply using the DE algorithm (i.e., the DE method) can also achieve parameter identification. However, the estimation accuracy of the DE method is considerably lower than the aforementioned sequential parameter identification algorithm. In addition, the SDE algorithm has the shortest computation time. The DE consumes twice time as much as the SDE method does, while SPSO consumes the longest time, which is nearly 46 times that of the SDE method.To further test the validity of the algorithm, and considering the effect of noise in practical application results, we add a Gaussian white noise with a standard deviation of 0.05 into (3) and recalculate the average identification results via an independent numerical simulation using the three aforementioned methods after 20 iterations, as shown in Table 4. The RMSE of all the parameters of all the three methods are also given in Table 5. As can be seen from Tables 4 and 5, in the presence of noise, the proposed SDE method and the SPSO method can still achieve parameter identification with high accuracy. The average values of all the parameters obtained by the proposed SDE method are closer to the true values than those obtained by the SPSO method. Compared with the SPSO method, the SDE method can achieve a litter bigger RMSE for the order p0, but lower RMSEs for r0, ω0, and a0,b0,k0, which shows that the stability of PSO is lower than that of the DE algorithm with the increase of identification dimension. However, the SPSO method takes too long a computation time, which is nearly 27 times that of the SDE method. Meanwhile, the accuracy of the identification provided by DE method is considerably lower than those of the SDE and SPSO methods.Table 4
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (with noise,D=0.05).
True values
SDE
SPSO
DE
p
0
0.96
0.960046346271175
0.959492674736478
0.955663495569419
r
0
1
1.002265003070312
1.010175006124644
1.174963453288619
ω
0
1
1.002160216098992
1.004045047242484
0.957180452883223
a
0
1
0.993668488134535
0.983730930793959
1.030124861207625
b
0
1
0.996382591525655
0.984664492992783
1.024923059605036
k
0
1
0.996729257854148
1.029990427269699
1.221970425815988
Cost time (seconds)
182
5086
454Table 5
RMSE comparisons of the calculation results of SDE, SPSO, and DE (with noise,D=0.05).
SDE
SPSO
DE
p
0
0.001584953335274
0.001293305038638
0.023784962594894
r
0
0.017375130670679
0.020862720773417
0.440801581318434
ω
0
0.009230747360439
0.012166424039078
0.120480160254245
a
0
0.037996942936807
0.055568589135587
0.199587022659009
b
0
0.032648336359859
0.048481351301828
0.202914504433050
k
0
0.050175281400924
0.057163873938459
0.748553299332483These results (Tables2–5) also suggest that compared with the SPSO and DE methods, the proposed SDE method has faster convergence speed and better identification accuracy, both in the presence and absence of noise.


## 6. Conclusion

This study first analyzes the stability of equilibrium points in a fractional-order Duffing system and then combines the dynamic characteristics of this equation to establish a sequential parameter identification method based on the DE algorithm. First, on the basis of the relationship between the two output states of the fractional-order Duffing system, identification of fractional-order is estimated using DE algorithm. Then, the relationship between system output oscillation variance and amplitude of the periodic driving force is considered and combined with the DE algorithm to identify unknown parameters of the periodic driving force. The remaining parameters are then identified using DE algorithm. The introduction of the sequential identification method gradually decreases identification dimension (in the three-step process, identification dimensions are 1, 2, and 3, resp., and the total identification dimension is 6), which significantly improves identification speed and accuracy.Through a simulation experimental analysis, the SDE method with successive step by step identification steps is proposed to obtain optimal performance. This method is stable. The identification accuracy of the SPSO method with successive step by step identification steps ranks second. However, the SPSO method costs a lot of computation time because of the PSO searching algorithm. The identification accuracy of the DE method which directly applies the parameter identification is quite lower than those of the other two methods. In addition, the comparison of the simulations of the DE and PSO algorithm shows that although PSO demonstrates rapid convergence during the early stage, it takes obviously more computation time to obtain the same identification accuracy as DE algorithm.



---

*Source: 3572365-2017-09-28.xml*