# Recursive Identification for Fractional Order Hammerstein Model Based on ADELS

Authors: Qibing Jin; Youliang Ye; Wu Cai; Zeyu Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

This paper deals with the identification of the fractional order Hammerstein model by using proposed adaptive differential evolution with the Local search strategy (ADELS) algorithm with the steepest descent method and the overparameterization based auxiliary model recursive least squares (OAMRLS) algorithm. The parameters of the static nonlinear block and the dynamic linear block of the model are all unknown, including the fractional order. The initial value of the parameter is obtained by the proposed ADELS algorithm. The main innovation of ADELS is to adaptively generate the next generation based on the fitness function value within the population through scoring rules and introduce Chebyshev mapping into the newly generated population for local search. Based on the steepest descent method, the fractional order identification using initial values is derived. The remaining parameters are derived through the OAMRLS algorithm. With the initial value obtained by ADELS, the identification result of the algorithm is more accurate. The simulation results illustrate the significance of the proposed algorithm.


---

## Body

## 1. Introduction

Currently, systems in industrial processes have become more and more increasingly complex, such as chemical plants, robotic arms, etc. This requires us to provide a more accurate mathematical model. Therefore, many parameter identification methods have been studied for system modelling and identification, such as bilinear systems, linear systems, and nonlinear systems [1–4]. Chen et al. proposed a recursive identification method to solve the problem of bilinear-parameter model identification [5]. An adaptive data filtering-based gradient iterative algorithm is presented for the system parameter identification [6]. In [7], based on Bayesian optimization, a serial two-level decomposition structure is introduced into the deep hybrid model.As a widely represented mathematical model, the nonlinear model has attracted more and more attention. The Hammerstein model is a typical nonlinear model composed of static nonlinear links and dynamic linear links. It can reflect the characteristics of process characteristics and describe a series of nonlinear processes such as neutralization processes [8,9], lithium-ion battery thermal systems [10,11], heat dissipation systems [12], and so on.Due to the wide application of the Hammerstein model, how to identify accurate mathematical models has become the research direction of many researchers. At present, some methods that can effectively identify the Hammerstein model are proposed, such as the least squares algorithm [13], maximum likelihood algorithm [14,15], etc. There are also some improved algorithms based on traditional algorithms, such as support vector machines [16,17] and multi-innovative stochastic gradient [18,19].In recent years, many scholars have researched the static nonlinear module of the Hammerstein model, hoping to obtain a mathematical model with higher accuracy and wider applicability. Some methods to describe the nonlinear part are proposed, such as radial basis functions [20], neural fuzzy networks [21,22], polynomials [23], etc. However, these studies rarely focus on dynamic linear modules, which are based on the integer order. Compared with integer-order differential operators, fractional order operators have more accurate results in the expression of heredity, memory, and other phenomena. Fractional calculus extends the order of calculus from the domain of integers to real numbers and even complex numbers. The increased degrees of freedom of fractional order operators can enable fractional order systems to describe some complex physical phenomena more accurately, such as the fuel cell physical phenomena [24], voltammetric E-tongue system [25], electric vehicle Li-ion batteries [26], etc. Experiments show that fractional order systems can more accurately model actual systems. Therefore, the research on the identification of fractional order Hammerstein model becomes very meaningful.At present, some identification methods of fractional order systems have been studied. The enhanced response sensitivity approach can reduce the sensitivity of the identification parameter results in measurement noise [27]. Block pulse functions can identify fractional order systems with initial conditions [28]. The recursive instrumental variable algorithm has acceptable accuracy for fractional order models in the absence of noise [29]. However, there are few studies about the identification problem of fractional order Hammerstein model. In [30], parameters of the fractional Hammerstein model are identified by using particle swarm optimization algorithm, which does not consider the internal relationship between system parameters. In [31], the Levenberg–Marquardt Algorithm combined with two decomposition principles is developed to identify parameters of the fractional Hammerstein model. However, when the system has noise interference, the effect of identification is not satisfactory.The intelligent optimization algorithm has attracted increasingly scholars’ attention because of its generalization, simple parameter setting, and easy programming. It has been applied to the field of parameter optimization. In [32], the modified particle swarm optimization algorithm was presented and applied to the parameter identification of fractional order PID controllers. In [33], an innovative metaheuristic algorithm called Hybrid Gray Wolf Optimization Algorithm was proposed, and it is applied to search the best parameters of the fuel cell. However, the results of the intelligent optimization algorithm may be inaccurate because of falling into local optimization.Based on the above background, this paper will adopt a new method to parameterize the nonlinear and linear coefficients and fractional order of the fractional Hammerstein model, which has rarely been considered before. The identification of each parameter has a corresponding mathematical derivation. In the proposed method, the ADELS algorithm adds a domain search strategy and improves the setting of algorithm parameters. Parameters of the fractional order Hammerstein model including the fractional order are identified. The identification result provides a relatively accurate initial value for the subsequent algorithms and solves the problem that most algorithms rely on the initial value. Then, the fractional order identification method is proposed by using the principle of steepest gradient descent, and the auxiliary model recursive least squares algorithm is used to estimate the coefficients of the fractional order Hammerstein model. Numerical simulations prove the effectiveness of the proposed methods.The main contributions are to propose an adaptive differential evolution with local search strategy (ADELS) algorithm for identifying the initial parameters of the fractional order Hammerstein system and develop new recursive identification methods for identifying the fractional order and coefficients of the fractional order nonlinear system by using the auxiliary model. Compared with the classic DE, the proposed ADELS algorithm has higher estimation accuracy because of the adaptive strategy and the introduction of Chebyshev sequence. The self-adaptive operator and mutation strategy make the individual take the corresponding search strategy according to the fitness function value, and the Chebyshev sequence has the characteristics of periodicity, randomness, and ergodicity, which makes up for the defect that DE is easy to fall into local optimum. Based on the auxiliary model, the recursive algorithm of fractional order Hammerstein model is given by using the overparameterized least squares and gradient descent algorithm. Both fractional orders and coefficients of the system are considered. The proposed ADELS algorithm can be applied to many fields like robust control, multi-objective optimization, and economic dispatch problems [34–36]. The proposed recursive identification can be extended to some other systems such as linear, bilinear, and nonlinear systems [37–39].In this paper, the mathematical background of identifying the fractional Hammerstein model is discussed in Section2. In Section 3, an improved DE algorithm is proposed. The overparameterized auxiliary model recursive least squares algorithm and fractional order identification algorithm are discussed in Section 4. In Section 5, numerical simulations prove the significance of the proposed method. Finally, Section 6 gives some conclusions.


## 2. Mathematical Background

### 2.1. Fractional Differentiation

Fractional calculus is the general form of integer calculus. So far, there are many definitions of fractional calculus. In some cases, different definitions of fractional calculus are completely equivalent. The rationality of these definitions has been proven in the literature. Among them, the definitions of three types of fractional calculus are widely used: Grünwald–Letnikov (GL), Riemann–Liouville (RL), and Caputo definitions [40].In this article, considering the ease of calculation, the GL definition will be used [41]. The α fractional derivative of ft is defined as follows:(1)Datαft=limh⟶01hα∑j=0t−a/h−1jαjft−jh,where a and t are the limits, h is the sampling time, and [.] stands for the integer part. Newton’s binomial coefficient can be calculated by using the Euler’s Gamma function as follows:(2)αj=Γα+1Γj+1Γα−j+1.By usingwjα to replace the binomial coefficient as follows:(3)wjα=−1jαj.Theα fractional derivative of ft can be approximated as [42](4)Dαft=1hα∑j=0t−a/hwjαft−jh,where Dαft denotes Datαft when a = 0.


