An Asymmetric Hysteresis Model and Parameter Identification Method for Piezoelectric Actuator

Hysteresis behaviour degrades the positioning accuracy of PZT actuator for ultrahigh-precision positioning applications. In this paper, a corrected hysteresis model based on Bouc-Wen model for modelling the asymmetric hysteresis behaviour of PZT actuator is established by introducing an input bias φ and an asymmetric factor Δ Φ into the standard Bouc-Wen hysteresis model. A modified particle swarmoptimization (MPSO) algorithm is established and realized to identify and optimize themodel parameters. Feasibility and effectiveness of MPSO are proved by experiment and numerical simulation. The research results show that the corrected hysteresis model can represent the asymmetric hysteresis behaviour of the PZT actuator more accurately than the noncorrected hysteresis model based on the Bouc-Wen model.TheMPSO parameter identification method can effectively identify the parameters of the corrected and noncorrected hysteresis models. Some cases demonstrate the corrected hysteresis model and the MPSO parameter identification method can be used to model smart materials and structure systems with the asymmetric hysteresis behaviour.


Introduction
In recent years, as the requirements of the high positioning accuracy increase in the semiconductor and precision manufacturing industry, high-precision positioning devices are necessary in the applications. Piezoelectric (PZT) actuator is based on the converse piezoelectric effect [1] of piezoceramic materials and has been popularly applied as actuator in precision positioning systems [2,3], microelectronic mechanical systems [4], micro-/nanomanufacturing systems [5,6], nanobioengineering [7], and flexible electronics manufacturing [8]. The advantages of PZT actuator are fast frequency response, high positioning precision, small size, high speed, high bandwidth, high electrical-mechanical transformation efficiency, and small thermal expansion. An easy and standard way to drive a PZT actuator is to use a voltage input [9], which does not reduce the operating range and bandwidth of PZT actuator. Nevertheless, a memory dependent and multivalued relation between the output displacement and the input voltage, that is, hysteresis, is often observed in the applications. The hysteresis is referred to a complex input/output multiloop behavior; that is, the future value of output depends not only on the instantaneous value of input but also on the history of its operation, especially on the extremum. Because of hysteresis, the response (output displacement) of a PZT actuator becomes unpredictable, which drastically degrades their performance and endangers system stability in precision positioning. In order to inherit the advantages and minimize the influence of hysteresis, it is important to build the mathematical model, which can describe hysteresis accurately.
There are several mathematical models proposed to describe the hysteresis of piezoceramic materials, which can be classified into two conceptually different types. One consists of the constitutive approaches [10] that are inspired from the underlying physics of the phenomenon and are derived based on the empirical observations. The other is phenomenological approaches, which essentially employ mathematical equations to describe the phenomenon without considering the physics mechanism. The phenomenological model includes Preisach model [11,12], Maxwell model [13], neural network model [14], Prandtl-Ishlinskii model [15][16][17], and modified Prandtl-Ishlinskii model [1]. Specially, One of the more widely used phenomenological models for 2 Mathematical Problems in Engineering hysteresis systems is the Bouc-Wen model because it can capture many commonly observed types of hysteresis loops which match the behaviour of a wide class of hysteretic systems, such as piezoelectric elements, magnetorheological dampers isolation mounting system, and wood joints. Compared with other models, it has the advantages of few parameters, easy numerical simulation, and being easy to apply. However, there are two problems in previous studies on modelling hysteresis with Bouc-Wen model. First, the standard Bouc-Wen model described a symmetric hysteresis behaviour, but the PZT actuator possesses asymmetric hysteresis [18]. The symmetric hysteresis model will result in large modelling error. Second, it is not easy to determine the parameters because of the nonlinearity of the Bouc-Wen model. In order to correct the modelling error caused by asymmetric hysteresis, the Bouc-Wen model should be modified. The parameters are mostly determined within the following black-box approach: given a set of experimental input-output data, adjust the parameters so that the output of the model matches the experimental data [19]. There are several system identification methods usually employed to determine the parameters of the Bouc-Wen model, including least square method [20], genetic algorithm (GA) [21], particle swarm optimization (PSO) [22]. The PSO algorithm is simplified and observed to be performing optimization. It has been in widespread use as optimization tool in manufacture operation optimization [23], image processing [24] and robot operating system [25]. PSO algorithm is similar to GA, which also refers to the theory of evolutionary algorithm. Its idea is from random solution through iterative methods to get the best solution. Without the operator of crossover and mutation of GA, its rules are simpler and it has better convergence. Compared with other evolutionary algorithms, the main advantages of PSO are robustness in controlling parameters and high computational efficiency. Recently, it is usually employed to solve the problems of parameter optimization and get favourable effects.
In this paper, a corrected hysteresis model based on Bouc-Wen model is proposed and established in order to solve the problem of asymmetric hysteresis for the PZT actuator. Moreover, a modified particle swarm optimization (MPSO) algorithm is established and realized to identify and optimize the model parameters. The plan of this paper is organized in the following manner. In Section 2, the experimental setup is presented. In Section 3, a mathematical model is introduced to describe the hysteresis behaviour based on the Bouc-Wen model and discuss the impact of the model parameters on the hysteresis loops. The MPSO algorithm is described in detail in Section 4, discussing how to identify the system parameters of the PZT actuator. In Section 5, the performance of the corrected hysteresis model with the corresponding parameter identification method is experimentally verified and compared with the noncorrected model. Finally, conclusions are drawn in Section 6.

