Rate-Dependent Modeling of Piezoelectric Actuators for Nano Manipulation Based on Fractional Hammerstein Model

Piezoelectric actuators (PEAs), as a smart material with excellent characteristics, are increasingly used in high-precision and high-speed nano-positioning systems. Different from the usual positioning control or fixed frequency tracking control, the more accurate rate-dependent PEA nonlinear model is needed in random signal dynamic tracking control systems such as active vibration control. In response to this problem, this paper proposes a Hammerstein model based on fractional order rate correlation. The improved Bouc-Wen model is used to describe the asymmetric hysteresis characteristics of PEA, and the fractional order model is used to describe the dynamic characteristics of PEA. The nonlinear rate-dependent hysteresis model can be used to accurately describe the dynamic characteristics of PEA. Compared with the integer order model or linear autoregressive model to describe the dynamic characteristics of the PEA Hammerstein model, the modeling accuracy is higher. Moreover, an artificial bee colony algorithm (DE-ABC) based on differential evolution was proposed to identify model parameters. By adding the mutation strategy and chaos search of the genetic algorithm into the previous ABC, the convergence speed of the algorithm is faster and the identification accuracy is higher, and the simultaneous identification of order and coefficient of the fractional model is realized. Finally, by comparing the simulation and experimental data of multiple sets of sinusoidal excitation with different frequencies, the effectiveness of the proposed modeling method and the accuracy and rapidity of the identification algorithm are verified. The results show that, in the wide frequency range of 1–100 Hz, the proposed method can obtain more accurate rate-correlation models than the Bouc-Wen model, the Hammerstein model based on integer order or the linear autoregressive model to describe dynamic characteristics. The maximum error (Max error) is 0.0915 μm, and the maximum mean square error (RMSE) is 0.0244.


Introduction
PEA is a smart material with excellent performance. It has outstanding advantages such as large force, high rigidity, high control accuracy, low power consumption, and fast response speed. Therefore, it has been widely used, such as in micro-manipulation [1], a micro-mechanical arm [2], active optical components [3], active vibration control [4], biomedical engineering [5], etc. In the application of PEA for active vibration control, it is mainly aimed at the low-frequency vibration with a frequency below 100 Hz, at which time the traditional passive vibration isolation method is difficult to work [6]. PEA is adopted to produce motion with the same size and opposite direction as the vibration wave, thus weakening the vibration amplitude. In this process, PEA is required to accurately track the reference signal with random and continuous changes in amplitude and frequency. Therefore, the modeling accuracy of its nonlinear characteristics such as hysteresis and Bouc-Wen model [25]. Although the single identification algorithm mentioned above can identify the parameters of the model, its premature phenomenon, which fast convergence to the local optimal solution rather than the global optimal solution, swing near the optimal solution when approaching the optimal solution, and slow convergence speed, will affect the identification accuracy of the model, and then the accurate model cannot be obtained.
In response to the above problems, the main research content and contributions of this article can be summarized as follows: 1.
In terms of PEA rate-dependent modeling, this paper proposes a modeling method based on the improved Bouc-Wen model and fractional-order model consisting of a separate Hammerstein model, which achieves a more accurate description of the dynamic characteristics of PEA rate-dependent hysteresis. At present, most Hammerstein models adopt a combination of traditional hysteresis models and integer-order or linear autoregressive models to model PEA. When using integer-order or linear autoregressive models to describe the dynamic characteristics of PEA, the accuracy is far inferior to that of fractional-order models. Compared with the integer-order or linear autoregressive models, it is more in line with engineering reality and can contain richer amplitude-frequency information.

2.
In terms of model parameter identification, this paper proposes an artificial bee colony algorithm based on differential evolution (DE-ABC), and uses it for the first time for the parameter identification of the fractional Hammerstein model. By adding the mutation strategy and chaos search for genetic algorithm into the previous ABC, the convergence speed of the algorithm is faster and the identification accuracy is higher, and the simultaneous identification of order and coefficient of fractional model is realized.