### 2.2. Differential Evolution

Differential evolution (DE) algorithm has better stability and global search ability. The population of each generation of DE is composed of NP parameter vectors with dimension D. Each individual is represented as(5)xiG=xi,1G,xi,2G,…,xi,DG,i=1,2,…,NP,where i means the sequence of individuals in the population, G denotes the number of evolutionary generations, and NP is the size of the population.In DE, random initial populations are generated by a uniform probability distribution. After initialization, the population evolves by using three steps: mutation, crossover, and selection.


#### 2.2.1. Mutation

For each target vectorxiG at generation G, the mutation vector is generated as follows:(6)viG=xr1G+F⋅xr2G−xr3G,r1≠r2≠r3≠i∈1,NP,where r1,r2, and r3 are different individuals randomly selected from the population except the individual i and F is the mutation operator that defines the amplitude of the difference vector. When the generated mutation vector exceeds the boundary value, it will be replaced by the newly generated vector according to the boundary update rule. Some other mutation strategies have also been widely used in the DE algorithm [43–45].


#### 2.2.2. Crossover

To increase the diversity of trial vectors, crossover operations are introduced as follows:(7)uj,iG=vjiG,ifrandj≤CR orj=rk,xjiG,otherwise.CR is the crossover operator that controls the probability of accepting the mutation in each dimension, randj is a random number generated by a uniform distribution in [0,1], j represents the dimension of the individual, and rk∈1,2,…,D is a uniformly distributed random integer.


#### 2.2.3. Selection

The trial vectorui,G and the target vector xiG are compared according to the greedy criterion, and the vector with the best fitness value will appear in the next generation population. The selection operation to minimize the fitness function can be expressed as(8)xiG+1=uiG,iffuiG≤fxiGxiG,x≥0.Repeat the three steps above until the termination condition is met, such as the total number of iterations or the accuracy requirement of the algorithm.


### 2.3. Fractional Order Hammerstein Model

The fractional order Hammerstein model is shown in Figure1.Figure 1
The fractional order Hammerstein model.The system inputs and outputs areut and yt, vt is the Gaussian white noise, and y′t is the measured output. It consists of a static nonlinear part and a dynamic linear part. The fractional order differential equation of the linear part can be expressed as(9)anDnαyt+an−1Dn−1αyt+⋯+a0yt=bmDmαu′t+bm−1Dm−1αu′t+⋯+b0u′t.The fractional order transfer function of the linear part is defined as follows:(10)Gs=bmsmα+bm−1sm−1α+⋯+b0ansnα+an−1sn−1α+⋯+a0,where a1,…,an,b0,b1,…,bm are the coefficients and α is the fractional order. They are unknown and need to be identified. Without loss of generality, assume a0=1. The input u′t of the linear link is generated by the system input ut through the nonlinear block, which can be expressed by the following formula:(11)u′t=∑k=1ncck⋅fkut,where c1,…,cnc are unknown coefficients to be identified and f1⋅,…,fnc⋅ are a series of basis functions.In summary, the fractional Hammerstein system discussed in this article can be expressed as follows:(12)yt+∑i=1naaiDiαyt=∑j=0nbbjDjαu′t,y′t=yt+vt.


## 2.1. Fractional Differentiation

Fractional calculus is the general form of integer calculus. So far, there are many definitions of fractional calculus. In some cases, different definitions of fractional calculus are completely equivalent. The rationality of these definitions has been proven in the literature. Among them, the definitions of three types of fractional calculus are widely used: Grünwald–Letnikov (GL), Riemann–Liouville (RL), and Caputo definitions [40].In this article, considering the ease of calculation, the GL definition will be used [41]. The α fractional derivative of ft is defined as follows:(1)Datαft=limh⟶01hα∑j=0t−a/h−1jαjft−jh,where a and t are the limits, h is the sampling time, and [.] stands for the integer part. Newton’s binomial coefficient can be calculated by using the Euler’s Gamma function as follows:(2)αj=Γα+1Γj+1Γα−j+1.By usingwjα to replace the binomial coefficient as follows:(3)wjα=−1jαj.Theα fractional derivative of ft can be approximated as [42](4)Dαft=1hα∑j=0t−a/hwjαft−jh,where Dαft denotes Datαft when a = 0.


## 2.2. Differential Evolution

Differential evolution (DE) algorithm has better stability and global search ability. The population of each generation of DE is composed of NP parameter vectors with dimension D. Each individual is represented as(5)xiG=xi,1G,xi,2G,…,xi,DG,i=1,2,…,NP,where i means the sequence of individuals in the population, G denotes the number of evolutionary generations, and NP is the size of the population.In DE, random initial populations are generated by a uniform probability distribution. After initialization, the population evolves by using three steps: mutation, crossover, and selection.


### 2.2.1. Mutation

For each target vectorxiG at generation G, the mutation vector is generated as follows:(6)viG=xr1G+F⋅xr2G−xr3G,r1≠r2≠r3≠i∈1,NP,where r1,r2, and r3 are different individuals randomly selected from the population except the individual i and F is the mutation operator that defines the amplitude of the difference vector. When the generated mutation vector exceeds the boundary value, it will be replaced by the newly generated vector according to the boundary update rule. Some other mutation strategies have also been widely used in the DE algorithm [43–45].


### 2.2.2. Crossover

To increase the diversity of trial vectors, crossover operations are introduced as follows:(7)uj,iG=vjiG,ifrandj≤CR orj=rk,xjiG,otherwise.CR is the crossover operator that controls the probability of accepting the mutation in each dimension, randj is a random number generated by a uniform distribution in [0,1], j represents the dimension of the individual, and rk∈1,2,…,D is a uniformly distributed random integer.


### 2.2.3. Selection

The trial vectorui,G and the target vector xiG are compared according to the greedy criterion, and the vector with the best fitness value will appear in the next generation population. The selection operation to minimize the fitness function can be expressed as(8)xiG+1=uiG,iffuiG≤fxiGxiG,x≥0.Repeat the three steps above until the termination condition is met, such as the total number of iterations or the accuracy requirement of the algorithm.


## 2.2.1. Mutation

For each target vectorxiG at generation G, the mutation vector is generated as follows:(6)viG=xr1G+F⋅xr2G−xr3G,r1≠r2≠r3≠i∈1,NP,where r1,r2, and r3 are different individuals randomly selected from the population except the individual i and F is the mutation operator that defines the amplitude of the difference vector. When the generated mutation vector exceeds the boundary value, it will be replaced by the newly generated vector according to the boundary update rule. Some other mutation strategies have also been widely used in the DE algorithm [43–45].


## 2.2.2. Crossover

To increase the diversity of trial vectors, crossover operations are introduced as follows:(7)uj,iG=vjiG,ifrandj≤CR orj=rk,xjiG,otherwise.CR is the crossover operator that controls the probability of accepting the mutation in each dimension, randj is a random number generated by a uniform distribution in [0,1], j represents the dimension of the individual, and rk∈1,2,…,D is a uniformly distributed random integer.


## 2.2.3. Selection

The trial vectorui,G and the target vector xiG are compared according to the greedy criterion, and the vector with the best fitness value will appear in the next generation population. The selection operation to minimize the fitness function can be expressed as(8)xiG+1=uiG,iffuiG≤fxiGxiG,x≥0.Repeat the three steps above until the termination condition is met, such as the total number of iterations or the accuracy requirement of the algorithm.


## 2.3. Fractional Order Hammerstein Model

