Genetic Algorithm-Based Identification of Fractional-Order Systems

Fractional calculus has become an increasingly popular tool for modeling the complex behaviors of physical systems from diverse domains. One of the key issues to apply fractional calculus to engineering problems is to achieve the parameter identification of fractional-order systems. A time-domain identification algorithm based on a genetic algorithm (GA) is proposed in this paper. The multi-variable parameter identification is converted into a parameter optimization by applying GA to the identification of fractional-order systems. To evaluate the identification accuracy and stability, the time-domain output error considering the condition variation is designed as the fitness function for parameter optimization. The identification process is established under various noise levels and excitation levels. The effects of external excitation and the noise level on the identification accuracy are analyzed in detail. The simulation results show that the proposed method could identify the parameters of both commensurate rate and non-commensurate rate fractional-order systems from the data with noise. It is also observed that excitation signal is an important factor influencing the identification accuracy of fractional-order systems.


Introduction
Fractional calculus [1][2][3][4] is the general expression of differential calculus.In recent years, researchers and engineers have increasingly used fractional-order dynamic models to model real physical systems that have independent frequency-domain and long memory transients [5][6][7][8][9][10][11][12][13].Some systems may have fractional-order dynamic characteristics, even if each unit has integer-order dynamic characteristics [14].What's more, applying fractional calculus to entropy theory has become a hotspot research domain [15][16][17][18][19][20]; the fractional entropy could be used in the formulation of algorithms for image segmentation where traditional Shannon entropy has presented limitations [16].In an analysis of the past ten years of trends and results in the fractional calculus application to dynamic problems of solid mechanics, the method of mechanical system dynamics analysis based on fractional calculus has gradually become one of main methods in the dynamics analysis of engineering [21].Fractional calculus has been introduced into the various engineering and science domains [22,23], including image processing [24][25][26], thermal systems identification [27,28], biological tissues identification [29][30][31], control theory and application [32][33][34][35][36], signal processing [37,38], path planning [39] and path tracking [40,41], robotics [42,43], mechanical damping [10,44], battery [45,46], mechanics [47,48], diffusion [49,50], chaos [51,52], and others.Therefore, the application of fractional calculus has become a focus of international academic research.Fractional-order system identification is a basic issue of application of fractional calculus [53][54][55][56][57][58].Several researchers have reported their work on identifying the fractional-order model in the time-domain and frequency-domain.Poinot and Trigeassou [53] proposed a time-domain method using the state-space equation, successfully obtaining the dynamical model of a heat transfer system.Cois et al. [54] modeled non-integer systems using the non-integer state-space representation, the modal coefficients, the eigenvalue, the differentiation order, and Marquardt algorithm.Lin et al. [55] used the least squares method to investigate the frequency response identification technique.Valério and Costa [56] demonstrated the fractional transfer function approximation based on phase characteristics in the frequency domain.Through the detailed analysis of these studies, it could be observed that the time domain identification proposed by Poinot and Trigeassou [53] and Cois et al. [54] could approximate most system parameters including the fractional order, but the solution of the derivation and inverse matrix is difficult and requires heavy computation.In comparison to the time-domain method, the identification methods derived by Lin et al. [55] and Valério and Costa [56] required simple calculation, but the fractional order could not be solved directly.
This paper presents an identification algorithm based on GA in the time domain.The identification process of the method is the process of parameter optimization, and the matrix inversion and differential coefficient are not needed in the method.Firstly, the effective fitness function based on output-error is put forward in the time domain, and then the multivariable parameters identification is converted into parameters optimization using GA.Secondly, the excitation signals used in parameter identification are demonstrated to indicate their effect on identification accuracy.In addition, to testify the effectiveness of the proposed methods, the identified data used in this paper adopted the real output of benchmark model coupled with various noise levels.
The organization of this paper is as follows: In Section 2, the fractional calculus, fractional-order systems and problem statement are introduced.In Section 3, the identification method based on GA is proposed.And corresponding numerical simulations and analysis are provided in Section 4. Finally, conclusions are made in Section 5.

The Definition of Fractional Calculus
There are several commonly used definitions for the general fractional differentiation and integration, such as the Grunwald-Letnikov  (GL) definition, the Riemann-Liouville (RL) definition and Caputo definition [3,21].The GL fractional derivative of continuous function f(t) is given by Grunwald  [57]: where h is the sampling period and ) denotes fractional-order differential operator.