3.
Through experimental data collection and analysis, the proposed model has a good frequency generalization ability within 1-100 Hz, and can better reflect the true characteristics of piezoelectric ceramic actuators than the traditional static hysteresis model.
The paper's structure is arranged as follows. Section 2 depicts the model structure, consisting of an improved Bouc-Wen model, a fractional-order model, and a rate-dependent Hammerstein model. Section 3 proposes to use the artificial bee colony algorithm based on differential evolution for the identification of the improved Bouc-Wen static hysteresis nonlinear model and the fractional dynamic model. Finally, Section 4 introduces the experimental design, compares the identification results with different algorithms, and validates the Hammerstein model and piezoelectric actuator. The staple conclusions are summed up in Section 5.

Fractional Hammerstein Model of PEA
The rate-dependent hysteresis system refers to a system whose output is related to the current and previous inputted signals and their frequency. As shown in Figure 1, PEA have rate-dependent hysteresis nonlinearity, that is, when the piezoelectric actuator changes greatly in the frequency of the input voltage, the relationship between its input and output changes greatly.
At present, there is no unified mathematical model for nonlinear systems. As a nonlinear model based on module connection, the Hammerstein model can more precisely depict the rate-dependent characteristics of PEA. The classic Hammerstein model [24,26] contains static nonlinearity and linear dynamic model [27,28]. Its structure diagram is shown in Figure 2 At present, there is no unified mathematical model for nonlinear systems. As a nonlinear model based on module connection, the Hammerstein model can more precisely depict the rate-dependent characteristics of PEA. The classic Hammerstein model [24,26] contains static nonlinearity and linear dynamic model [27,28]. Its structure diagram is shown in Figure 2.

Static nonlinearity N(· )
Linear dynamics G(s) Experiments show that when the frequency of the drive signal is low, the hysteresis curve of the PEA shows rate-independence; when the frequency is high, the hysteresis curve shows rate-dependence. So as to bewrite the rate-dependent hysteresis nonlinearity of PEA, this paper proposes a rate-dependent hysteresis nonlinear model based on the Bouc-Wen model and fractional-order model. Among them, the linear dynamic part of the Hammerstein model is described by the fractional-order model, and the static nonlinear part is described by the Bouc-Wen model. The structure is shown in Figure 3.   At present, there is no unified mathematical model for nonlinear systems. As a nonlinear model based on module connection, the Hammerstein model can more precisely depict the rate-dependent characteristics of PEA. The classic Hammerstein model [24,26] contains static nonlinearity and linear dynamic model [27,28]. Its structure diagram is shown in Figure 2.

Static nonlinearity N(· )
Linear dynamics G(s) Experiments show that when the frequency of the drive signal is low, the hysteresis curve of the PEA shows rate-independence; when the frequency is high, the hysteresis curve shows rate-dependence. So as to bewrite the rate-dependent hysteresis nonlinearity of PEA, this paper proposes a rate-dependent hysteresis nonlinear model based on the Bouc-Wen model and fractional-order model. Among them, the linear dynamic part of the Hammerstein model is described by the fractional-order model, and the static nonlinear part is described by the Bouc-Wen model. The structure is shown in Figure 3.  Experiments show that when the frequency of the drive signal is low, the hysteresis curve of the PEA shows rate-independence; when the frequency is high, the hysteresis curve shows rate-dependence. So as to bewrite the rate-dependent hysteresis nonlinearity of PEA, this paper proposes a rate-dependent hysteresis nonlinear model based on the Bouc-Wen model and fractional-order model. Among them, the linear dynamic part of the Hammerstein model is described by the fractional-order model, and the static nonlinear part is described by the Bouc-Wen model. The structure is shown in Figure 3. At present, there is no unified mathematical model for nonlinear systems. As a nonlinear model based on module connection, the Hammerstein model can more precisely depict the rate-dependent characteristics of PEA. The classic Hammerstein model [24,26] contains static nonlinearity and linear dynamic model [27,28]. Its structure diagram is shown in Figure 2.

Static nonlinearity N(· )
Linear dynamics G(s) Experiments show that when the frequency of the drive signal is low, the hysteresis curve of the PEA shows rate-independence; when the frequency is high, the hysteresis curve shows rate-dependence. So as to bewrite the rate-dependent hysteresis nonlinearity of PEA, this paper proposes a rate-dependent hysteresis nonlinear model based on the Bouc-Wen model and fractional-order model. Among them, the linear dynamic part of the Hammerstein model is described by the fractional-order model, and the static nonlinear part is described by the Bouc-Wen model. The structure is shown in Figure 3.