The fractional order Hammerstein model is shown in Figure1.Figure 1
The fractional order Hammerstein model.The system inputs and outputs areut and yt, vt is the Gaussian white noise, and y′t is the measured output. It consists of a static nonlinear part and a dynamic linear part. The fractional order differential equation of the linear part can be expressed as(9)anDnαyt+an−1Dn−1αyt+⋯+a0yt=bmDmαu′t+bm−1Dm−1αu′t+⋯+b0u′t.The fractional order transfer function of the linear part is defined as follows:(10)Gs=bmsmα+bm−1sm−1α+⋯+b0ansnα+an−1sn−1α+⋯+a0,where a1,…,an,b0,b1,…,bm are the coefficients and α is the fractional order. They are unknown and need to be identified. Without loss of generality, assume a0=1. The input u′t of the linear link is generated by the system input ut through the nonlinear block, which can be expressed by the following formula:(11)u′t=∑k=1ncck⋅fkut,where c1,…,cnc are unknown coefficients to be identified and f1⋅,…,fnc⋅ are a series of basis functions.In summary, the fractional Hammerstein system discussed in this article can be expressed as follows:(12)yt+∑i=1naaiDiαyt=∑j=0nbbjDjαu′t,y′t=yt+vt.


## 3. Adaptive Differential Evolution with Local Search Strategy

The DE algorithm is an intelligent search algorithm worked by cooperation and competition among individuals in the population. It still remains the population-based global search strategy in the evolutionary algorithm. As an efficient parallel search algorithm, the DE algorithm has strong global convergence and robustness, which is worthy of theoretical and applied research.The adaptive differential evolution with the local search strategy (ADELS) algorithm proposed in this paper has been innovatively improved in three places. Mutation and crossover operators in the DE algorithm control the enlargement ratio of the deviation vector and the probability that the trial vector comes from the mutation vector, which plays a pivotal role. Therefore, the operators in the proposed algorithm adopt the strategy of adaptive change according to the individual fitness value, which improves the optimization efficiency of the algorithm. Different mutation strategies also affect the accuracy of parameter identification. In this paper, different mutation strategies are selected according to the fitness value, which further improves the efficiency of the algorithm. Then, the local search ability of the algorithm is improved by using the Chebyshev chaotic sequence, which makes up for the shortcoming that the DE algorithm is easy to fall into the local optimum.


### 3.1. Parameters Adaption

Mutation operator and crossover operator affect the optimization ability and the convergence speed of the algorithm. Generally, these operators are fixed constants. In this paper, the individual score is obtained based on the fitness value of each individual in each generation by introducing a nonlinear scoring method, and the score ranges from 0 to 1. The operator will be determined based on the score of each individual, which makes the operator more adaptable. The scoring rules based on the minimized objective function are as follows:(13)scoreiG=0.5+0.5⋅favg−fiGfavg−fmin,fiG<favg,0.5⋅fmax−fiGfmax−favg,otherwise,where fmin,fmax,andfavg mean the minimum, maximum, and average fitness function values in the current generation population, respectively, where fiG means the value of the i individual fitness function in the G generation. It can be seen from the above formula that individuals with fitness function values below the average get higher scores. Then, the formulas of mutation operator and crossover operator can be expressed as follows:(14)FiG=Fmax−scoreiG⋅Fmax−Fmin,(15)CRiG=CRmax−scoreiG⋅CRmax−CRmin,where FiG and CRiG are the i-th mutation operator and the i-th crossover operator at generation G; Fmax and Fmin represent the limits of the mutation operator; and CRmax and CRmin are the limits of the crossover operator.It can be seen from the above formulas that the value of the operator changes with the individual’s fitness value score. When the individual’s fitness value is low, the score will be high according to equation (13), and the value of FiG and CRiG will decrease. When the individual’s fitness value is higher than the average, the score is low, and the value of FiG and CRiG will increase accordingly. The individual is close to the optimal solution when the individual’s fitness value is low, and the required mutation range and crossover probability are reduced. Then, the values of the mutation operator and crossover operator need to be reduced accordingly to find the optimal value. On the contrary, when the individual’s fitness value is high, the value of the operator needs to be increased to search for the optimal solution in a larger range.


### 3.2. Mutation Strategy

In each iteration, individuals are sorted according to the fitness value, and the population will be divided into two subpopulations with different individual fitness function. Two subpopulations adopt different mutation strategies as follows:(16)viG=xBestG+FiG⋅xr1G−xr2G,fiG<favg,viG=xr1G+FiG⋅xr2G−xr3G,otherwise,where xr1G,xr2G, and xr3G are different individuals randomly selected from the population except for the individual xiG and xBestG represents the best individual of the previous generation population. From equation (16), it is worth noticing that at least four individuals in the population are needed.


### 3.3. Local Search Strategy

To avoid the algorithm falling into local optimum, a Chebyshev chaotic sequence is introduced to perform a local search for the optimal individual. Chaotic systems have the characteristics of periodicity, randomness, and ergodicity, which can increase the diversity of individual populations. The family of Chebyshev polynomial can be expressed as(17)ρmh=cosm⋅cos−1h,m=0,1,2,….Then, the Chebyshev map can be expressed as(18)hn+1=ρmhn,hn∈−1,1,where m is the order of the polynomial. In this paper, to get better ergodicity of the mapping, the parameter of the Chebyshev map is set to m = 4 and h1=0.7.Then, a dual search strategy is introduced to search for the best gene for each dimension of the optimal individual. Due to the randomness and ergodicity of chaos, chaotic mapping has a better performance in improving the diversity of the population compared to a uniform probability distribution. The local search strategy is as follows:(19)xi,j′=xBest,j+xjmax−xjmin⋅hn+1G2,xi,j″=xBest,j+xjmax−xjmin⋅hn+1GG,where xi′ and xi″are obtained by the local search, and they will be compared with the current optimal individual based on the fitness value. If a new optimal individual is obtained, it will be randomly copied to ten individuals of the next generation population.Because of adaptive operators, different mutation strategies, and local search based on the Chebyshev chaotic map, the algorithm’s optimization ability and convergence speed have been significantly improved according to the following numerical simulations. The effectiveness of the ADELS algorithm will be proved in the numerical simulation inSection 4. The proposed ADELS Algorithm 1 is shown as follows.Algorithm 1: ADELS.
(1)
Define the objective functionfX,X=x1,x2,…,xDT;(2)
Initialize parameters of the Chebyshev map:m and h1;(3)
Initialize individualsX;(4)
Evaluate all the individuals in the population by the objective functionfX;(5)
Initialize the number of iterationk = 1;(6)
While (k < max number of iterations N)(7)
For each individual(8)
Update operators adaptively (equations (13)–(15));(9)
The mutation vectorviG is obtained by mutation (equations (16));(10)
If the generated mutation vector exceeds the boundary, a new mutation vector is generated randomly, until it is within the boundary;(11)
The trial vectorsuiG is obtained by equations (7);(12)
The best individualsX is obtained by greedy selection (equations (8));(13)
Find the current bestX according to the local search strategy (equations (17)–(19));(14)
If a new optimal individualX is obtained, it will be randomly copied to ten individuals of the next generation population;(15)
End(16)
k = k + 1;(17)
End while;(18)
Postprocess results and visualization.


## 3.1. Parameters Adaption