Experimental System
In this section, the modelling approach for hysteresis is validated experimentally on a set of data measured from a piezoelectric actuator (P-563.3CD, Physik Instrumente Co.), as shown in Figure 1. The experimental configuration is composed of the following two parts.
(2) The industrial computer with dynamic signal acquisition card (PCI-4461, National Instrument Co.) with 24-bit / converter and 24-bit / converter is used to build data acquisition system. Labview 2011 is used to process control and displacement signal in real time, and analysis of the experimental data is accomplished in MATLAB 2011b.

Standard Bouc-Wen Model.
The Bouc-Wen model is first proposed by Bouc in 1967 [26], where a function that describes the hysteresis behaviour between the displacement and restoring force was proposed, and it was then generalized by Wen in 1976 [27]. This model consists of a first-order nonlinear state equation and an output equation where the input and state signals appear linearly. Through appropriate choices of parameters in the model, it can represent wide variety of hysteresis behaviours. Considering a physical system with a hysteresis behaviour, the Bouc-Wen model represents the hysteresis between the output Φ BW ( , ) and input ( ) which is shown in the following expressions [28]: The standard Bouc-Wen model represents the output Φ BW ( , ) by the superposition of a linear component ( ) and a nonlinear hysteresis component (1− ) ℎ( ). Where is the stiffness coefficient; 0 < < 1 is a weighting parameter; , , , , and are the model parameters; ℎ is the hysteretic variable; andḣ denotes the time derivative.
Hysteresis loops with different parameters value of , , , and under the input signal ( ) = Amp ⋅ sin( ) at Amp = 2 are shown in Figure 2. The relationship between and and their effects on hysteresis with different amplitudes input signal ( ) = Amp ⋅ sin( ) at several values of Amp, 0.5, 1, and 2, are shown in Figure 3, and the Bouc-Wen model parameters are listed in Table 1.
According to Figures 2 and 3, the Bouc-Wen model can represent a wide variety of hysteresis loops by appropriate choices of different parameters. Actually, the shape of the hysteresis loop depends only on the parameters , , , and .     Figure 3.   Thereby, , , and are the parameters for the control of basic hysteresis shape and for sharpness to influence the loop smoothness. However, these parameters are functionally redundant in some applications; researchers usually fix some parameters to constant value [29,30]. But, the parameters value of , , , , and cannot be set up arbitrary in order to model a real physical system with hysteresis. The restrictions for the parameters in the standard Bouc-Wen model are shown in Table 2 [31].

