1 Introduction

A fractional order model is a system that is characterized by a fractional differential equation containing derivatives of non-integer order. Several engineering problems, such as transmission lines [¹], electrical circuits [²] and control systems [³], can be more accurately described by fractional differential equations than integer order schemes. For this reason, in the last decade, the fractional order systems [⁴–⁸] have attracted the interests of several research communities.

System identification is a practical way to model a fractional order system. However, because the mathematical interpretation of fractional calculus is lightly distinct to integer calculus, it is difficult to model real fractional order systems directly based on analytic mechanisms [⁹]. For classical integer order system, once the maximum order of the system has been defined, the parameters of the model can be identified directly. However, for a fractional order system, because identification requires the choice of the fractional order of the operators, and the systematic parameters, the identification process of such systems is more complex than that of the integer order models [¹⁰]. Under such conditions, most of the classical identification methods cannot directly applied to identification of a fractional order systems [¹¹].

The problem of estimating the parameters of fractional order systems has been commonly solved through the use of deterministic methods such as non-linear optimization techniques [¹²], input output frequency contents [¹³] or operational matrix [¹⁴]. These methods have been exhaustively analyzed and represent the most consolidated available tools. The interested reader in such approaches can be referred to [¹⁵] for a recent survey on the state-of-the-art.

As an alternative to classical techniques, the problem of identification in fractional order systems has also been handled through evolutionary methods. In general, they have demonstrated, under several circumstances, to deliver better results than those based on deterministic approaches in terms of accuracy and robustness [¹⁶]. Under these methods, an individual is represented by a candidate model. Just as the evolution process unfolds, a set of evolutionary operators are applied in order to produce better individuals.

The quality of each candidate solution is evaluated through an objective function whose final result represents the affinity between the estimated model and the actual one. Some examples of these approaches used in the identification of fractional order systems involve methods such as Genetic Algorithms (GA) [¹⁷], Artificial Bee Colony (ABC) [¹⁸], Differential Evolution (DE) [¹⁹] and Particle Swarm Optimization (PSO) [²⁰]. Although these algorithms present interesting results, they have an important limitation:

They frequently obtain sub-optimal solutions as a consequence of the limited balance between exploration and exploitation in their search strategies. This limitation is associated to their evolutionary operators employed to modify the individual positions. In such algorithms, during their operation, the position of each individual for the next iteration is updated producing an attraction towards the position of the best particle seen so-far or towards other promising individuals. Therefore, as the algorithm evolves, such behaviors cause that the entire population rapidly concentrates around the best particles, favoring the premature convergence and damaging the appropriate exploration of the search space [²¹, ²²].

This paper presents an algorithm for parameter identification of fractional-order chaotic systems. In order to determine the parameters, the proposed method uses a novel evolutionary method called Locust Search (LS) [²³, ³¹, ³²] which is based on the behavior of swarms of locusts. In the proposed algorithm, individuals emulate a group of locusts which interact to each other based on the biological laws of the cooperative swarm.

The algorithm considers two different behaviors: solitary and social. Depending on the behavior, each individual is conducted by a set of evolutionary operators which mimics different cooperative conducts that are typically found in the swarm. Different to most of existent evolutionary algorithms, the behavioral model in the proposed approach explicitly avoids the concentration of individuals in the current best positions. Such fact allows avoiding critical flaws such as the premature convergence to sub-optimal solutions and the incorrect exploration-exploitation balance. Numerical simulations have been conducted on the fractional-Order Van der Pol oscillator to show the effectiveness of the proposed scheme.

The paper is organized as follows. In Section 2, the concepts of fractional calculus are introduced. Section 3 gives a description for the Locust Search algorithm. Section 4 gives a brief description of the fractional-order Van der Pol Oscillator. Section 5 formulates the parameter estimation problem. Section 6 shows the experimental results. Finally some conclusions are discussed in Section 7.

2 Fractional Calculus

Fractional calculus is a generalization of integration and differentiation to non-integer order fundamental operator. The differential-integral operator, denoted by takes both the fractional derivative and the fractional integral in a single expression which is defined as:

(1)

where a and t represents the operation bounds whereas q ∈ ℜ. The commonly used definitions for fractional derivatives are the Grünwald-Letnikov, Riemann-Liouville [²⁴] and Caputo [²⁵]. According to the Grünwald-Letnikov approximation, the fractional-order derivative of order q is defined as follows:

(2)

In the numerical calculation of fractional-order derivatives, the explicit numerical approximation of the q-th derivative at the points kh, (k = 1, 2,...) maintains the following form [²⁶]:

(3)

where L_m is the memory length t_k = kh, h, is the time step and are the binomial coefficients. For their calculation we can use the following expression:

(4)

Then, the general numerical solution of the fractional differential equation is defined as follows:

(5)

3 Locust Search (LS) Algorithm

In the operation of LS [²³], a population of N locusts (individuals) is processed from the initial stage (k=0) to a total gen number iterations (k= gen). Each individual symbolizes an n-dimensional vector where each dimension represents a domain variable of the optimization problem to be solved. The set of variables represents the valid search space , where lb_d and ub_d represents the lower and upper bounds for the d dimension, respectively. The quality of each element (candidate solution) is evaluated by using the objective function . In LS, at each iteration consists of two operators: (A) solitary and (B) social.

3.1 Solitary Operation (A)

In the solitary operation, a new location p_i (i ∈ [1,..., N]) is generated by modifying the current element location with a change of position . is the result of the individual interactions experimented by as a consequence of its biological behavior. Such interactions are pairwise computed among and the other N-1 individuals in the swarm. Therefore, the final force exerted between and is computed by considering the following model:

(6)

where is the unit-vector, pointing from to . Furthermore, rand(1,-1) is an number randomly produced between 1 and -1. The factor s(r_ij) represents the social relation between and , which is calculated as follows:

(7)

Here, r_ij is the distance between and , F represents the strength of attraction whereas L is the attractive length factor. It is assumed that F<1 and L>1 so that repulsion is stronger in a shorter-scale, while attraction is applied in a weaker and longer-scale. is a function that calculates the dominance value of the best element between and . In order to operate , all the individuals from are arranged in terms of their fitness values. Therefore, a rank is assigned to each element, so that the best individual obtains the rank 0 (zero) whereas the worst individual receives the rank N-1. Under such conditions, the function is defined as follows:

(8)

where rank(α) delivers the rank of the α-element. According to Eq. 8, gives a value within [0,1]. Fig. 1 shows the behavior of considering 100 elements. In the Figure, it is assumed that represents one of the 99 individuals with ranks among 0 and 98 whereas is fixed to the worst individual (rank 99).

Fig. 1 Behavior of considering 100 individuals 

Then, the resultant force on each element is computed as the superposition of all of the pairwise interactions exerted on it:

(9)

Finally, Δl_i is assumed similar to the social force experimented by as the superposition of all of the pairwise reciprocal forces. Consequently, Δl_i is represented as follows:

(10)

After calculating the new locations of the population , the final locations must be computed. This procedure can be summarized by the following formulation (in terms of a minimization problem):

(11)

3.2 Social Operation (B)

The social operation is a discriminating operation which considers only to a subset E of the final positions F (where E ⊆ F). In the process first is necessary to order F in terms of their fitness values and collect the individuals in a temporal population B = {b₁, b₂,..., b_N}. The individuals of B are arranged so that the best element is located in the first position whereas the worst individual is situated in the last location b_N. Under such conditions, E is composed by the first g position of B (the best elements). Then, a subspace C_j is defined around each selected element f_j ∈ E. The size of C_j depends on the distance e_d which is determined as follows:

(12)

where ub_q and lb_q are the upper and lower limits of the q-th dimension, n is the number of dimensions of the optimization problem, whereas β ∈ [0,1] is a tuning factor. Therefore, the bounds of C_j are modeled as follows:

(13)

where and are the upper and lower limits of the q-th-dimension for the subspace C_j, respectively. Once creating the subspace C_j in the neighborhood of the element fj ∈ E, a set of h new elements are randomly produced within the limits defined by Eq. 13. Considering the h samples, the new individual of the next population L^k+1 must be extracted. In order to select , the best element , in terms of fitness value from the h samples (where ), is examined. If is better than f_j according to their fitness values, is updated with , otherwise the position of f_j is assigned to . The elements of F that have not been considered by the procedure (f_w ∉ E) transport their corresponding values to L^k+1 without variation.

The social operation is used to exploit only favorable solutions. According to the social operation, inside each subspace C_j, h random samples are produced. Since the number of selected elements in each subspace is very small (typically h < 4), the use of this operator cannot be considered computational expensive.

4 Fractional-order Van der Pol Oscillator

The Van der Pol Oscillator model has been extensively studied as a complex example of non-linear system. It provides important models for a wide range of dynamic behaviors for several engineering applications [²⁶, ²⁷]. The classical integer-order Van der Pol Oscillator is described by a second-order non-linear differential equation as follows:

(14)

where ε is a control parameter that reflects the nonlinearity degree of the system. On the other hand, the fractional-order Van der Pol Oscillator model of order q is defined by the following formulation [²⁸]:

(15)

Considering the Grünwald-Letnikov approximation (see Eq. 5), the numerical solution for the fractional-order Van der Pol Oscillator is given by:

(16)

5 Problem Formulation

In the proposed approach, the identification process is considered as a multidimensional optimization problem. In the optimization process, the parameters of a new fractional-order chaotic system FOC_E are determined by using the LS method from the operation of the original fractional-order chaotic system FOC_O.

The idea is that FOC_E presents the best possible parametric affinity with FOC_O. Under such circumstances, the original fractional-order chaotic system FOC_O can be defined as follows:

(17)

here Y = [y₁, y₂, ..., y_m]^T denotes the state vector of the system, Y₀ symbolizes the initial state vector, θ = [θ₁, θ₂, ..., θ_m]^T represents the original systematic parameter set, q = [q₁, q₂, ..., q_m]^T for 0 < q_i < 1 (i ∈ [1,..., m]) corresponds to the fractional derivative orders and F is a generic non-linear function. On the other hand, the estimated fractional-order chaotic system FOC_E can be modeled as follows:

(18)

where , and denotes the estimated state system, the estimated systematic parameter vector and the estimated fractional orders, respectively. Since the goal is that FOC_E presents the best possible parametric affinity with FOC_O, the problem can be approached as an optimization problem described by the following formulation:

(19)

where denotes the best possible parametric values obtained by the optimization process, Ω symbolizes the search space admitted for parameters ( and ) whereas J represents the objective function that evaluates the parametric affinity between FOC_O and FOC_E. This affinity can be computed as follows:

(20)

where and represent the state values produced by the original and estimated systems, respectively. On the other hand, k denotes the sampling time point and M represents the length of data used for parameter estimation. According to the optimization problem formulated in Eq. 19, the parameter identification can be achieved by searching suitable values of and within the searching space Ω, such that the objective function has been minimized.

Fig. 2 shows the graphic representation of the identification process. Since the fractional-order Van der Pol oscillator has been chosen to test the performance of the proposed approach, the fractional-order system maintain two different fractional derivative orders q = [q₁, q₂]^T (m=2) and one systematic parameter ε.

Fig. 2 Evolutionary Algorithm for Fractional-order System Parameter Estimation 

6 Experimental Results

To verify the effectiveness and robustness of the proposed approach, the fractional-order Van der Pol oscillator is chosen to test its performance. The simulations has been conducted by using MATLAB (Version 7.1, MathWorks, Natick, MA, USA) on an Intel(R) Core(TM) i7-3470 CPU, 3.2 GHz with 4 GB of RAM. In order to calculate the objective function, the number of samples is set as 300 and the step size is 0.01.

In this section, the results of the LS algorithm have been compared to those produced by the Genetic Algorithms (GA) [¹⁷], Particle Swarm Optimization (PSO) method [²⁰], the Differential Evolution (DE) [¹⁹], and the proposed method. In all comparisons, the population has been set to 40 (N=40) individuals. The maximum iteration number for all functions has been set to 100. Such stop criterion has been selected to maintain compatibility to similar works reported in the literature [¹⁶].

The parameter setting for each of the algorithms in the comparison is described as follows:

In the experiments, the fractional-order Van der Pol Oscillator to be estimated has been configured such that q₁ = 1.2, q₂ = 0.8 and ε = 1. Similarly, the initial state has been set to [0.02- 0.2].

The statistical results of the best, the mean and the worst estimated parameters with the corresponding relative error values over 100 independent runs are shown in Table 1. From Table 1, it can be easily seen that the estimated values generated by the proposed LS algorithm are closer to the actual parameter values, which means that it is more accurate than the standard GA, PSO and DE algorithms.

Table 1 Simulation result of the algorithms GA, PSO, DE and LS 

Likewise, it can also be clearly found that the relative error values obtained by the LS algorithm are all smaller than those of the standard GA, PSO and DE algorithms, which can also prove that the LS algorithm has a higher performance in terms of accuracy. Therefore, the estimated parameters can be closer to the true values than the GA, PSO and DE algorithms.

With this evidence, it can be concluded that the LS algorithm can more efficiently identify a fractional-order systems than the other algorithms used in the comparisons. Therefore, the estimated parameters can be closer to the true values than the GA, PSO and DE algorithms.