Mutation operator and crossover operator affect the optimization ability and the convergence speed of the algorithm. Generally, these operators are fixed constants. In this paper, the individual score is obtained based on the fitness value of each individual in each generation by introducing a nonlinear scoring method, and the score ranges from 0 to 1. The operator will be determined based on the score of each individual, which makes the operator more adaptable. The scoring rules based on the minimized objective function are as follows:(13)scoreiG=0.5+0.5⋅favg−fiGfavg−fmin,fiG<favg,0.5⋅fmax−fiGfmax−favg,otherwise,where fmin,fmax,andfavg mean the minimum, maximum, and average fitness function values in the current generation population, respectively, where fiG means the value of the i individual fitness function in the G generation. It can be seen from the above formula that individuals with fitness function values below the average get higher scores. Then, the formulas of mutation operator and crossover operator can be expressed as follows:(14)FiG=Fmax−scoreiG⋅Fmax−Fmin,(15)CRiG=CRmax−scoreiG⋅CRmax−CRmin,where FiG and CRiG are the i-th mutation operator and the i-th crossover operator at generation G; Fmax and Fmin represent the limits of the mutation operator; and CRmax and CRmin are the limits of the crossover operator.It can be seen from the above formulas that the value of the operator changes with the individual’s fitness value score. When the individual’s fitness value is low, the score will be high according to equation (13), and the value of FiG and CRiG will decrease. When the individual’s fitness value is higher than the average, the score is low, and the value of FiG and CRiG will increase accordingly. The individual is close to the optimal solution when the individual’s fitness value is low, and the required mutation range and crossover probability are reduced. Then, the values of the mutation operator and crossover operator need to be reduced accordingly to find the optimal value. On the contrary, when the individual’s fitness value is high, the value of the operator needs to be increased to search for the optimal solution in a larger range.


## 3.2. Mutation Strategy

In each iteration, individuals are sorted according to the fitness value, and the population will be divided into two subpopulations with different individual fitness function. Two subpopulations adopt different mutation strategies as follows:(16)viG=xBestG+FiG⋅xr1G−xr2G,fiG<favg,viG=xr1G+FiG⋅xr2G−xr3G,otherwise,where xr1G,xr2G, and xr3G are different individuals randomly selected from the population except for the individual xiG and xBestG represents the best individual of the previous generation population. From equation (16), it is worth noticing that at least four individuals in the population are needed.


## 3.3. Local Search Strategy

To avoid the algorithm falling into local optimum, a Chebyshev chaotic sequence is introduced to perform a local search for the optimal individual. Chaotic systems have the characteristics of periodicity, randomness, and ergodicity, which can increase the diversity of individual populations. The family of Chebyshev polynomial can be expressed as(17)ρmh=cosm⋅cos−1h,m=0,1,2,….Then, the Chebyshev map can be expressed as(18)hn+1=ρmhn,hn∈−1,1,where m is the order of the polynomial. In this paper, to get better ergodicity of the mapping, the parameter of the Chebyshev map is set to m = 4 and h1=0.7.Then, a dual search strategy is introduced to search for the best gene for each dimension of the optimal individual. Due to the randomness and ergodicity of chaos, chaotic mapping has a better performance in improving the diversity of the population compared to a uniform probability distribution. The local search strategy is as follows:(19)xi,j′=xBest,j+xjmax−xjmin⋅hn+1G2,xi,j″=xBest,j+xjmax−xjmin⋅hn+1GG,where xi′ and xi″are obtained by the local search, and they will be compared with the current optimal individual based on the fitness value. If a new optimal individual is obtained, it will be randomly copied to ten individuals of the next generation population.Because of adaptive operators, different mutation strategies, and local search based on the Chebyshev chaotic map, the algorithm’s optimization ability and convergence speed have been significantly improved according to the following numerical simulations. The effectiveness of the ADELS algorithm will be proved in the numerical simulation inSection 4. The proposed ADELS Algorithm 1 is shown as follows.Algorithm 1: ADELS.
(1)
Define the objective functionfX,X=x1,x2,…,xDT;(2)
Initialize parameters of the Chebyshev map:m and h1;(3)
Initialize individualsX;(4)
Evaluate all the individuals in the population by the objective functionfX;(5)
Initialize the number of iterationk = 1;(6)
While (k < max number of iterations N)(7)
For each individual(8)
Update operators adaptively (equations (13)–(15));(9)
The mutation vectorviG is obtained by mutation (equations (16));(10)
If the generated mutation vector exceeds the boundary, a new mutation vector is generated randomly, until it is within the boundary;(11)
The trial vectorsuiG is obtained by equations (7);(12)
The best individualsX is obtained by greedy selection (equations (8));(13)
Find the current bestX according to the local search strategy (equations (17)–(19));(14)
If a new optimal individualX is obtained, it will be randomly copied to ten individuals of the next generation population;(15)
End(16)
k = k + 1;(17)
End while;(18)
Postprocess results and visualization.


## 4. The Identification Algorithm of Fractional Order Hammerstein Model

In the fractional order Hammerstein model, the polynomial coefficients of the static nonlinear block, the coefficients of the dynamic linear block, and the fractional order are all needed to be determined, which is rarely concerned before. In this paper, a set of input and output data will be first used to obtain the initial values of coefficients and the fractional order through the proposed ADELS algorithm. Then, the initial value will be used to get the final parameter identification result of the fractional order Hammerstein model through an over-parameterized based auxiliary model recursive least squares (OAMRLS) algorithm and the steepest descent method.


### 4.1. Parameter Estimation Using Overparameterized Auxiliary Model Recursive Least Squares (OAMRLS) Algorithm

The OAMRLS algorithm is required to get the initial values of the parameters, which have been calculated in the ADELS algorithm. According to equations (1) and (2), the fractional Hammerstein system can be expressed as follows:(20)y′t=11+∑i=1naai/hiα⋅∑i=0nb∑k=1ncbihiαck⋅∑j=0t/hwjiαfkut−jh−1hiα∑j=1t/hyt−jh−∑i=1naaihiα∑j=1t/hwjiαyt−jh+vt.The input-output relations can be written in the regression form as(21)y′t=yt+vt=φTt⋅θ+vt,where yt is the noise-free output of the system and φ is the variable vector including input-output data, which can be expressed as(22)φt=φat,φb0t,φb1t,…,φbnbtT,Ymt=∑j=1t/hwjmαyt−jh,φat=−Y0t,−Y1t,…,−Ynat,Ui,pk=∑j=0t/hwjiαfput−jh,φbit=Ui,1k,Ui,2k,…,Ui,nck.The estimated vectorθ is given by(23)θ=θa,θb0,…,θbnbT,θa=a0′,a1′,…,ana′,θbj=bj′c1,bj′c2,…,bj′cnc,where(24)a0′=11+∑i=1naai/hiα,ai′=ai/hiα1+∑i=1naai/hiα,bj′=bj/hiα1+∑i=1naai/hiα,i=1,2,…,na,j=0,1,2,…,nb.It should be noted that the vectorφat contains an unknown variable yt−jh, which makes it impossible to directly use the least squares algorithm or gradient descent method to solve the problem. In view of this, the auxiliary model is used to estimate the unknown variable yt−jh. It can be seen from Figure 2 that ut is used as an input to build an auxiliary model. The main idea of the auxiliary model is that the real output of the system yt is replaced by the output of the auxiliary model yamt; then, the parameter identification problem can be solved by ut and yamt. The parameters of the auxiliary model will eventually approximate the real ones after iteration.Figure 2
The fractional order Hammerstein based on the auxiliary model.According to Figure2, Fam⋅ is the nonlinear block of the auxiliary system, and Gam⋅ represents the linear block of the auxiliary system. The expression of the auxiliary model regression form can be written as follows:(25)yamt=φamTt⋅θam,where the estimate of φ^t can be used as the value of the auxiliary model information vector φamt and the parameter estimate of θ^ can be used as the value of the auxiliary model parameter vector θam:(26)φamt=φ^t,θam=θ^.Then, the output of the auxiliary modelyamt is used to represent the estimate of φa^t:(27)φ^t=φ^at,φ^b0t,φ^b1t,…,φ^bnbtT,Ym′t=∑j=1t/hwjmαyamt−jh,φ^at=−Y0′,−Y1′,…,−Yna′,φ^bit=φbit,θ^=θ^a,θ^b0,…,θ^bnbT.Define the criterion function as(28)Jθ^=12∑i=1ty′i−φ^Tiθ^2.Then, the estimatedθ^ can be obtained by minimizing the criterion function as the following formula:(29)∂Jθ^∂θ^=−∑i=1tφ^Ti⋅y′i−φ^Tiθ^=0.Provided the inverse∑i=1tφ^i−1⋅φ^Ti−1−1 exists, then θ^t can be obtained by using least squares estimation:(30)θ^t=∑i=1tφ^i−1⋅φ^Ti−1−1⋅∑i=1tφ^Ti−1y′i.Then, a recursive version is given as the following formula:(31)θ^t=θ^t−1+Lty′t−φ^Ttθ^t−1,(32)Lt=Pt−1φ^t1+φ^TtPt−1φ^t−1,(33)Pt=In−Ltφ^TtPt−1,(34)yamt=φ^Tt⋅θ^t,where the adaptation gain matrix P is generally chosen to be started by(35)P0=1p0⋅In,0<p0≪1.For the problem of overparameterization, to ensure the uniqueness of parameter identification, the value ofc1 can be assumed to be 1 without loss of generality. Then, the value of ck and bj′ can be calculated as follows:(36)bj′=θbj1,ck=∑j=0nbθbjtbj′.Based on the identification result of the ADELS algorithm, the OAMRLS algorithm cooperates with the following order identification method to complete the identification of the system parameters.