Fractional-Order Systems
The fractional-order system is a more general expression than the integer order system; and the fractional-order system is a mathematical model based on fractional calculus.Due to the continuous order, fractional-order systems have independent frequency-domain and long memory transients [2,[5][6][7][8][9][10][11][12][13]59], which can describe complex physical system more accurately.
The single-input-single-output linear fractional-order differential equation is shown by Equation (3): In the zero initial condition, transfer function expression in the s domain of Equation ( 3) is obtained as follows: Equation ( 5) is the commensurate rate fractional-order system which is the common fractional-order system studied at this stage, and its expression is shown as follows:

Common Parameter Identification Methods
Commonly used methods of parameter estimation include least square method, maximum likelihood, correlation identification and others [60].Taking Equation (5) as an example, this article deduces an identification method based on least square method.The fractional-order differential equation of Equation ( 5) is given by: Selecting h as the sampling period, we may solve the order through using the step-by-step method, and and L+1 is total calculation number of times.Each step is , then a group of optimal coefficients are produced; we can calculate the optimal order and coefficients by using the error function J.

[ ( ), (2 ), , ( )] T y h y h y Nh
Then, the linear equation is gotten in Equation ( 12): The least square principle is introduced and its expression is shown in Equation ( 13): where e(k) is equation error and its expression is shown as follows: The residual standard is introduced by Equation ( 15): The partial derivative of J is deduced in Equation ( 16): (16) where: The estimated value of the parameters based on the least square method is obtained: The above method may figure out the coefficients and fractional order of Equation ( 5), however, Equation ( 18) could be solved only if   T is a nonsingular matrix.This might result in limitation of its application to some issues.Moreover, the above algorithm could not identify the fractional order directly, and the accuracy of identification result largely relies on step length   .Therefore, this paper introduces another identification method based on genetic algorithm in the next section.

Fractional-Order Benchmark Model
Transforming from s domain to the z domain is the discrete process of continuous transfer function, where different methods perform differently.This paper introduces first-order backwards finite difference formula [65] used as the discrete method and expands into 1,000-item truncated MacLaurin series, in order to approach the true fractional-order model.
First-order backwards finite difference formula is shown as follows: whose MacLaurin series is given in Equation (20).
where T is the sampling period, v is the fractional order, and N is power series of expansion equation.This formula is equivalent to Grunwald-Letnikov  .The classical fractional model [53,58] is taken as the benchmark model (21) for identification in Sections 4.2 and 4.3, as follows: where Continuous transfer function of the benchmark model in the s domain is described as: It is obvious that a ˆ, b ˆ, v ˆ in Equation ( 22) are parameters to be identified, and the Equation ( 22) is transformed into differential equation as follows: Then the expression of Equation (22) in the time domain is obtained as: where N is the memory length of the Equation (24).
It is assumption that i is the memory length of the Equation (24) (the maximum of i is 1,000 for this study), which is consistent with the number of simulation data.The total time t is iT.Therefore, we use y i expresses y(t), y i-k expresses y(t-kT), u i expresses u(t), the Equation ( 24) can be turned into the final expression (26):

Fitness Function of Optimization
Fitness function plays a key role in the accuracy of results.This paper takes the weighted value of output-error as the fitness function in the time domain.Therefore, parameters identification is converted into parameters optimization.
The fitness function is defined as:

Evaluation Index of Fitness Function
A series of results are obtained through the simulation and identification, the average of parameters may be regarded as net result to clear up accidental factors.Standard deviations of results can be used as the criterion of the algorithm's stability, where standard deviation is smaller usually means the algorithm is more stable.
However, because of coupling among parameters, there will be plenty of similar identification results that have the same accuracy but relevant parameter among results does not have the same value, when accuracy is not particularly high.For instance, fractional-order model (21) based on the parameters (a = 0.986, b = 0.995, v = 0.697) nearly has the same dynamic characteristics with that based on another parameters (a = 1.006, b = 0.9998, v = 0.701), although the error of each parameter (true value: a = 1, b = 1, v = 0.7) of first ones is larger than that of the latter.Therefore, it is not enough to only use each parameter's precision to evaluate the results.This paper separately takes Magnitude (dB) and Phase (degree) in the frequency domain and approximation fit in the time domain as evaluation indices to evaluate identification results.
The larger fit is, the more precise identification results are.The definition is provided as follows: where y is the real output of benchmark model without any noise.y e is the estimated output based on relevant final identification results without any noise.In order to evaluate identification results in unified standard, both y and y e are the output under the same frequency signal (VFS) excitation.