The Improved Bouc-Wen Model
The classic Bouc-Wen model was introduced by Bouc [29], and Wen [30] extended it. The Bouc-Wen model contains fewer parameters and can be recognized by less data than other models depicting hysteresis characteristics. But in reality, the PEA input has an inherently multi-value memory dependent hysteresis loop asymmetry. Therefore, the improved Bouc-Wen model is adopted in this paper. By defining l[u(t)] as a generalized input function, that is, the function of input voltage u(t), it is used to describe the asymmetric property of PEA input and output about origin asymmetry. The improved Bouc-Wen model expression is: In the formula, x(t) is the output displacement; l[u(t)] is a nonsingular input function; z(t) is the hysteresis displacement component; . z(t) is the derivative of z(t) with respect to time; u(t) represents the input voltage; . u(t) is the derivative of u(t) with respect to time; x and y control the asymmetry of input signals; α controls the size of the hysteresis loop; β and γ respectively control the shape of the hysteresis loop; n represents the smoothness of the transition from the elastic part to the sculpted part. By appropriately selecting model parameters, it can express various hysteresis in various shapes.
However, The improved Bouc-Wen model still does not have the rate-independent features. Considering these limitations, modelling errors cannot be ignored in practical applications [31].

Fractional Dynamic Model
Most Hammerstein models use integer-order calculus models or linear autoregressive model (ARX) models to describe the dynamic characteristics of PEA, which often overlook some real phenomena and properties with fractional-order characteristics. While using the ARX model to describe dynamic characteristics, the output at the current moment is not only determined by the input at the current moment, but also by the input and output at all previous moments. The fractional-order model can be a good substitute for the ARX model due to its long memory characteristics. When using an integer-order model to describe the dynamic characteristics, it is necessary to analyze the dynamic characteristics of the PEA and to estimate the redundant parameters that affect the robustness of the control application. However, the fractional-order system is considered to have improved robustness in the control design. Moreover, due to the mathematical characteristics of fractional globality and long memory, it not only has the advantages of the ARX model but also can use fewer parameters to model complex systems. It has been widely used in physics and control [32,33]. The fractional Hammerstein model avoids complicated internal mechanism analysis on the basis of obtaining good modelling accuracy.
Fractional calculus is the generalization of integer order integration and differentiation to all real numbers. The basic operator t0 D α t is defined as [34]: Among them, t 0 is the lower limit and t is the upper limit of the integral, α is the fractional-order, which can be a complex number, and Re(α) is the real part of α.
There are many definitions of fractional calculus, among which the most commonly used are the Grunwald-Letnikov (G-L) definition. This paper uses the approximate calculation defined by G-L to carry out the numerical simulation of fractional operators. The definition of G-L is given as [34]: Fractional systems can be represented by fractional linear differential equations, which have the form: Among them, u(t) and x(t) are system input and system output.
If the initial condition is zero, the Laplace transform of D α x(t) is: Take the Laplace to transform on both ends of the above Equation (4). The fractionalorder linear time-invariant system can be rewritten as the below transfer function form: where α 0 < α 1 < · · · < α n , and β 0 < β 1 < · · · < β m .

Artificial Bee Colony Algorithm Based on Differential Evolution (DE-ABC)
Due to the complexity of the fractional systems, the dynamic model using fractionalorder description is not easy to estimate. In addition, due to the existence of multiple variables in the problem, there are multiple local search optimal solutions in the objective function, which is effortless to fall into the local optimal solution, and the amount of calculation is large. Practice results show that the traditional differential evolution algorithm or artificial bee colony algorithm is prone to problems such as premature maturity or slow convergence rate [35]. In addition, as the parameters of the model that needs to be identified increases, the algorithm will deteriorate the search space, and the modal error will increase with the increase in complexity. Therefore, it is difficult to efficiently and accurately search for the global optimal solution using traditional general methods. Thus, to solve this problem, this paper uses an effective artificial bee colony algorithm based on the differential evolution strategy (DE-ABC).The flowchart of DE-ABC is shown in Figure 4.
Fractional systems can be represented by fractional linear differential equations, which have the form: Among them, ( ) and ( ) are system input and system output.
If the initial condition is zero, the Laplace transform of ( ) is: Take the Laplace to transform on both ends of the above Equation (4). The fractionalorder linear time-invariant system can be rewritten as the below transfer function form: where < <···< , and < <···< .