### 4.2. Fractional Order Identification Based on Steepest Descent Method

The initial value of the fractional order has been given by the ADELS algorithm; then, the iterative optimization of the fractional order is performed according to the steepest descent method. Define the criterion function as(37)Jα=12∑t=1ky′t−y^t2.The criterion function can be minimized by the steepest descent method:(38)αn+1=αn−γn∇Jαn,where γn satisfies Jαn−γn∇Jαn<Jαn and ∇Jαn is the partial derivative of J to αn which can be expressed as follows:(39)∇Jαn=∑i=1k∂y′t−y^t∂αy′t−y^t,where(40)∂y′t−y^t∂α=−∂∂αG^sαu′t.According to formula (10), G^sα is defined as(41)G^sα=bnbsnbα+bnb−1snb−1α+⋯+b0anasnaα+ana−1sna−1α+⋯+a0,(42)−∂∂αG^sαu′t=−bnbsnbα+⋯+b0anasnaα+⋯+a02⋅na⋅anasnaα+⋯+a1sα+nb⋅bnbsnbα+⋯+b1sαanasnaα+⋯+a0⋅lns⋅u′t,where a0=1. As for the polynomial lns⋅u′t, it can be replaced by sα⋅lns/sα⋅u′t [46,47], the inverse Laplace transform of lns/sα is(43)tα−1ΓαΨα−lnt,where Ψα can be expressed as follows:(44)Ψα=1ΓαdΓαdα.Then,lns⋅u′t can be expressed as follows:(45)Dαtα−1ΓαΨα−lnt∗u′t=DαΨαDαu′t−1Γα∫0tt−τα−1lnt−τu′τdτ.Consequently,∇Jαn can be obtained by equations (39)–(45).Combined with ADELS and OAMRLS algorithms, the identification process of the fractional Hammerstein model is shown as follows.


## 4.1. Parameter Estimation Using Overparameterized Auxiliary Model Recursive Least Squares (OAMRLS) Algorithm

The OAMRLS algorithm is required to get the initial values of the parameters, which have been calculated in the ADELS algorithm. According to equations (1) and (2), the fractional Hammerstein system can be expressed as follows:(20)y′t=11+∑i=1naai/hiα⋅∑i=0nb∑k=1ncbihiαck⋅∑j=0t/hwjiαfkut−jh−1hiα∑j=1t/hyt−jh−∑i=1naaihiα∑j=1t/hwjiαyt−jh+vt.The input-output relations can be written in the regression form as(21)y′t=yt+vt=φTt⋅θ+vt,where yt is the noise-free output of the system and φ is the variable vector including input-output data, which can be expressed as(22)φt=φat,φb0t,φb1t,…,φbnbtT,Ymt=∑j=1t/hwjmαyt−jh,φat=−Y0t,−Y1t,…,−Ynat,Ui,pk=∑j=0t/hwjiαfput−jh,φbit=Ui,1k,Ui,2k,…,Ui,nck.The estimated vectorθ is given by(23)θ=θa,θb0,…,θbnbT,θa=a0′,a1′,…,ana′,θbj=bj′c1,bj′c2,…,bj′cnc,where(24)a0′=11+∑i=1naai/hiα,ai′=ai/hiα1+∑i=1naai/hiα,bj′=bj/hiα1+∑i=1naai/hiα,i=1,2,…,na,j=0,1,2,…,nb.It should be noted that the vectorφat contains an unknown variable yt−jh, which makes it impossible to directly use the least squares algorithm or gradient descent method to solve the problem. In view of this, the auxiliary model is used to estimate the unknown variable yt−jh. It can be seen from Figure 2 that ut is used as an input to build an auxiliary model. The main idea of the auxiliary model is that the real output of the system yt is replaced by the output of the auxiliary model yamt; then, the parameter identification problem can be solved by ut and yamt. The parameters of the auxiliary model will eventually approximate the real ones after iteration.Figure 2
The fractional order Hammerstein based on the auxiliary model.According to Figure2, Fam⋅ is the nonlinear block of the auxiliary system, and Gam⋅ represents the linear block of the auxiliary system. The expression of the auxiliary model regression form can be written as follows:(25)yamt=φamTt⋅θam,where the estimate of φ^t can be used as the value of the auxiliary model information vector φamt and the parameter estimate of θ^ can be used as the value of the auxiliary model parameter vector θam:(26)φamt=φ^t,θam=θ^.Then, the output of the auxiliary modelyamt is used to represent the estimate of φa^t:(27)φ^t=φ^at,φ^b0t,φ^b1t,…,φ^bnbtT,Ym′t=∑j=1t/hwjmαyamt−jh,φ^at=−Y0′,−Y1′,…,−Yna′,φ^bit=φbit,θ^=θ^a,θ^b0,…,θ^bnbT.Define the criterion function as(28)Jθ^=12∑i=1ty′i−φ^Tiθ^2.Then, the estimatedθ^ can be obtained by minimizing the criterion function as the following formula:(29)∂Jθ^∂θ^=−∑i=1tφ^Ti⋅y′i−φ^Tiθ^=0.Provided the inverse∑i=1tφ^i−1⋅φ^Ti−1−1 exists, then θ^t can be obtained by using least squares estimation:(30)θ^t=∑i=1tφ^i−1⋅φ^Ti−1−1⋅∑i=1tφ^Ti−1y′i.Then, a recursive version is given as the following formula:(31)θ^t=θ^t−1+Lty′t−φ^Ttθ^t−1,(32)Lt=Pt−1φ^t1+φ^TtPt−1φ^t−1,(33)Pt=In−Ltφ^TtPt−1,(34)yamt=φ^Tt⋅θ^t,where the adaptation gain matrix P is generally chosen to be started by(35)P0=1p0⋅In,0<p0≪1.For the problem of overparameterization, to ensure the uniqueness of parameter identification, the value ofc1 can be assumed to be 1 without loss of generality. Then, the value of ck and bj′ can be calculated as follows:(36)bj′=θbj1,ck=∑j=0nbθbjtbj′.Based on the identification result of the ADELS algorithm, the OAMRLS algorithm cooperates with the following order identification method to complete the identification of the system parameters.