The Identification Process Based on Genetic Algorithm
GA [60][61][62][63][64] is an optimization method that models natural selection mechanism in the biological evolution process, and it has been investigated by John Holland and his students in 1975.This algorithm has global and parallel search ability and is appropriate for solving complex nonlinear problems.
In a genetic algorithm, a population of certain solutions (called individuals) to an optimization problem is evolved toward better solutions.Each potential solution has a lot of properties (its chromosomes or genotype) which can be varied and altered.Traditionally, solutions are represented in binary, but other encodings such as decimal, octal and other codes are also possible.
The evolution commonly begins with a population of randomly generated individuals and is an iterative process, with the population in each iteration called a generation.In every generation, the fitness of every individual in the population is evaluated, the more fit individuals are stochastically selected from the current population, and each individual's genome is modified (recombined and possibly randomly mutated) to form the population of next generation.This new population is then used in the next iteration of the algorithm.Usually, the algorithm would terminate when either a maximum number of generations has been produced, or other conditions have reached our request.A typical genetic algorithm requires: a genetic representation of the solution domain and an efficient fitness function to evaluate the solution domain.
Basic operation procedure of GA: Step 1: Initialization: Set the counter of evolution t = 0, the maximum number of generation T, an initial population P(0), and set other termination conditions.
Step 2: Individual evaluation: Calculate the fitness value of each individual in population P(t).
Step 3: Selection operation: Apply the selection operation to P(t) based on the fitness value.
Step 4: Crossover operation: Apply the crossover operation to P(t).
Step 5: Mutation operation: Apply the mutation operation to P(t).After Step 3, Step 4, Step 5, next generation population P(t + 1) will be gotten.
Step 6: Termination condition judgment: If t = T or meet other termination conditions, the individual which has the most suitable fitness value in the processing will be selected as the optimal solution.Otherwise, back to Step 2.
In this paper, the identification process is illustrated in Figure 1 based on the preceding analysis and basic operation process of GA.It can be seen from Figure 1 that the individual of population is composed of parameters to be identified.In the identification process, the binary encoding type and the initial population of 120 are selected.The optimization process will be stopped if one of three terminal conditions, which are 1,800 consecutive generations, 1,000 seconds runtime and 0.000001 fitness value, is satisfied.

Excitation Signals
The input signal is vital to system identification because it controls the output characteristics of the model [60,66].Taken in this sense, it also determines the accuracy of identification results and whether the system is cognizable.Different systems need different optimal excitation signals to get more built-in features.Several input signals are selected to simulate the fractional-order system, including pseudo-random binary sequence (PRBS), sawtooth wave signal, sin-swept signal, and variable frequency signal (VFS).In this paper, PRBS and VFS are a periodic square wave signals; the frequency of sawtooth wave signal is 0.05 Hz; the sin-swept signal's frequency increases from 0 to 200 Hz in one second linearly.The signals are generated by using signal functions in MATLAB.In addition, the amplitude of all the excitation signals is 1.00 and the length of data of each signal is 1,000.In engineering practice, the obtained data are contaminated by noise more or less because of the sensor precision and interference factors.To approach the actual situation, the output is joined by Gaussian white noise, and the input data have no noise.
Signal to Noise Ratio (SNR) is defined with the most commonly used method, as follows: where YW, NW express the powers of signal and noise respectively.