Artificial Bee Colony Algorithm Based on Differential Evolution (DE-ABC)
Due to the complexity of the fractional systems, the dynamic model using fractionalorder description is not easy to estimate. In addition, due to the existence of multiple variables in the problem, there are multiple local search optimal solutions in the objective function, which is effortless to fall into the local optimal solution, and the amount of calculation is large. Practice results show that the traditional differential evolution algorithm or artificial bee colony algorithm is prone to problems such as premature maturity or slow convergence rate [35]. In addition, as the parameters of the model that needs to be identified increases, the algorithm will deteriorate the search space, and the modal error will increase with the increase in complexity. Therefore, it is difficult to efficiently and accurately search for the global optimal solution using traditional general methods. Thus, to solve this problem, this paper uses an effective artificial bee colony algorithm based on the differential evolution strategy (DE-ABC).The flowchart of DE-ABC is shown in Figure  4.  The calculation process of ABC mainly includes three stages: employed bees, onlooker bees and scout bees. When using ABC to identify model parameters, each food source represents a feasible solution to the model parameters, the amount of nectar(fitness) represents the quality of the solution, and the number of solutions is equal to the number of leading bees. In this paper, the root mean square error (RMSE) is introduced as the fitness function of DE-ABC to reflect the modeling error, as shown in Formula (7): Among them y exp is the sampled value of the experimental data, is the sampled value of the model data, N is the number of sampling points, each solution x i (I = 1, 2, . . . ) is represented by a D-dimensional vector x i = (x i1 , x i2 , . . . , x iD ) T , where D is the number of model parameters to be identified.
The specific process of DE-ABC identification is as follows: 1.
Obtain the input and output data of the experiment.

2.
Set the initial conditions and the number of parameters.

3.
Parameter identification based on DE-ABC.

4.
Verification. If the fitness does not meet the requirements, return to step 2 to continue identification.
DE-ABC is to introduce the mutation strategy of the differential evolution algorithm into ABC. In the search of the lead bee, the search of the lead bee is carried out according to the Formula (8). x Here, x best,j is the best individual of the previous generation, j ∈ {1, 2, . . . , D}, k ∈ {1, 2, . . . , N}, j and k are selected randomly, but k = j, F i is no longer a constant used in traditional ABC, but as shown in the Formula (9), it is an adaptive dynamic adjustment variable. Among them, Iter is the maximum number of iterations, and iter is the current number of iterations.
In the early stage of algorithm evolution, iter is smaller, F i is larger, and the algorithm mutation intensity is larger, so that it can evolve to the optimal value more effectively and quickly. As the evolution continues, to the later stage of the algorithm evolution, iter becomes larger and F i becomes smaller. As the individual evolves toward the optimal value, the function can quickly and stably converge to the optimal value.
In traditional ABC, if a certain solution x i does not improve after L cycles, this solution will be abandoned by the employed bees, and the employed bees will become scout bees and randomly generate a new solution instead. In the DE-ABC used in this article, chaotic search is introduced into the identification algorithm.
When a solution is still not improved after L cycles, it may fall into a local optimum, and the scout bees will perform a chaotic search to jump out of the local optimum. The chaotic search here uses the chaotic sequence generated by the Logistic chaotic map instead of the random number in the traditional ABC formula.
The logistic chaotic mapping equation is as follows: Among them, 0 < C n < 1. Suppose the solution of search stagnation is x i = (x i1 , x i2 , . . . , x iD ), the main steps of chaos search of the scout bees are as follows: 1.
The initial value of the chaotic sequence generated according to Formula (11); 2. Generate chaotic sequence according to Formula (10); 4.
If the maximum number of chaotic iterations is reached, the search ends, otherwise go to step 2.

Model Verification
The below paper is about that the Hammerstein model is applied to model the ratedependent hysteresis characteristics of the PEA under input signals of different frequencies.