## 4.2. Fractional Order Identification Based on Steepest Descent Method

The initial value of the fractional order has been given by the ADELS algorithm; then, the iterative optimization of the fractional order is performed according to the steepest descent method. Define the criterion function as(37)Jα=12∑t=1ky′t−y^t2.The criterion function can be minimized by the steepest descent method:(38)αn+1=αn−γn∇Jαn,where γn satisfies Jαn−γn∇Jαn<Jαn and ∇Jαn is the partial derivative of J to αn which can be expressed as follows:(39)∇Jαn=∑i=1k∂y′t−y^t∂αy′t−y^t,where(40)∂y′t−y^t∂α=−∂∂αG^sαu′t.According to formula (10), G^sα is defined as(41)G^sα=bnbsnbα+bnb−1snb−1α+⋯+b0anasnaα+ana−1sna−1α+⋯+a0,(42)−∂∂αG^sαu′t=−bnbsnbα+⋯+b0anasnaα+⋯+a02⋅na⋅anasnaα+⋯+a1sα+nb⋅bnbsnbα+⋯+b1sαanasnaα+⋯+a0⋅lns⋅u′t,where a0=1. As for the polynomial lns⋅u′t, it can be replaced by sα⋅lns/sα⋅u′t [46,47], the inverse Laplace transform of lns/sα is(43)tα−1ΓαΨα−lnt,where Ψα can be expressed as follows:(44)Ψα=1ΓαdΓαdα.Then,lns⋅u′t can be expressed as follows:(45)Dαtα−1ΓαΨα−lnt∗u′t=DαΨαDαu′t−1Γα∫0tt−τα−1lnt−τu′τdτ.Consequently,∇Jαn can be obtained by equations (39)–(45).Combined with ADELS and OAMRLS algorithms, the identification process of the fractional Hammerstein model is shown as follows.


## 5. Experimental Results

### 5.1. Benchmark Functions

Firstly, 15 benchmark functions will be used to verify the effectiveness of the ADELS algorithm, and the results will be compared with other intelligent optimization algorithms proposed in recent years. These functions include unimodal and multimodal functions. The function expression and search range are shown in Table1.Table 1
Benchmark functions and search range.
FunctionRangeF1x=∑i=1Dxi2[−100,100]F2x=∑i=1Dxi+∏i=1Dxi[−10,10]F3x=∑i=1Dxi+0.52[−100,100]F4x=∑i=1D−xi⋅sinxi[−500,500]F5x=∑i=1Di⋅xi2[−10,10]F6x=x12+106∑i=2Dxi2[−100,100]F7x=∑i=1D−1100⋅xi2−xi+12+xi−12[−100,100]F8x=sin2πw1+wD−121+sin22πwD+∑i=1D−1wi−12⋅1+10sin2πwi+1,wi=1+xi−1/4,i=1,…,D[−100,100]F9x=∑i=1D106i−1/D−1⋅xi2[−100,100]F10x=106⋅x12+∑i=2Dxi2[−100,100]F11x=−20exp−0.21/D∑i=1Dxi2−exp1/D∑i=1Dcos2πxi+20+e[−32,32]F12x=∑i=1D∑k=020akcos2πbkxi+0.5−D∑k=020akcos2πbk⋅0.5,a=0.5,b=3[−100,100]F13x=∑i=1Dxi2/4000−∏i=1Dcosxi/i+1[−100,100]F14x=10/D2∏iD1+i∑j=1322jxi−round2jxi/2j10/D1.2−10/D2[−100,100]F15x=∑i=1Dxi2−10cos2πxi+10[−5.12,5.12]In this paper, the ADELS algorithm will be compared with three algorithms including JADE [48], GWO [49], and WOA [50]. Their parameter settings are given in Table 2.Table 2
Parameter settings.
AlgorithmParameter settingsJADECR = 0.5,F = 0.5, size = 100, P0=0.05GWOSize = 50WOASize = 50ADELSSize = 50,Fmin=0.1,Fmax=0.9,CRmin=0.1,CRmax=0.9In this section, the population size of the four algorithms is set to 50, the number of iterations is 500, and the simulation dimensions of the benchmark function are 10 and 30. The simulation results are shown in Tables3 and 4.Table 3
Comparison of optimization results obtained for the 30-dimensional benchmark functions.
FJADEGWOWOAADELSAveStdAveStdAveStdAveStdF12.55E–222.43E–225.46E–334.78E–333.35E–862.52E–8300F27.65E–115.09E–117.26E–206.93E–202.13E–545.88E–521.80E–1900F32.01E–221.98E–220.45870.200.0866.7E–0200F4−11684.98142.47−6281.5861.64−11535.271685.05−12391.8376.85F51.96E–214.00E–213.46E–327.43E–323.12E–852.68E–8400F61.01E–159.11E–161.14E–272.34E–273.15E–787.32E–7900F717.692.0626.670.8727.540.330.239119.20F82.3E–020.100.940.220.490.281.50E–325.62E–48F91.78E–176.95E–183.52E–301.63E–291.79E–808.55E–8200F105.25E–215.43E–217.65E–334.43E–332.44E–831.12E–8000F113.05E–124.89E–124.42E–144.53E–153.73E–152.55E–151.15E–153.63E–15F12−675.4910.04−476.5166.40−870.001.02E–12−870.003.61E–13F139.60E–104.28E–091.4E–037.8E–035.1E–030.007800F149.11E–126.81E–133.68E–115.30E–120000F1526.513.412.693.380000Table 4
Comparison of optimization results obtained for the 10-dimensional benchmark functions.
FJADEGWOWOAADELSAveStdAveStdAveStdAveStdF19.16E–863.47E–858.14E–703.06E–698.55E–892.07E–8800F26.74E–562.22E–552.88E–402.65E–402.43E–579.34E–575.77E–1910F31.46E–041.78E–042.56E–067.58E–071.99E–042.2E–0400F4−3542.12672.55−2806.47322.09−3439.87567.21−4189.8360.45F52.74E–881.22E–871.01E–684.31E–681.34E–894.04E–8900F61.42E–806.29E–801.94E–644.90E–641.79E–845.36E–8400F76.676.38E–016.535.31E–016.604.79E–011.99E–032.61E–03F86.81E–021.01E–011.02E–017.3E–024.77E–028.06E–021.50E–325.62E–48F95.18E–852.15E–842.08E–674.60E–671.58E–843.50E–8400F106.55E–872.32E–867.91E–691.83E–686.55E–872.91E–8600F113.02E–152.13E–156.75E–151.74E–153.73E–152.47E–154.80E–141.09E–14F12−90.000−82.637.16−90.000−90.000F136.80E–021.02E–011.8E–022.7E–025.59E–028.67E–023.11E–022.67E–02F14002.16E–111.85E–110000F15005.99E–011.850000The convergence curves of ADELS, GWO, WOA, and JADE are as shown in Figures3 and 4; ADELS is very competitive with other metaheuristic algorithms due to its adaptive operators and local search with the Chebyshev chaotic map.Figure 3
Comparison of convergence curves of WOA and literature algorithms obtained in some of the 30-dimensional benchmark problems.Figure 4
Comparison of convergence curves of WOA and literature algorithms obtained in some of the 10-dimensional benchmark problems.The unimodal function can be used to evaluate the exploitation capability of the optimization algorithm, whereas the multimodal function is very effective to evaluate the exploration capability of the algorithm. The results reported in Tables3 and 4 indicate that the exploitation and exploration capabilities of ADELS are significantly improved compared to other algorithms.