The Effect of Noise Level
It is well known that the noise often influences the accuracy of different identification methods.To estimate the sensitive extent of the proposed method, various noise levels of output signal is obtained by adding the Gaussian white noise, where the system input excitation is VFS without noise.The numerical simulations are carried out in the different noise level conditions.The external excitation and system response output with different SNR of model (21), and estimated output based on relevant final identification results are indicated in Figure 2 which includes output data with no noise, 28.4, 20, 16, 14 and 12 dB.
In order to reduce the stochastic error of each identification process, the mean and standard deviation of optimized parameters are adopted to evaluate the identification results.The identification result of each run is different.Therefore, in order to eliminate random error, the number of runs (five) was selected to calculate the statistical characteristics, such as the mean and standard deviation.The parameters listed in the paper are statistical values of multiple runs.For example, a, b, v are means of result of five times identifications; a  are standard deviations of parameters.At the same time, the estimation index fit is introduced to evaluate the identification results fairly under different conditions.It can be viewed from Table 1 that, in case of on noise, the identified parameters are consistent with the true values and the fit is very close to 100% and the error could be ignored.In case of 20 dB noise level, the identification result presented in Table 1 shows that the means of identified parameters a, b, v are 0.976071, 0.988563 and 0.696444, respectively.In this condition, the maximum standard deviation is 2.33917 × 10 −5 and the identification accuracy is 99.35%.While the SNR is 28.4 dB, the maximum standard deviation is 2.725825 × 10 −5 and the parameter fit could reach to 99.58%.When the SNR is 16 dB, the fit is 99.18%.Even if the SNR is equal to 14 dB, the identification accuracy keeps 98.76% and its maximum standard deviation is 0.695666 × 10 −5 .However, when the noise continues to increase, the identification accuracy will become worse.The parameter fit is only 97.09% in case of 12 dB.   2 shows maximum errors of Magnitude and Phase between benchmark model and identified models in Figure 3.It can be seen from Figure 3 that estimated models have almost the same dynamic characteristics compared to the true one.In details, in case of on noise, the maximum errors of Magnitude and Phase are 5.8848 × 10 −5 (dB) and 4.3329 × 10 −4 (degree) orderly, and the error could be ignored.When the SNR is 28.4 dB, the maximum errors of Magnitude and Phase are merely −0.0205 and 0.0793 in turn.While the SNR is between 20 dB and 14 dB, the errors are closed to each other.The maximum errors of Magnitude and Phase are 0.1131, 0.3190 in case of 20 dB, 0.1376 and 0.3404 in case of 16 dB, 0.1269 and 0.3416 in case of 14 dB.However, when the SNR is 12 dB, the errors increase apparently, and they are 2.0961 dB and 3.3482 degrees.Therefore, from the above simulation and analysis results, it is obvious that the proposed method is insensitive to noise above 14 dB.

The Effect of Excitation Signals
In the system identification, the external excitation is also an important factor to affect the identification accuracy.For the purpose of investigating the effect of excitation on identifying fractional-order system, this work designs various excitation signals, which include PRBS, VFS, sawtooth wave signal, and sin-swept signal.Meanwhile, the condition of numerical simulation is selected as system response with 28.4 dB noise level in order to keep the simulation conformable to reality.Different excitations and response outputs for the benchmark model (21), and estimated output based on relevant final identification results are indicated in Figure 2, Figure 4, Figure 5 and Figure 6.It is obvious that Figure 2 shows the system input and output under VFS excitation.Figure 4 describes the PRBS excitation.The system responses under sawtooth wave excitation and sin-swept signal excitation are exhibited in Figure 5 and Figure 6, respectively.The identification results using the proposed method are listed in Table 3.It can be seen that the higher identification accuracy is obtained using PRBS and VFS excitation.The maximum accuracy reaches to 99.58% and its corresponding estimated parameters is 1.006228, 0.99538 and 0.701493.While adopting sawtooth wave excitation, the estimated values of a, b, v is 0.996426, 0992339 and 0.711274 and the fit is 98.67%.The worse results are obtained by using sin-swept signal excitation.In this case, the identification accuracy is 81.59% and its corresponding estimated parameters are 1.999999, 1.427971 and 0.817042.Figure 7 shows the Bode diagram of benchmark model and identified models under various excitations (SNR = 28.4dB), and Table 4 shows the maximum errors of Magnitude and Phase between benchmark model and identified models in Figure 7.It can be seen from Figure 7 that estimated models under VFS and PRBS have almost the same dynamic characteristics compared to the true one.In addition, the maximum errors of Magnitude and Phase under VFS (−0.0205 dB, 0.0793 degree) are close to sawtooth wave excitation (0.0232 dB, 0.1052 degree), and are smaller than PRBS (0.1669 dB, 0.3146 degree).However, the errors are very big under sin-swept signal excitation, which are 0.9075 dB and 9.2118 degrees.Therefore, there is large difference between benchmark model and estimated model under sin-sweep excitation, which confirms the results in Table 3 and Table 4.In this paper, identification results of the benchmark model ( 21) are more accurate based on square wave excitations such as PRBS and VFS excitation than sin-swept signal excitation.Therefore, it can be viewed that the external excitation plays a significant role in parameter identification of fractional-order systems.Because optimal excitation signal might arouse the most characters of the system, different systems need different optimal excitation signals to embody more features.