Experimental Setup
Experimental equipment for the data acquisition experiment of PEA is shown in Figure 5. This equipment consists of a computer, a data acquisition card, a drive power supply, a piezoelectric micro-positioning platform, a piezoelectric amplification module, a piezoelectric control module, and a displacement sensor. The data acquisition card is the USB-6346(BNC) produced by NI. The piezoelectric micro-positioning platform is P733.2DD [36] produced by the PI company. The platform comes with a displacement sensor, the piezoelectric amplifier module E-505.00 produced by PI company, and the piezoelectric control module is an E-509.C2A produced by PI company. In order to obtain the experimental data, the voltage signal generated by the computer MATLAB software is outputted through USB-6346(BNC) and transmitted to the E-509.C2A piezoelectric control module. The control voltage generated is amplified by the E-505.00 amplification module as the driving voltage of the piezoelectric micro-positioning platform. At that moment, the corresponding output displacement of the platform is gauged by the displacement sensor of the P733.2DD platform and then transmitted back to Matlab via USB-6346(BNC) for storage.
2. Generate chaotic sequence according to Formula (10); 3. Generate a new solution according to Formula (12), calculate its fitness value, compare it with the original solution, and keep the best solution; 4. If the maximum number of chaotic iterations is reached, the search ends, otherwise go to step 2.

Model Verification
The below paper is about that the Hammerstein model is applied to model the ratedependent hysteresis characteristics of the PEA under input signals of different frequencies.

Experimental Setup
Experimental equipment for the data acquisition experiment of PEA is shown in Figure 5. This equipment consists of a computer, a data acquisition card, a drive power supply, a piezoelectric micro-positioning platform, a piezoelectric amplification module, a piezoelectric control module, and a displacement sensor. The data acquisition card is the USB-6346(BNC) produced by NI. The piezoelectric micro-positioning platform is P733.2DD [36] produced by the PI company. The platform comes with a displacement sensor, the piezoelectric amplifier module E-505.00 produced by PI company, and the piezoelectric control module is an E-509.C2A produced by PI company. In order to obtain the experimental data, the voltage signal generated by the computer MATLAB software is outputted through USB-6346(BNC) and transmitted to the E-509.C2A piezoelectric control module. The control voltage generated is amplified by the E-505.00 amplification module as the driving voltage of the piezoelectric micro-positioning platform. At that moment, the corresponding output displacement of the platform is gauged by the displacement sensor of the P733.2DD platform and then transmitted back to Matlab via USB-6346(BNC) for storage.

Comparison of Model Identification Effects
While the input signal frequency is lower than 5 Hz, the experimental results show that the hysteresis loop of the PEA hardly changes. Therefore, the static Bouc-Wen model is established by taking the input and output data of the PEA at an input frequency of 1 Hz, which reflects the static hysteresis characteristics of the piezoelectric actuator. In order to compare the effects of different identification algorithms, the traditional differential evolution algorithm (DE), ABC and DE-ABC are employed to discern the parameters of the improved Bouc-Wen model.
In order to realize the improved Bouc-Wen model identification process, the initial values of relevant parameters must firstly be determined. Through trial and error, the upper and lower circumscriptions of the explore scope are set to MinX = [2 0 1 1 0], MaxX = [1 −2 0 0 −1], and the number of iterations is set to 150 generations. Compared with the experimental data, the DE-ABC is superior to the DE and the ABC in terms of convergence rate and identification accuracy. The effects of the above three algorithms are shown in Figure 6 and Table 1. While the input signal frequency is lower than 5 Hz, the experimental results show that the hysteresis loop of the PEA hardly changes. Therefore, the static Bouc-Wen model is established by taking the input and output data of the PEA at an input frequency of 1 Hz, which reflects the static hysteresis characteristics of the piezoelectric actuator. In order to compare the effects of different identification algorithms, the traditional differential evolution algorithm (DE), ABC and DE-ABC are employed to discern the parameters of the improved Bouc-Wen model.
In order to realize the improved Bouc-Wen model identification process, the initial values of relevant parameters must firstly be determined. Through trial and error, the upper and lower circumscriptions of the explore scope are set to MinX = [2 0 1 1 0], MaxX = [1 −2 0 0 −1], and the number of iterations is set to 150 generations. Compared with the experimental data, the DE-ABC is superior to the DE and the ABC in terms of convergence rate and identification accuracy. The effects of the above three algorithms are shown in Figure 6 and Table 1.