Hysteresis
Model of PZT Actuator. The nanopositioning platform driven by a PZT actuator was constructed as shown in Figure 1; the physical model of the one-axes can be equivalent to the following spring-mass-damper system.
The one-axes of PZT actuator can be depicted in Figure 4. Applying an input voltage , an elongation is produced, which results in a force PZT acting on the sliding block and causes the displacement . The electromechanical equation can be obtained according to the Newton's laws of motion: where is the equivalent mass, is relative position with respect to the capacitive displacement sensor, 0 is a constant that depends on the choice of the origin. The relationship between the actuated force PZT and the applied voltage is not linear but nonlinear due to the hysteresis behaviour of PZT actuator; it can be formulated by the Bouc-Wen model: There are eleven parameters in (3) and (4). In order to simplify the calculation, it is necessary to fix some parameters to constants and simplify these equations without deviating from physical reality. First, the constant 0 is an initial displacement that can be eliminated by zero point calibrating, so it can be set to zero. Second, in the applications of precision measurement and positioning, the input voltage is usually using a periodic voltage with a low frequency in order to prevent the mechanical vibration and overheating. In this situation, the terms̈anḋcan be neglected. Then putting the simplified equation (4) into (3) can result in Because the parameters , , , and are constant, (5) can be further simplified. Let Because the parameter > 0, the sign of 1 , 1 , and 1 is the same as , , and and is bound by constraint condition of Table 2. The physical description for the hysteresis system of PZT actuator can be written in the following forms, representing the relationship of the output displacement to its input voltage: Through the simplification and constraints, the total number of unknown parameters has reduced to six, respectively, being 1 , 2 , 1 , 1 , 1 , and .

Modified Particle Swarm Optimization (MPSO).
PSO is an optimization approach based on stochastic population and was originally attributed to Kennedy and Eberhart [32]. And then, a modified particle swarm optimization (MPSO) algorithm with inertia weight factor was proposed by Shi and Eberhart [33] in order to improve the performance of algorithm.
In a -dimensional parametric search space, the population size of particle is defined as Size; each particle represents a candidate solution in the solution space; the position and velocity of the particle no. ( = 1, 2, . . . , Size) are defined as and V . All the particles are evaluated by the predefined fitness function ( ). Comparing the fitness function, each particle records its individual best position as best and the global best position as best. The maximum number of iteration is defined as , and ( = 1, 2, . . . , ) is the current number of iteration. Then the velocity and position of the particle no. at iteration are updated by the following functions: where is the inertia weighting factor, 1 and 2 are local acceleration coefficients and global acceleration coefficients, respectively, and rand 1 () and rand 2 () are random numbers in the range of [0, 1].
For the particle , the best personal and global positions are calculated by the following equations (on the basis of minimum value of fitness function, maximum value in the same way) Equation (9) indicates that the velocity of a particle is updated according to three terms. The first is previous velocity V , scaled by , and it is often known as habitual behaviour. The second is previous best position best , scaled by the 1 and rand 1 (), and it is often known as cognition model or self-memory model. The third is the global best position best , scaled by the 2 and rand 2 (); this term is often known as social model. The optimization mechanism of MPSO is shown in Figure 5. For the particle , path A of Figure 5 reflects the impact of current velocity in MPSO; path B reflects the impact of cognition model or self-memory model; path C reflects the impact of the swarm or social model; and these correspond to the first term, second term, and third term in (9), respectively.

Hysteresis Model Parameters Identification with MPSO.
The principle of parameters identification of the hysteresis  model for PZT actuator using MPSO algorithm can be described in Figure 6.
The hysteresis model as described in (7) and (8), = [ 1 , 2 , 1 , 1 , 1 , ] is a set of estimated parameters. Suppose that the search space is -dimensional, the particle (1 ≤ ≤ Size) of the swarm can be represented by a -dimensional vector , = ( ,1 , ,2 , . . . , , ), each dimension of particle swarm represents an unknown parameter. The velocity of particle can represent V , = (V ,1 , V ,2 , . . . , V , ), the individual best position best , = ( ,1 , ,2 , . . . , , ), the global best position best = ( 1 , 2 , . . . , ). The fitness is evaluated as the root-mean-square error between the actual measured displacement and the model output; that is, where is number of samples, ( ) is the real measured value, and̂( ) is the model predicted output value. The velocity and position of the particle no. at iteration are updated by (9) and (10). The inertia weighting factor is responsible for the definition of the searching area, which provides a balance between the global and local explorations. Usually, a bigger is advisable at the beginning of the searching procedure, rather than after some iterations. In this paper, the inertia weighting factor is decreased linearly Mathematical Problems in Engineering during the optimization, which is defined as the following function: where ( ) is the inertia weighting factor at the th iteration and max and min are the maximum and minimum value. If ( ) = 0, the velocity depends only on best and best, with no memory effect. Else if ( ) ̸ = 0, the value of decreases linearly from max to min with the increasing of iteration, which makes the MPSO algorithm explore new areas in the initial stage. In the later stage, it makes the MPSO algorithm conduct refined searching around the optimum solution. Figure 7.