With this evidence, it can be concluded that the LS algorithm can more efficiently identify a fractional-order systems than the other algorithms used in the comparisons. In order to show the proficiency, of the proposed approach, Figure 3 presents the phase diagrams of the Van der Pol Oscillator by using the mean estimated parameters for each method.

Fig. 3 Phase diagrams of the Van der Pol Oscillator by using the mean estimated parameters for (a) GA, (b) PSO, (c) DE and (d) the proposed approach 

The convergence curves of the parameters and fitness values estimated by the set of algorithms are shown in Figures 4-6 in a single execution. From Figures 4-6, it can be clearly observed that convergence processes of the parameters and fitness values of LS algorithm are better than other algorithms. Additionally, the estimated parameter values obtained by the LS algorithm fall faster than the other algorithms.

Fig. 4 Estimated parameter q₁ (fractional order). 

Fig. 5 Estimated parameter q₂ (fractional order). 

Fig. 6 Estimated systematic parameter ε. 

Furthermore, Table 2 shows the average best solution obtained by each algorithm. The average best solution (ABS) expresses the average value of the best function evaluations that have been obtained from 100 independent executions. A non-parametric statistical significance test known as the Wilcoxon’s rank sum test for independent samples [³⁰, ³¹] has been conducted with an 5% significance level, over the “average best-solution” data of Table 2.

Table 2 Average best solution obtained by each algorithm GA, PSO, DE and LS 

Table 3 reports the p-values produced by Wilcoxon’s test for the pair-wise comparison of the “average best-solution” of two groups. Such groups are formed by LS vs. GA, LS vs. PSO and LS vs. DE. As a null hypothesis, it is assumed that there is no significant difference between mean values of the two algorithms.

Table 3 p-values produced by Wilcoxon’s test that compares LS vs GA, LS vs PSO and DE over the “average best-solution” values from Table 3  

The alternative hypothesis considers a significant difference between the “average best-solution” values of both approaches. All p-values reported in the table are less than 0.05 (5% significance level) which is a strong evidence against the null hypothesis, indicating that the LS results are statistically significant and that it has not occurred by coincidence (i.e. due to the normal noise contained in the process).

On the other hand, the time spent by all methods in the estimation of the parameter set is evaluated. Evolutionary methods are, in general, complex pieces of software with several operators and stochastic branches. Therefore, it is difficult to conduct a complexity analysis from a deterministic perspective. Therefore, the number of function evaluations (NFE) is commonly used in order to evaluate the computational effort. Table 4 presents the NFE values obtained by each approach. The Table reports the averaged value considering 100 different executions. From the table, it can be seen that the proposed LS method present the best performance.

Table 4 Average Number of Function Evaluations (NFE) obtained by each algorithm GA, PSO, DE and LS 

7 Conclusions

Due to its multiple applications, parameter identification for fractional-order chaotic systems has attracted the interests of several research communities. In the identification, the parameter estimation process is transformed into a multidimensional optimization problem where fractional orders, as well as functional parameters of the chaotic system are considered the decision variables. Under this approach, the complexity of fractional-order chaotic systems tends to produce multimodal error surfaces for which their cost functions are significantly difficult to minimize. Several algorithms based on evolutionary computation principles have been successfully applied to identify the parameters of fractional-order chaotic systems. However, most of them maintain an important limitation, they frequently obtain sub-optimal results as a consequence of an inappropriate balance between exploration and exploitation in their search strategies.

In this paper, an algorithm for parameter identification of fractional-order chaotic systems has been presented. In order to determine the parameters, the proposed method uses a novel evolutionary method called Locust Search (LS) [R1] which is based on the behavior of swarms of locusts. In the proposed algorithm, individuals emulate a group of locusts which interact to each other based on the biological laws of the cooperative swarm.

The algorithm considers two different behaviors: solitary and social. Depending on the behavior, each individual is conducted by a set of evolutionary operators which mimics different cooperative conducts that are typically found in the swarm. Different to most of existent evolutionary algorithms, the behavioral model in the proposed approach explicitly avoids the concentration of individuals in the current best positions. Such fact allows to avoid critical flaws such as the premature convergence to sub-optimal solutions and the incorrect exploration-exploitation balance.

In order to test the proficiency and robustness of the presented method, it has been compared to other algorithms based on evolutionary principles such as GA, PSO and DE. The comparison examines the identification of the fractional Van der Pol Oscillator. The results show a high performance of the proposed descriptor in terms of precision and robustness.