Identification of General Non-Commensurate Rate Fractional-Order System
In order to verify that the proposed method is effective for general model, the general non-commensurate rate fractional-order model (model (30)) is used as benchmark model for identification, as follows: where q = 0.5, v = 0.7.
The external excitation (VFS) and system response output [model (30)] with different SNR which includes output data with 28.4 dB and 16 dB, and estimated output based on relevant final identification results are indicated in Figure 8.The Bode diagrams of benchmark model and estimated models are shown in Figure 9.When the SNR is 28.4 dB, the fit is 99.62%, and maximum errors of Magnitude and Phase are 0.0337 dB and 0.3844 degree orderly.Estimated q is 0.503864 and v is 0.703319, and relevant standard deviations are 0.237492 × 10 −5 and 0.319154 × 10 −5 .In case of 16 dB, estimated q is 0.497678 and v is 0.691340, and relevant standard deviations are 0.545049 × 10 −5 and 0.716980 × 10 −5 .The fit is 99.26%, and maximum errors of Magnitude and Phase are 0.3775 dB and 0.6656 degree.The identified results prove that the proposed method is also suitable for general fractional-order systems.

Conclusions
This paper proposes an identification algorithm based on GA in the time domain with the weighted value of output error for fractional-order systems.The results verify that this algorithm can precisely identify the coefficients and fractional-order, even when the output mixed with noise.Taking the effective fitness function, GA can do the global search and solve the parameter identification issue for fractional-order systems.In addition, it is not enough to only use each parameter's precision to evaluate the results on account of coupling among parameters.Taking errors between benchmark mode and estimated model both in the frequency domain and in the time domain as evaluation indices might be a good choice.
Excitation signal is of great importance for fractional-order systems identification; the best result might be found if the input is the optimal excitation signal of the system.However, it is very difficult to find the best signal, so it is important to select an input signal in accordance with specific conditions, which may base on antecedent analysis and experiences.What's more, the results demonstrate that this method could identify the parameters of both commensurate rate and non-commensurate rate fractional-order systems from the data with noise.Exploring the application fractional-order systems identification to engineering practice is the further work in the future.
where N is the number of data.* i y is the real output of benchmark model under different excitation signals and noise level, i y is the estimated output without any noise.as the 2-norm of an N-dimensional error vector.* i y  is the standard deviation of * i y .To eliminate the influence of the excitation signal's waveform, the error pseudo distance is divided by * i y  .

Figure 1 .
Figure 1.Flow chart of fractional-order systems identification based on genetic algorithm.

Figure 2 .
Figure 2. VFS excitation based on different noise levels of model (21).

Figure 3 .
Figure 3. Identification results under different noise levels of model (21).

Figure 3
Figure 3 shows the Bode diagram of benchmark model (Equation (21), a = 1, b = 1, v = 0.7) and identified models under various noise levels, and Table2shows maximum errors of Magnitude and Phase between benchmark model and identified models in Figure3.It can be seen from Figure3that estimated models have almost the same dynamic characteristics compared to the true one.In details, in case of on noise, the maximum errors of Magnitude and Phase are 5.8848 × 10 −5 (dB) and 4.3329 × 10 −4 (degree) orderly, and the error could be ignored.When the SNR is 28.4 dB, the

Figure 7 .
Figure 7. Identification results based on different excitations of model (21).

Table 2 .
Errors in frequency domain based on different SNR of model (21).

Table 4 .
(21)rs in frequency domain based on different excitations of model(21).