### 5.2. Identification of Fractional Hammerstein Model

In this section, several different algorithms will be compared to prove the effectiveness of the proposed algorithm. The fractional Hammerstein model is presented as follows:(46)a2D2αyt+a1Dαyt+yt=b1Dαu′t+b0u′t,u′t=ut+c2u2t+c3u3t+c4u4t,y′t=yt+vt.The parameter vector to be identified is(47)θ=a1,a2,b0,b1,c2,c3,c4,αT=1.3,1.5,1.9,1.7,0.5,0.3,0.1,0.3T.The inputu is a persistent excitation signal sequence with zero mean and unit variance. v(t) is the stochastic Gaussian noise with zero mean and σ2=0.001. Then, the output y(t) is generated by the corresponding fractional Hammerstein model.The ADELS algorithm is used to estimate initial parameter values of the nonlinear and linear parts of the system. The population size of the algorithm is set to 40, the minimum and maximum values of the mutation operator are 0.1 and 0.9, the minimum and maximum values of the crossover operator are 0.1 and 0.9, and the number of iterations of the algorithm is set to 200. The parameters estimated by the optimization algorithm converge to the true value, and the relative quadratic error is as follows:(48)RQE=θ^−θ2θ2,where θ is the true parameter vector andθ^ is the estimated parameter vector. The simulation results compared with other approaches are shown in Table 5. It reveals that more accurate initial values can be obtained by ADELS than other heuristic algorithms.Table 5
Initial values obtained by three different heuristic algorithms.
ApproachParameterRqea1a2b0b1c2c3c4α1GWO1.08631.77691.82391.86150.51590.29470.09320.28410.14592JADE0.38281.90531.63461.59200.53990.27280.08170.26030.38753ADELS1.17301.53981.85161.72890.50650.28820.09500.29700.0539Simulation results obtained by different algorithms are given in Table6. In the case of the same initial value obtained by ADELS, the accuracy of OAMRLS compared to KTRLS has been improved effectively.Table 6
Simulation results obtained by different algorithms.
ApproachParameterRqea1a2b0b1c2c3c4α1GWO & OAMRLS1.32331.48211.89951.71730.50390.29120.09530.30220.01332JADE & OAMRLS1.36991.45821.90441.73390.50650.29140.09460.30630.03323ADELS & KTSG1.19321.56741.88671.72890.50650.28820.09500.29690.04864ADELS & OAMRLS1.30601.48991.89771.71000.50290.29120.09560.30070.0069The performance of the intelligent optimization algorithm to search for the initial value can be measured by Root Mean Square Error (RMSE) and Mean Square Error (MSE), which are defined as follows:(49)RMSE=MSE,MSE=∑ni=1yi−y^i2n,where yi is the system output and y^i is the estimated system output. The output estimation errors of different intelligent optimization algorithms are depicted in Figure 5.Figure 5
The output estimation errors of different intelligent optimization algorithms.Using the initial estimation value, OAMRLS and the steepest descent method are combined to estimate the accurate parameter vector. Figure6 shows the actual outputs and the identified outputs regarding Algorithm 2, and the output error of the system identified by Algorithm 2 is shown in Figure 7. It can be seen from Figure 8 that the identified model can stably track the real output of the system under the step response. Figure 9 illustrates the trend of estimated fractional order convergence curve.Figure 6
The actual outputs and the identified outputs regarding Algorithm2.Algorithm 2: Identification process of the fractional Hammerstein model.
(1)
Collect the input/output data {ut,yt};(2)
Obtain the initial of unknown parameters by using Algorithm1;(3)
While (t < max number of iterations N)(4)
Estimate the value of system coefficient according to the equation (31);(5)
Update the fractional order according to the equation (38);(6)
If the criterion function valueJ satisfies the error accuracy(7)
Break;(8)
End;(9)
t = t + 1;(10)
End while;(11)
Postprocess results and visualization.Figure 7
The error between the estimated and the actual outputs regarding Algorithm2.Figure 8
Step response of real system and estimation model regarding Algorithm2.Figure 9
The estimated fractional order regarding Algorithm2.


## 5.1. Benchmark Functions

Firstly, 15 benchmark functions will be used to verify the effectiveness of the ADELS algorithm, and the results will be compared with other intelligent optimization algorithms proposed in recent years. These functions include unimodal and multimodal functions. The function expression and search range are shown in Table1.Table 1
Benchmark functions and search range.
FunctionRangeF1x=∑i=1Dxi2[−100,100]F2x=∑i=1Dxi+∏i=1Dxi[−10,10]F3x=∑i=1Dxi+0.52[−100,100]F4x=∑i=1D−xi⋅sinxi[−500,500]F5x=∑i=1Di⋅xi2[−10,10]F6x=x12+106∑i=2Dxi2[−100,100]F7x=∑i=1D−1100⋅xi2−xi+12+xi−12[−100,100]F8x=sin2πw1+wD−121+sin22πwD+∑i=1D−1wi−12⋅1+10sin2πwi+1,wi=1+xi−1/4,i=1,…,D[−100,100]F9x=∑i=1D106i−1/D−1⋅xi2[−100,100]F10x=106⋅x12+∑i=2Dxi2[−100,100]F11x=−20exp−0.21/D∑i=1Dxi2−exp1/D∑i=1Dcos2πxi+20+e[−32,32]F12x=∑i=1D∑k=020akcos2πbkxi+0.5−D∑k=020akcos2πbk⋅0.5,a=0.5,b=3[−100,100]F13x=∑i=1Dxi2/4000−∏i=1Dcosxi/i+1[−100,100]F14x=10/D2∏iD1+i∑j=1322jxi−round2jxi/2j10/D1.2−10/D2[−100,100]F15x=∑i=1Dxi2−10cos2πxi+10[−5.12,5.12]In this paper, the ADELS algorithm will be compared with three algorithms including JADE [48], GWO [49], and WOA [50]. Their parameter settings are given in Table 2.Table 2
Parameter settings.
AlgorithmParameter settingsJADECR = 0.5,F = 0.5, size = 100, P0=0.05GWOSize = 50WOASize = 50ADELSSize = 50,Fmin=0.1,Fmax=0.9,CRmin=0.1,CRmax=0.9In this section, the population size of the four algorithms is set to 50, the number of iterations is 500, and the simulation dimensions of the benchmark function are 10 and 30. The simulation results are shown in Tables3 and 4.Table 3
Comparison of optimization results obtained for the 30-dimensional benchmark functions.
FJADEGWOWOAADELSAveStdAveStdAveStdAveStdF12.55E–222.43E–225.46E–334.78E–333.35E–862.52E–8300F27.65E–115.09E–117.26E–206.93E–202.13E–545.88E–521.80E–1900F32.01E–221.98E–220.45870.200.0866.7E–0200F4−11684.98142.47−6281.5861.64−11535.271685.05−12391.8376.85F51.96E–214.00E–213.46E–327.43E–323.12E–852.68E–8400F61.01E–159.11E–161.14E–272.34E–273.15E–787.32E–7900F717.692.0626.670.8727.540.330.239119.20F82.3E–020.100.940.220.490.281.50E–325.62E–48F91.78E–176.95E–183.52E–301.63E–291.79E–808.55E–8200F105.25E–215.43E–217.65E–334.43E–332.44E–831.12E–8000F113.05E–124.89E–124.42E–144.53E–153.73E–152.55E–151.15E–153.63E–15F12−675.4910.04−476.5166.40−870.001.02E–12−870.003.61E–13F139.60E–104.28E–091.4E–037.8E–035.1E–030.007800F149.11E–126.81E–133.68E–115.30E–120000F1526.513.412.693.380000Table 4
Comparison of optimization results obtained for the 10-dimensional benchmark functions.
FJADEGWOWOAADELSAveStdAveStdAveStdAveStdF19.16E–863.47E–858.14E–703.06E–698.55E–892.07E–8800F26.74E–562.22E–552.88E–402.65E–402.43E–579.34E–575.77E–1910F31.46E–041.78E–042.56E–067.58E–071.99E–042.2E–0400F4−3542.12672.55−2806.47322.09−3439.87567.21−4189.8360.45F52.74E–881.22E–871.01E–684.31E–681.34E–894.04E–8900F61.42E–806.29E–801.94E–644.90E–641.79E–845.36E–8400F76.676.38E–016.535.31E–016.604.79E–011.99E–032.61E–03F86.81E–021.01E–011.02E–017.3E–024.77E–028.06E–021.50E–325.62E–48F95.18E–852.15E–842.08E–674.60E–671.58E–843.50E–8400F106.55E–872.32E–867.91E–691.83E–686.55E–872.91E–8600F113.02E–152.13E–156.75E–151.74E–153.73E–152.47E–154.80E–141.09E–14F12−90.000−82.637.16−90.000−90.000F136.80E–021.02E–011.8E–022.7E–025.59E–028.67E–023.11E–022.67E–02F14002.16E–111.85E–110000F15005.99E–011.850000The convergence curves of ADELS, GWO, WOA, and JADE are as shown in Figures3 and 4; ADELS is very competitive with other metaheuristic algorithms due to its adaptive operators and local search with the Chebyshev chaotic map.Figure 3
Comparison of convergence curves of WOA and literature algorithms obtained in some of the 30-dimensional benchmark problems.Figure 4
Comparison of convergence curves of WOA and literature algorithms obtained in some of the 10-dimensional benchmark problems.The unimodal function can be used to evaluate the exploitation capability of the optimization algorithm, whereas the multimodal function is very effective to evaluate the exploration capability of the algorithm. The results reported in Tables3 and 4 indicate that the exploitation and exploration capabilities of ADELS are significantly improved compared to other algorithms.


