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

Using the dynamic properties of fractional-orderDuffing system, a sequential parameter identificationmethod 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.


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 fractionalorder 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][4][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][8][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 integerorder 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][22][23][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.

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]: where   0  is the fractional-order derivative based on the definition of Caputo [2]: in which, Γ(⋅) is the gamma function and Γ() = ∫ +∞ 0  −1  −1 d.When  0 = 1, (1) is the classical Duffing equation.
Here,  0 ,  0 ,  0 , amplitude  0 , angular frequency  0 , and order  0 are unknown parameters to be identified.In this paper, we suppose these parameters are all positive real numbers and  0 ∈ (0, 1]. The numerical algorithm for solving the fractional-order Duffing equation ( 1) is based on the numerical algorithm given in [23]: Mathematical Problems in Engineering 3 where

Differential Evolution Algorithm
Consider a general optimization problem: where  is the objective function, Θ = [ 1 ,  2 , . . .,   ] is the -dimensional parameter vector, and 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 , and then the initial value of the th individual is in which,   (, 0) = Min   + [Max   − Min   ] × rand(1),  = 1, 2, . . ., , and rand(1) denotes random numbers generated by a uniform distribution.At the same time, set up the mutation factor  ∈ [0, 2], crossover probability CR, and maximum evolution times .
(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 Θ(, ), the mutation vector is produced by the following difference strategy: where the random selected numbers  1 ,  2 , and  3 are not the same,  1 ,  2 ,  3 , and the sequence number  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 th generation and its mutation population, th generation test vector was introduced by the crossover operation in which,

Parameter Identification of Fractional-Order Duffing System
Based on the dynamic characteristics of the fractionalorder Duffing system, the parameter identification method proposed in this paper will be divided into three steps: Step 1. Identification of order  0 ; Step 2. Identification of amplitude  0 and angular frequency  0 ; Step 3. Identification of  0 ,  0 , and  0 .
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 Order  0 .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  2 () of the final output vector Y() is the fractional derivative of the first component  1 () with order  0 , that is, Suppose    1 =   , then if and only if  =  0 , the energy difference of   () and  2 () reaches its minimum value.The estimated value of order  0 can be obtained by adjusting the order of  and searching for the minimum value of the energy difference: where ‖()‖ 2 denotes the energy of ().Therefore, in the first step of our identification method, the objective function  1 of DE algorithm is chosen as the energy difference of the fractional-order derivative    1 and  2 (): the order  0 is then estimated by using the DE algorithm.For each  0 ∈ (0, 1], there are three equilibrium points of the fractional-order Duffing system (1), (0, 0), and (±√ 0 / 0 , 0).

Identification of
For the equilibrium point (±√ 0 / 0 , 0), its Jacobian matrix is given by Then there are two cases: (1) If  0 2 ≥ 8 0 , the characteristic values of  are two negative numbers: whose arguments are both .
(2) If  0 2 < 8 0 , the characteristic values of  are a pair of conjugate complex numbers: whose arguments are Since 0 <  ≤ 1, then both the two cases satisfy |arg(  )| > (/2),  = 1, 2. From the conclusions in [26], the equilibrium points (±√ 0 / 0 , 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 the characteristic values of  are one negative and one positive number  1,2 = [− 0 ± √  0 2 + 4 0 ]/2.Thus, (0, 0) is a saddle point; it is neither an attractor nor a repeller.Similar to the integer-order Duffing system, the output X() of the fractional-order Duffing equation (1) varies as the amplitude  0 of the periodic driving force varies.In the following section, we set the parameters as  0 = 1,  0 = 1,  0 = 1, and  0 = 1, and the fractional-order  0 is 0.96.We set the amplitude  0 as 0, 0.1, 0.7, 0.922, 0.933, and 1.8, respectively, and then compare and analyze the first component  1 () of the output X() of the system.
(1) When  0 = 0, for any given initial value, the output response  1 () 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.
When  0 > 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  0 cos( 0 ), there is a critical threshold   in the bistable system: if 0 <  0 <   , then the output response  1 () will produce local periodic motion near one steady attractor when the time is long enough.Then  1 () 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) When  0 >   ,  1 () will create a wide-range transition around the two steady attractors; as the value of  0 increases, the oscillation amplitude of  1 () 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.
It can be observed that the first component  1 () of system output around the equilibrium point has increasing oscillation amplitude as the amplitude  0 increases.Note that the horizontal coordinate of the focus point can be regarded as the mean value of  1 () at time domain, and the variance of  1 () 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  0 = ( 0 ,  0 ) is  1 (); we may assume that its steady state response starts at time  0 ; then the oscillation variance of  1 () on [ 0 , ] is defined by  With a fixed angular frequency  0 , the oscillation variance   is a function of the input amplitude  0 .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 the oscillation variance   of  1 () to estimate the amplitude of the input periodic driving force.
Figure 3 shows the oscillation variance   of fractionalorder Duffing system with different orders varies as a function of  0 , in which, sampling interval ℎ = 0.01 s and sampling time  = 20 s.
From Figure 3 we can see, for any fractional-order Duffing system, the oscillation variance   is an increasing function of absolute value of the amplitude  0 and reaches its minimum value at  0 = 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 of  0 and  0 Based on the Oscillation Variance.Input  = −cos() into the fractional-order Duffing system (1); we obtained the following controlled system: Suppose  0 + Δ = ,  0 + Δ = , then Since 1 ≥ cos(Δ) ≥ −1 then we have Thus,  = 0 if and only if For any , , the oscillation variance of the first output component  1 () of ( 20) is a function of , : 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."Then from the above discussion in Section 4.2.1, we know that  = 0 if and only if  =  0 ,  =  0 , while  = 0 if and only if the output is in the steady state, and in this case, (, ) reaches at its minimum.Thus, the amplitude  0 and angular frequency  0 can be estimated by searching the minimum value of the output oscillation variance of  1 (): r0 , ω0 = arg min ,  (, ) , in which,  max and  max are the search upper bounds of  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): the amplitude  0 and angular frequency  0 are then estimated by using the DE algorithm.

4.3.
Identification of  0 ,  0 , and  0 .After identification of the order  0 , the amplitude  0 , and angular frequency  0 , we obtained the estimated values p0 , r0 , and ω0 , and then the remaining parameters to be identified are  0 ,  0 , and  0 .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 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 output  1 () in driving system (3) and the first component  1 () of the output in approximate response system (27): The remaining parameters  0 ,  0 , and  0 are then estimated by using the DE algorithm.

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.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.,  = ).If it is met, then output the optimal solution p0 ; otherwise, let  =  + 1 and go to Step 1.3.
Step 2. Identification of amplitude  0 and angular frequency  0 .
Step 3.2.For the th individual, solve the approximate response system (27) with estimated parameters p0 , r0 , and ω0 under initial state  0 and obtain output states, and then calculate its objective function according to (28).
Step 3.3.It is the same as Step 1.3.

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 = [ 1 (),  2 ()]  ; the estimated system corresponding to system (1) with noise is Thus, the state vector of the estimated system is  = [ 1 (),  2 ()]  ; the objective function of the DE method with 6 parameters is defined as follows: The parameters  0 ,  0 , and  0 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 Figure 4, respectively.
Following simulation experiments, the sampling step ℎ = 0.01.The real value of the order is  0 = 0.96, and the variation range is (0, 1].The real values of the other parameters are set as  0 = 1,  0 = 1,  0 = 1,  0 = 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.

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 value  0 .From the initial value, the evolution time is set to  = 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.
As can be seen from Figure 5, 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, Figure 5(a) shows the convergence curves of the LIE achieved by the SDE and SPSO methods to identify the order  0 .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. 0 ,  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  0 ,  0 , and  0 .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 .
Figure 6 shows the convergence curve of the LIE of DE method to identify all the 6 parameters ( 0 ,  0 ,  0 ,  0 ,  0 ,  0 ); 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 fractionalorder Duffing equation dynamic characteristics, thereby achieving successive step by step identification, which significantly reduces identification dimension and, consequently, improves identification speed and accuracy.

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 Table 2.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 where  0 is the true value of  and θ is the th estimated value of ,  = 1, 2, . . ., .
The 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 Tables 2 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  0 ) of all the parameters obtained by SDE method are the same as those obtained by the SPSO method, while the RMSE of  0 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  0 , but lower RMSEs for  0 ,  0 , and  0 ,  0 ,  0 , 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.
These results (Tables 2-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.

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

Figure 1 :
Figure 1: The phase trajectory and output time response of the fractional-order system with  0 = 0 and  0 = 0.1.

Figure 3 :
Figure 3: Curves of the oscillation variance of the output of the fractional-order Duffing system as functions of the amplitude  0 of the input periodic driving force.

Figure 4 :
Figure 4: Design steps for parameter identification using SDE, SPSO, and DE methods.

Figure 5 :
Figure 5: Evolutionary convergence curves of the logarithmic identification errors (LIE) of the SDE and SPSO methods.
Figure 5(b)  shows the convergence curve of the LIE of the SDE and SPSO methods to identify parameters = 100 DE, size = 50 DE, size = 30

Figure 6 :
Figure 6: Evolutionary convergence curves of the logarithmic identification errors (LIE) of DE method.

)
(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: Amplitude  0 and Angular Frequency  0 4.2.1.Oscillation Variance of the Output of the Fractional-Order Duffing System.

Table 1 :
Main parameters of the SDE, SPSO, and DE methods.

Table 2 :
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (without noise).

Table 3 :
RMSE comparisons of the calculation results of SDE, SPSO, and DE (without noise).

Table 4 :
Average value and cost time comparison of the calculation results of SDE, SPSO, and DE (with noise,  = 0.05).