Algorithm Implementation. The flowchart of the identification model parameters with MPSO is shown in
The parameter identification procedures consist of three parts: Part I, the particle swarm initialization; Part II, fitness evaluation; and Part III, swarm optimizing, which are expressed in detail as follows.

Part I: The Particle Swarm Initialization
Step 1. Set the appropriate population size Size, the maximum number of iteration , and the value of local acceleration coefficients 1 and global acceleration coefficients 2 .
Step 2. Determine the maximum and minimum value of inertia weighting factor max and min .
Step 3. Set the position boundary value [ min , max ] and velocity boundary value [V min , V max ]. The position boundary value is the range of each unknown parameter.
Step 4. Initialize a population of particles in -dimensions searching space. Randomly initialize the position vector and velocity vector of each particle according to the following equations:

Part II: Fitness Evaluation
Step 5. Acquire the output of PZT actuator and determine the number of samples of data . The higher the value , the higher the credibility of the fitness function, but the higher the amount of calculation.
Step 6. Evaluate the fitness function of each particle according to (12).

Part III: Swarm Optimizing
Step 7. Find the best position best and the global position best of individual point according to (11).
Step 8. Update the position and velocity using (9) and (10); check whether or not the position and velocity are beyond boundary.
Step 9. Check the exit condition If = , output the results of parameter identification, and terminate the MPSO algorithm; otherwise, jump to Step 6.

Experimental Design and MPSO Algorithm Parameters
Selection. The experimental system is depicted in Section 2. It considers that the input voltage of voltage amplifier is cosine voltage with the frequency by 0.1 Hz with the form of 1 ( ) = 2 − 2 cos(0.2 ) to drive the precision positioning platform; then the experiment data are obtained by the data sampling system. Because the general property of these phenomenological models is that the number of sample data determines the accuracy of the models, the sampling frequency is set to 1000 Hz; this means that the single-period sample size is  Table 3.
Equations (7) and (8) of the hysteresis model cannot be used directly in the simulation; the fourth-order Runge-Kutta method [34] is adopted by MATLAB 2011b, in order to obtain the numerical results.