## 5.2. Identification of Fractional Hammerstein Model

In this section, several different algorithms will be compared to prove the effectiveness of the proposed algorithm. The fractional Hammerstein model is presented as follows:(46)a2D2αyt+a1Dαyt+yt=b1Dαu′t+b0u′t,u′t=ut+c2u2t+c3u3t+c4u4t,y′t=yt+vt.The parameter vector to be identified is(47)θ=a1,a2,b0,b1,c2,c3,c4,αT=1.3,1.5,1.9,1.7,0.5,0.3,0.1,0.3T.The inputu is a persistent excitation signal sequence with zero mean and unit variance. v(t) is the stochastic Gaussian noise with zero mean and σ2=0.001. Then, the output y(t) is generated by the corresponding fractional Hammerstein model.The ADELS algorithm is used to estimate initial parameter values of the nonlinear and linear parts of the system. The population size of the algorithm is set to 40, the minimum and maximum values of the mutation operator are 0.1 and 0.9, the minimum and maximum values of the crossover operator are 0.1 and 0.9, and the number of iterations of the algorithm is set to 200. The parameters estimated by the optimization algorithm converge to the true value, and the relative quadratic error is as follows:(48)RQE=θ^−θ2θ2,where θ is the true parameter vector andθ^ is the estimated parameter vector. The simulation results compared with other approaches are shown in Table 5. It reveals that more accurate initial values can be obtained by ADELS than other heuristic algorithms.Table 5
Initial values obtained by three different heuristic algorithms.
ApproachParameterRqea1a2b0b1c2c3c4α1GWO1.08631.77691.82391.86150.51590.29470.09320.28410.14592JADE0.38281.90531.63461.59200.53990.27280.08170.26030.38753ADELS1.17301.53981.85161.72890.50650.28820.09500.29700.0539Simulation results obtained by different algorithms are given in Table6. In the case of the same initial value obtained by ADELS, the accuracy of OAMRLS compared to KTRLS has been improved effectively.Table 6
Simulation results obtained by different algorithms.
ApproachParameterRqea1a2b0b1c2c3c4α1GWO & OAMRLS1.32331.48211.89951.71730.50390.29120.09530.30220.01332JADE & OAMRLS1.36991.45821.90441.73390.50650.29140.09460.30630.03323ADELS & KTSG1.19321.56741.88671.72890.50650.28820.09500.29690.04864ADELS & OAMRLS1.30601.48991.89771.71000.50290.29120.09560.30070.0069The performance of the intelligent optimization algorithm to search for the initial value can be measured by Root Mean Square Error (RMSE) and Mean Square Error (MSE), which are defined as follows:(49)RMSE=MSE,MSE=∑ni=1yi−y^i2n,where yi is the system output and y^i is the estimated system output. The output estimation errors of different intelligent optimization algorithms are depicted in Figure 5.Figure 5
The output estimation errors of different intelligent optimization algorithms.Using the initial estimation value, OAMRLS and the steepest descent method are combined to estimate the accurate parameter vector. Figure6 shows the actual outputs and the identified outputs regarding Algorithm 2, and the output error of the system identified by Algorithm 2 is shown in Figure 7. It can be seen from Figure 8 that the identified model can stably track the real output of the system under the step response. Figure 9 illustrates the trend of estimated fractional order convergence curve.Figure 6
The actual outputs and the identified outputs regarding Algorithm2.Algorithm 2: Identification process of the fractional Hammerstein model.
(1)
Collect the input/output data {ut,yt};(2)
Obtain the initial of unknown parameters by using Algorithm1;(3)
While (t < max number of iterations N)(4)
Estimate the value of system coefficient according to the equation (31);(5)
Update the fractional order according to the equation (38);(6)
If the criterion function valueJ satisfies the error accuracy(7)
Break;(8)
End;(9)
t = t + 1;(10)
End while;(11)
Postprocess results and visualization.Figure 7
The error between the estimated and the actual outputs regarding Algorithm2.Figure 8
Step response of real system and estimation model regarding Algorithm2.Figure 9
The estimated fractional order regarding Algorithm2.


## 6. Summary

This paper outlines the identification methods for fractional order Hammerstein models. To improve the accuracy of identification, the heuristic algorithm is considered to search the initial value. Simulation results show that the heuristic algorithm is easy to fall into local optimum. To solve this problem, ADELS algorithm is proposed. In this algorithm, combined with Chebyshev mapping, genes are adaptively searched in each generation, which leads to faster convergence and more accurate results. The initial values of the coefficients of the linear part and the nonlinear part are obtained by the ADELS. Furthermore, the fractional order of the transfer function obtains the initial value through the ADELS. Then, the method of fractional order Hammerstein model coefficient estimation is proposed and compared with other algorithms. When the initial value is obtained, the parameters can be accurately identified by the corresponding OAMRLS algorithm. The fractional order is iterated by the steepest descent algorithm through derivation. Through comparison with other algorithms, the simulation results show the effectiveness of the proposed algorithm. The proposed methods in this paper can be applied to other literatures [51–56] such as parameter identification problems of different systems, engineering applications, fault diagnosis, and so on.



---

*Source: 6629820-2021-02-23.xml*