Fractional Hammerstein Model
Parameter identification based on the fractional-order Hammerstein rate-dependent nonlinear hysteresis model requires two steps to complete, that is: 1. Bouc-Wen model identification. This model mainly reflects the static nonlinear hysteresis characteristics of the PEA. The model parameters and identification results are shown in Table 2 and Figure 7.

Fractional Hammerstein Model
Parameter identification based on the fractional-order Hammerstein rate-dependent nonlinear hysteresis model requires two steps to complete, that is:

1.
Bouc-Wen model identification. This model mainly reflects the static nonlinear hysteresis characteristics of the PEA. The model parameters and identification results are shown in Table 2 and Figure 7.

2.
Identification of the fractional-order model, which mainly reflects the linear dynamic characteristics of the PEA. A sine frequency sweep signal with a frequency range of 1-100 Hz is generated by Matlab as the input of the Bouc-Wen model. The output data of the Bouc-Wen model is used as the input of the fractional-order model, and the collected output data of the experimental equipment of the frequency sweep signal is used as the output of the fractional-order model. Through the input and output data of the aforementioned fractional model, the DE-ABC is used to identify the parameters of the fractional model. Through the DE-ABC, the fractional linear dynamic model is: G(s) = 5.2589 × 10 6 s 2.0384 + 3081s 1.0523 + 5.4955 × 10 6 (13)

Model Examination
In this section, the effectiveness of the rate-dependent Hammerstein model is testified when the PEA receives input signals of different frequencies.
When the input signal frequencies are 10 Hz, 20 Hz, 50 Hz and 100 Hz, respectively, the comparison between the experimentally measured hysteresis shape and the hysteresis shape simulated by the classical Bouc-Wen model, the hysteresis shape simulated by the Hammerstein model based on the integer-order dynamic model and the hysteresis shape simulated by the Hammerstein model based on the fractional-order dynamic model is shown in Figures 8-11, respectively. Table 3 shows the hysteresis shape simulated by the classic Bouc-Wen model, the hysteresis shape simulated by the Hammerstein model based on the integer-order dynamic model, and the maximum error(Max error) and the mean square error(RMSE) of the Hammerstein model based on the fractional-order dynamic model.

Model Examination
In this section, the effectiveness of the rate-dependent Hammerstein model is testified when the PEA receives input signals of different frequencies.
When the input signal frequencies are 10 Hz, 20 Hz, 50 Hz and 100 Hz, respectively, the comparison between the experimentally measured hysteresis shape and the hysteresis shape simulated by the classical Bouc-Wen model, the hysteresis shape simulated by the Hammerstein model based on the integer-order dynamic model and the hysteresis shape simulated by the Hammerstein model based on the fractional-order dynamic model is shown in Figures 8-11, respectively. Table 3          It can be known by analyzing the experimental results. The traditional Bouc-Wen model can only describe the symmetrical and rate-independent hysteresis characteristics. However, in actual conditions, the hysteresis characteristics of experimental equipment are often asymmetric and rate-dependent. It can be known from Table 3 that as the voltage frequency increases, the error of the traditional Bouc-Wen model begins to increase. The Hammerstein model based on the integer-order dynamic model can reveal the asymmetric and rate-dependent hysteresis characteristics. However, because the fractional-order model is more in line with the engineering practice and can highlight the dynamic characteristics of piezoelectric actuator more accurately, the Hammerstein model using the fractional order model to describe the dynamic characteristics of PEA is more accurate.

Conclusions
This paper proposes a rate-dependent hysteresis model based on fractional Hammerstein. The Bouc-Wen model describes the nonlinear static features of the piezoelectric actuator, and the fractional model describes the dynamic features of the piezoelectric actuator. First, the DE-ABC is used to identify the Bouc-Wen model parameters with static hysteresis. Secondly, on the basis of this Bouc-Wen static hysteresis model, the fractionalorder dynamics model is obtained through the input and output data of the sweep signal with a frequency of 1-100 Hz. The DE-ABC is applied to the parameter identification of the Bouc-Wen model and fractional order model. The simulation results show that in the wide frequency range of 1-100 Hz, the Max error is about 0.0915 µm, and the RMSE is 0.0244. In summary, when the excitation voltage frequency of the piezoelectric actuator is at different frequencies, the Hammerstein model based on the fractional-order proposed in this paper can provide higher accuracy at different frequencies.