Simulation Result and Analysis.
In order to analyse the accuracy of hysteresis model and quantify modelling error, the maximum absolute deviation (MAD), the maximum relative deviation (MRD), the mean absolute error (MAE), and the normalized root-mean-square deviation (NRMSD) are defined as follows: where Max is the maximum value of the measured output displacement. Among them, the MAD and MRD are used to evaluate local error accuracy and the MAE and NRMSD are used to evaluate global error accuracy. The results of parameter identification of the of MPSO algorithm are listed in Table 4. The hysteresis loops of voltagedisplacement responses are shown in Figure 8 under the input of 1 ( ), within which, the simulation and experimental results are compared.
Observing Figure 8, the modelling error of Case 1 (Figure 8(a)) and Case 2 (Figure 8 than Case 3 (Figure 8(c)) and Case 4 (Figure 8(d)). With the increasing of iterations and the population size, the simulation output curves of hysteresis model approximate the experimental curves gradually. Furthermore, the MAD in Case 3 and Case 4 are 2.6782 um and 0.7726 um; the MRD are 2.2318% and 0.6438% of the full displacement range of 0∼120 um.
Simulation results show that the parameter identification method based on the MPSO algorithm is feasible and effective. But the modelling error apparently increases with the increase in the applied voltage, and the maximum error appears in the maximum applied voltage area. Although, the results of Case 3 and Case 4 are very close to the real hysteresis loop in shape, there are still modelling error, characterized Table 4: List of the identified parameters of hysteresis model.  by the phase deviation and asymmetric shape deviation. The mainly reason of modelling error is that the hysteresis behaviour of PZT actuator is asymmetric, but the Bouc-Wen model describes the symmetric hysteresis. It is necessary to correct the phase deviation and asymmetric shape deviation to improve the hysteresis model accuracy. (8) can be manipulated and rewritten as follows:

Error Correction of Hysteresis Model. Equation
The above equation is modified by adding input bias and the asymmetric factor Δ Φ to correct the phase deviation and asymmetric shape deviation, respectively, in order to eliminate the modelling error, hence called the corrected hysteresis model. Accordingly, the corrected hysteresis model of the PZT actuator can be rewritten as follows: where the is asymmetric coefficient. In the voltage rising process, the sign(( )) = 1, if > 0; the higher the applied voltage ( ), the more apparently the asymmetric factor Δ Φ positive effects on ℎ( ); and the < 0 causes the same trend of negative effects. In the voltage descending process, the sign(( )) = −1; these trends are exactly opposite. Moreover, the input bias can direct causes phase offset of output displacement ( ). The asymmetric factor Δ Φ and the input bias cause the model with the asymmetry and the deflection of hysteresis loop not to change its BIBO stability, so, the parameters 1 , 1 , 1 , and of (17) are still bounded by constraint condition of Table 2.
For verifying the accuracy of the corrected hysteresis model, Case 3 and Case 4 are still used. Firstly, fix parameters  Table 5.
Figures 9(a) and 9(b) show that the fitting precision of the corrected hysteresis model is significantly better than the previous noncorrected model. According   0.4303% of the full displacement range of 0∼120 um, and the MAE and the NRMSD by the corrected model are reduced to 0.2340 um, 0.0023 and 0.2200 um, 0.0023, respectively. Figure 10 shows that the maximum absolute error by the noncorrected model generally occurs in the peak values of the applied voltage, but the absolute error has been effectively inhibited by the corrected hysteresis model, which indicates the input bias and the asymmetric factor Δ Φ improves the model precision effectively. Consider a 0.1 HZ multiperiod input voltage with the form of 2 ( ) = 4 − 4 cos(0.2 ), and use the MPSO to identify the parameters [ 1 , 2 , 1 , 1 , 1 , , , ] for (17) with the Size = 300 and = 8000. The results of MPSO algorithm are 1 = 39.3383 * 10 −6 , 2 = 0.2412, 1 = −15.3487 * 10 −6 , 1 = 9.2620, 1 = −6.4160, = 2, = 0.8612, and = 1.05 × 10 −8 , and the single-period hysteresis loops of the PZT actuator obtained by experimental data and the model response are shown in Figure 11, the time histories of the output displacements are shown in Figure 12, the trends of absolute error are shown in Figure 13, and quantitative error of the output displacements predicted is shown in Table 6.
Then, consider a 0.1 Hz input voltage with decreasing amplitude; the time histories of the output displacements from the PZT actuator and the response of the two models are shown in Figure 14.
Observing Figures 11, 12, and 13, the modelling error can be effectively reduced by the corrected hysteresis model. The responses of the corrected hysteresis model are very close to the real system and are in good agreement with    Figure 14, under the applied voltage with decayed amplitude, the modelling error can be effectively reduced by the corrected hysteresis model, and the higher the amplitude is, the more evident the effect is. According to the above simulations and analyses, it can be determined that the corrected hysteresis model can achieve a higher accuracy when modelling the hysteresis of the PZT actuator due to introducing an input bias and an asymmetric factor Δ Φ to represent the asymmetric hysteresis behaviour for PZT actuator. However, experiments found that the strength of the hysteresis behaviour is largely related to the rate of change applied voltage, and the hysteresis loop will be clockwise rotated with increasing the frequency of applied voltage, that is, rate-dependent hysteresis. Therefore, the further work will focus on how to model the ratedependent hysteresis and constructing adaptive control system for the precise positioning control of the PZT actuator.

Conclusion
In this paper, a corrected hysteresis model based on Bouc-Wen model for modelling the hysteresis behaviour of PZT actuator is established by introducing an input bias and an asymmetric factor Δ Φ into the standard Bouc-Wen hysteresis model. Moreover, a modified particle swarm optimization (MPSO) algorithm was established and realized to identify and optimize the model parameters. Feasibility and effectiveness of MPSO are proved by experiment and numerical simulation. The research results have shown that the corrected hysteresis model can represent the asymmetric hysteresis behaviour of the PZT actuator more accurately than the general hysteresis based on the Bouc-Wen model. The input bias and the asymmetric factor Δ Φ correct the phase deviation and asymmetric shape deviation and improve the model precision effectively. The MPSO parameter identification method can effectively and quickly identify the parameters of the corrected and noncorrected hysteresis models for the PZT actuator. Although, only a few cases for modelling the PZT actuator are presented, the corrected hysteresis model with the MPSO parameter identification method can be also used to model smart materials and structure systems with the asymmetric hysteresis behaviour.