# An Asymmetric Hysteresis Model and Parameter Identification Method for Piezoelectric Actuator

Authors: Haichen Qin; Ningbin Bu; Wei Chen; Zhouping Yin
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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 swarm optimization (MPSO) algorithm is 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 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. The MPSO 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.


---

## Body

## 1. 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–17], and modified Prandtl-Ishlinskii model [1]. Specially, One of the more widely used phenomenological models for 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 Section2, 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.


## 2. 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 Figure1. The experimental configuration is composed of the following two parts.(1)
The piezoelectric actuator (P-563.3CD, Physik Instrumente Co.) isX-Y-Z high-precision nanopositioning stage with the nominal displacement expansion of 0~340 μm under the input voltage 0~120 V on each axis; the voltage amplifier (E503.00, Physik Instrumente Co.) with input signal voltage is 0~12 V and output is 0~120 V; capacitive sensor/controller module (E-509.C3A, Physik Instrumente Co.) is used in the experiments.
(2)
The industrial computer with dynamic signal acquisition card (PCI-4461, National Instrument Co.) with 24-bitA
/
D converter and 24-bit D
/
A 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.The experimental system: (a) structure of the experimental system and (b) the experimental equipment.
(a)
(b)


## 3. Hysteresis Modelling of PZT Actuator

### 3.1. 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
(
u
,
t
) and input u
(
t
) which is shown in the following expressions [28]:
(1)
Φ
BW
(
u
,
t
)
=
α
k
u
(
t
)
+
(
1
-
α
)
D
k
h
(
t
)
,
(2)
h
˙
(
t
)
=
D
-
1
(
A
u
˙
(
t
)
-
β
|
u
˙
(
t
)
|
|
h
(
t
)
|
n
-
1
h
(
t
)
-
γ
u
˙
(
t
)
|
h
(
t
)
|
n
)
.The standard Bouc-Wen model represents the outputΦ
BW
(
u
,
t
) by the superposition of a linear component α
k
u
(
t
) and a nonlinear hysteresis component (
1
-
α
)
D
k
h
(
t
). Where k is the stiffness coefficient; 0
<
α
<
1 is a weighting parameter; D, A, β, γ, and n are the model parameters; h is the hysteretic variable; and h
˙ denotes the time derivative.Hysteresis loops with different parameters value ofA, β, γ, and n under the input signal u
(
t
)
=
Amp
·
sin
(
ω
t
) at Amp
=
2 are shown in Figure 2. The relationship between β and γ and their effects on hysteresis with different amplitudes input signal u
(
t
)
=
Amp
·
sin
(
ω
t
) 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.Table 1
Bouc-Wen model parameters in Figure3.
D
A
β
γ
n
β
+
γ
β
-
γ
Figure3(a)
1
1
0.6
0.4
1
>0
>0
Figure3(b)
1
1
0.6
0.6
1
>0
=0
Figure3(c)
1
1
0.6
−0.6
1
=0
>0
Figure3(d)
1
1
0.4
0.6
1
>0
<0
Figure3(e)
1
1
0.4
−0.6
1
<0
>0
Figure3(f)
1
1
0.8
0.2
1
=1
>0
0.5
0.5
=0
0.4
0.6
<0Hysteresis loops via the differential equation (2) with different parameters value of  A, β, γ, and n; the input signal u
(
t
)
=
Amp
·
sin
⁡
(
ω
t
), Amp
=
2. (a) The shape of hysteresis loops with several values of A. (b) The shape of hysteresis loops with several values of n. (c) The shape of hysteresis loops with several values of β. (d) The shape of hysteresis loops with several values of γ.
(a)
(b)
(c)
(d)Hysteresis loops via (2) under different amplitudes input signal with fixed parameters D
=
1, A
=
1, and  n
=
1 and several values of β and γ; the input signal u
(
t
)
=
Amp
·
sin
⁡
(
ω
t
), Amp
=
[
0.5,1
,
2
].
(a)
(b)
(c)
(d)
(e)
(f)According to Figures2 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 A, β, γ, and n. Thereby, A, β, and γ are the parameters for the control of basic hysteresis shape and n 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 D, A, β, γ, and n 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].Table 2
Classification of stable Bouc-Wen model (D
>
0 and n
≥
1).
A
β and γ
A
>
0
β
+
γ
>
0 and β
-
γ
≥
0
β
-
γ
<
0 and β
≥
0
A
<
0
β
-
γ
>
0 and β
+
γ
≥
0
β
+
γ
<
0 and β
≥
0


### 3.2. Hysteresis Model of PZT Actuator

The nanopositioning platform driven by a PZT actuator was constructed as shown in Figure1; 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 Figure4. Applying an input voltage u, an elongation is produced, which results in a force F
PZT acting on the sliding block m and causes the displacement  y. The electromechanical equation can be obtained according to the Newton’s laws of motion:
(3)
m
y
¨
+
c
y
˙
+
k
s
(
y
-
y
0
)
=
F
PZT
,
where m is the equivalent mass, y is relative position with respect to the capacitive displacement sensor, y
0 is a constant that depends on the choice of the origin. The relationship between the actuated force F
PZT and the applied voltage u is not linear but nonlinear due to the hysteresis behaviour of PZT actuator; it can be formulated by the Bouc-Wen model:
(4)
F
PZT
=
α
k
u
+
(
1
-
α
)
D
k
h
,
h
˙
=
D
-
1
(
A
u
˙
-
β
|
u
˙
|
h
|
h
|
n
-
1
-
γ
u
˙
|
h
|
n
)
.Figure 4
The schematic diagram of precision moving stage driven by a PZT actuator.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 y
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 m
y
¨ and c
y
˙ can be neglected. Then putting the simplified equation (4) into (3) can result in(5)
y
=
α
k
u
k
s
+
(
1
-
α
)
k
D
k
s
h
,
h
˙
=
A
D
u
˙
-
β
D
|
u
˙
|
h
|
h
|
n
-
1
-
γ
D
u
˙
|
h
|
n
.Because the parametersα, k, k
s, and D are constant, (5) can be further simplified.Let(6)
k
1
=
α
k
k
s
>
0
,
k
2
=
(
1
-
α
)
k
D
k
s
>
0
,
A
1
=
A
D
,
β
1
=
β
D
,
γ
1
=
γ
D
.Because the parameterD
>
0, the sign of A
1, β
1, and γ
1 is the same as A, β, 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:(7)
y
=
k
1
u
(
t
)
+
k
2
h
(
t
)
,
(8)
h
˙
(
t
)
=
A
1
u
˙
(
t
)
-
β
1
|
u
˙
(
t
)
|
h
(
t
)
|
h
(
t
)
|
n
-
1
-
γ
1
u
˙
(
t
)
|
h
(
t
)
|
n
.Through the simplification and constraints, the total number of unknown parameters has reduced to six, respectively, beingk
1, k
2, A
1, β
1, γ
1, and n.


## 3.1. 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
(
u
,
t
) and input u
(
t
) which is shown in the following expressions [28]:
(1)
Φ
BW
(
u
,
t
)
=
α
k
u
(
t
)
+
(
1
-
α
)
D
k
h
(
t
)
,
(2)
h
˙
(
t
)
=
D
-
1
(
A
u
˙
(
t
)
-
β
|
u
˙
(
t
)
|
|
h
(
t
)
|
n
-
1
h
(
t
)
-
γ
u
˙
(
t
)
|
h
(
t
)
|
n
)
.The standard Bouc-Wen model represents the outputΦ
BW
(
u
,
t
) by the superposition of a linear component α
k
u
(
t
) and a nonlinear hysteresis component (
1
-
α
)
D
k
h
(
t
). Where k is the stiffness coefficient; 0
<
α
<
1 is a weighting parameter; D, A, β, γ, and n are the model parameters; h is the hysteretic variable; and h
˙ denotes the time derivative.Hysteresis loops with different parameters value ofA, β, γ, and n under the input signal u
(
t
)
=
Amp
·
sin
(
ω
t
) at Amp
=
2 are shown in Figure 2. The relationship between β and γ and their effects on hysteresis with different amplitudes input signal u
(
t
)
=
Amp
·
sin
(
ω
t
) 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.Table 1
Bouc-Wen model parameters in Figure3.
D
A
β
γ
n
β
+
γ
β
-
γ
Figure3(a)
1
1
0.6
0.4
1
>0
>0
Figure3(b)
1
1
0.6
0.6
1
>0
=0
Figure3(c)
1
1
0.6
−0.6
1
=0
>0
Figure3(d)
1
1
0.4
0.6
1
>0
<0
Figure3(e)
1
1
0.4
−0.6
1
<0
>0
Figure3(f)
1
1
0.8
0.2
1
=1
>0
0.5
0.5
=0
0.4
0.6
<0Hysteresis loops via the differential equation (2) with different parameters value of  A, β, γ, and n; the input signal u
(
t
)
=
Amp
·
sin
⁡
(
ω
t
), Amp
=
2. (a) The shape of hysteresis loops with several values of A. (b) The shape of hysteresis loops with several values of n. (c) The shape of hysteresis loops with several values of β. (d) The shape of hysteresis loops with several values of γ.
(a)
(b)
(c)
(d)Hysteresis loops via (2) under different amplitudes input signal with fixed parameters D
=
1, A
=
1, and  n
=
1 and several values of β and γ; the input signal u
(
t
)
=
Amp
·
sin
⁡
(
ω
t
), Amp
=
[
0.5,1
,
2
].
(a)
(b)
(c)
(d)
(e)
(f)According to Figures2 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 A, β, γ, and n. Thereby, A, β, and γ are the parameters for the control of basic hysteresis shape and n 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 D, A, β, γ, and n 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].Table 2
Classification of stable Bouc-Wen model (D
>
0 and n
≥
1).
A
β and γ
A
>
0
β
+
γ
>
0 and β
-
γ
≥
0
β
-
γ
<
0 and β
≥
0
A
<
0
β
-
γ
>
0 and β
+
γ
≥
0
β
+
γ
<
0 and β
≥
0


## 3.2. Hysteresis Model of PZT Actuator

The nanopositioning platform driven by a PZT actuator was constructed as shown in Figure1; 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 Figure4. Applying an input voltage u, an elongation is produced, which results in a force F
PZT acting on the sliding block m and causes the displacement  y. The electromechanical equation can be obtained according to the Newton’s laws of motion:
(3)
m
y
¨
+
c
y
˙
+
k
s
(
y
-
y
0
)
=
F
PZT
,
where m is the equivalent mass, y is relative position with respect to the capacitive displacement sensor, y
0 is a constant that depends on the choice of the origin. The relationship between the actuated force F
PZT and the applied voltage u is not linear but nonlinear due to the hysteresis behaviour of PZT actuator; it can be formulated by the Bouc-Wen model:
(4)
F
PZT
=
α
k
u
+
(
1
-
α
)
D
k
h
,
h
˙
=
D
-
1
(
A
u
˙
-
β
|
u
˙
|
h
|
h
|
n
-
1
-
γ
u
˙
|
h
|
n
)
.Figure 4
The schematic diagram of precision moving stage driven by a PZT actuator.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 y
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 m
y
¨ and c
y
˙ can be neglected. Then putting the simplified equation (4) into (3) can result in(5)
y
=
α
k
u
k
s
+
(
1
-
α
)
k
D
k
s
h
,
h
˙
=
A
D
u
˙
-
β
D
|
u
˙
|
h
|
h
|
n
-
1
-
γ
D
u
˙
|
h
|
n
.Because the parametersα, k, k
s, and D are constant, (5) can be further simplified.Let(6)
k
1
=
α
k
k
s
>
0
,
k
2
=
(
1
-
α
)
k
D
k
s
>
0
,
A
1
=
A
D
,
β
1
=
β
D
,
γ
1
=
γ
D
.Because the parameterD
>
0, the sign of A
1, β
1, and γ
1 is the same as A, β, 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:(7)
y
=
k
1
u
(
t
)
+
k
2
h
(
t
)
,
(8)
h
˙
(
t
)
=
A
1
u
˙
(
t
)
-
β
1
|
u
˙
(
t
)
|
h
(
t
)
|
h
(
t
)
|
n
-
1
-
γ
1
u
˙
(
t
)
|
h
(
t
)
|
n
.Through the simplification and constraints, the total number of unknown parameters has reduced to six, respectively, beingk
1, k
2, A
1, β
1, γ
1, and n.


## 4. Model Parameters Identification

### 4.1. 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 ad-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. i   (i
=
1,2
,
…
,
Size) are defined as x
i and v
i. All the particles are evaluated by the predefined fitness function f
(
x
i
). Comparing the fitness function, each particle records its individual best position as p
best
i and the global best position as g
best. The maximum number of iteration is defined as G, and k
g   (k
g
=
1,2
,
…
,
G) is the current number of iteration. Then the velocity and position of the particle no. i at iteration k
g are updated by the following functions:(9)
v
i
k
g
+
1
=
w
·
v
i
k
g
+
c
1
·
rand
1
(
)
·
(
p
best
i
k
g
-
x
i
k
g
)
+
c
2
·
rand
2
(
)
·
(
g
best
k
g
-
x
i
k
g
)
,
(10)
x
i
k
g
+
1
=
x
i
k
g
+
v
i
k
g
+
1
,
where w is the inertia weighting factor, c
1 and c
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 particlei, 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)
(11)
p
best
i
k
g
+
1
=
{
p
best
i
k
g
,
f
(
x
i
k
g
+
1
)
≥
f
(
p
best
i
k
g
)
,
x
i
k
g
+
1
,
f
(
x
i
k
g
+
1
)
<
f
(
p
best
i
k
g
)
,
g
best
k
g
+
1
=
min
⁡
[
f
(
p
best
i
k
g
)
]
.Equation (9) indicates that the velocity of a particle is updated according to three terms. The first is previous velocity v
i
k
g, scaled by w, and it is often known as habitual behaviour. The second is previous best position p
best
i
k
g, scaled by the c
1 and rand
1
(
), and it is often known as cognition model or self-memory model. The third is the global best position g
best
k
g, scaled by the c
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 x
i, 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.Figure 5
The optimization mechanism of PSO.


### 4.2. 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 Figure6.Figure 6
Parameters identification of Hysteresis model.The hysteresis model as described in (7) and (8), θ
=
[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
] is a set of estimated parameters. Suppose that the search space is d-dimensional, the particle i (1
≤
i
≤
Size) of the swarm can be represented by a d-dimensional vector x
i
,
j
=
(
x
i
,
1
,
x
i
,
2
,
…
,
x
i
,
d
), each dimension of particle swarm represents an unknown parameter. The velocity of particle i can represent v
i
,
j
=
(
v
i
,
1
,
v
i
,
2
,
…
,
v
i
,
d
), the individual best position p
best
i
,
j
=
(
p
i
,
1
,
p
i
,
2
,
…
,
p
i
,
d
), the global best position g
best
j
=
(
g
1
,
g
2
,
…
,
g
d
). The fitness is evaluated as the root-mean-square error between the actual measured displacement and the model output; that is,(12)
f
(
x
i
)
=
1
N
∑
l
=
1
N
[
y
(
l
)
-
y
^
(
l
)
]
2
,
where N is number of samples, y
(
l
) is the real measured value, and y
^
(
l
) is the model predicted output value.The velocity and position of the particle no.i at iteration k
g are updated by (9) and (10). The inertia weighting factor w is responsible for the definition of the searching area, which provides a balance between the global and local explorations. Usually, a bigger w is advisable at the beginning of the searching procedure, rather than after some iterations. In this paper, the inertia weighting factor is decreased linearly during the optimization, which is defined as the following function:
(13)
w
(
k
g
)
=
w
max
⁡
-
k
g
G
(
w
max
⁡
-
w
min
⁡
)
,
where w
(
k
g
) is the inertia weighting factor at the k
gth iteration and w
max
⁡ and w
min
⁡ are the maximum and minimum value. If w
(
k
g
)
=
0, the velocity depends only on pbest and gbest, with no memory effect. Else if w
(
k
g
)
≠
0, the value of w decreases linearly from w
max
⁡ to w
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.


### 4.3. Algorithm Implementation

The flowchart of the identification model parameters with MPSO is shown in Figure7.Figure 7
Flowchart of the identification hysteresis model parameters with MPSO.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.


#### 4.3.1. Part I: The Particle Swarm Initialization

Step 1.
Set the appropriate population size Size, the maximum number of iterationG, and the value of local acceleration coefficients c
1 and global acceleration coefficients c
2.Step 2.
Determine the maximum and minimum value of inertia weighting factorw
max
⁡ and w
min
⁡.Step 3.
Set the position boundary value[
x
min
⁡
,
x
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 ind-dimensions searching space. Randomly initialize the position vector and velocity vector of each particle according to the following equations:
(14)
x
(
i
,
j
)
=
x
(
i
,
j
)
min
⁡
+
(
x
(
i
,
j
)
max
⁡
-
x
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
v
(
i
,
j
)
=
v
(
i
,
j
)
min
⁡
+
(
v
(
i
,
j
)
max
⁡
-
v
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
p
p
p
p
p
p
p
i
=
[
1,2
,
…
,
Size
]
,
j
=
[
1,2
,
…
,
d
]
.


#### 4.3.2. Part II: Fitness Evaluation

Step 5.
Acquire the output of PZT actuator and determine the number of samples of dataN. The higher the value N, 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).


#### 4.3.3. Part III: Swarm Optimizing

Step 7.
Find the best positionpbest and the global position gbest 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 Ifk
g
=
G, output the results of parameter identification, and terminate the MPSO algorithm; otherwise, jump to Step 6.


## 4.1. 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 ad-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. i   (i
=
1,2
,
…
,
Size) are defined as x
i and v
i. All the particles are evaluated by the predefined fitness function f
(
x
i
). Comparing the fitness function, each particle records its individual best position as p
best
i and the global best position as g
best. The maximum number of iteration is defined as G, and k
g   (k
g
=
1,2
,
…
,
G) is the current number of iteration. Then the velocity and position of the particle no. i at iteration k
g are updated by the following functions:(9)
v
i
k
g
+
1
=
w
·
v
i
k
g
+
c
1
·
rand
1
(
)
·
(
p
best
i
k
g
-
x
i
k
g
)
+
c
2
·
rand
2
(
)
·
(
g
best
k
g
-
x
i
k
g
)
,
(10)
x
i
k
g
+
1
=
x
i
k
g
+
v
i
k
g
+
1
,
where w is the inertia weighting factor, c
1 and c
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 particlei, 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)
(11)
p
best
i
k
g
+
1
=
{
p
best
i
k
g
,
f
(
x
i
k
g
+
1
)
≥
f
(
p
best
i
k
g
)
,
x
i
k
g
+
1
,
f
(
x
i
k
g
+
1
)
<
f
(
p
best
i
k
g
)
,
g
best
k
g
+
1
=
min
⁡
[
f
(
p
best
i
k
g
)
]
.Equation (9) indicates that the velocity of a particle is updated according to three terms. The first is previous velocity v
i
k
g, scaled by w, and it is often known as habitual behaviour. The second is previous best position p
best
i
k
g, scaled by the c
1 and rand
1
(
), and it is often known as cognition model or self-memory model. The third is the global best position g
best
k
g, scaled by the c
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 x
i, 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.Figure 5
The optimization mechanism of PSO.


## 4.2. 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 Figure6.Figure 6
Parameters identification of Hysteresis model.The hysteresis model as described in (7) and (8), θ
=
[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
] is a set of estimated parameters. Suppose that the search space is d-dimensional, the particle i (1
≤
i
≤
Size) of the swarm can be represented by a d-dimensional vector x
i
,
j
=
(
x
i
,
1
,
x
i
,
2
,
…
,
x
i
,
d
), each dimension of particle swarm represents an unknown parameter. The velocity of particle i can represent v
i
,
j
=
(
v
i
,
1
,
v
i
,
2
,
…
,
v
i
,
d
), the individual best position p
best
i
,
j
=
(
p
i
,
1
,
p
i
,
2
,
…
,
p
i
,
d
), the global best position g
best
j
=
(
g
1
,
g
2
,
…
,
g
d
). The fitness is evaluated as the root-mean-square error between the actual measured displacement and the model output; that is,(12)
f
(
x
i
)
=
1
N
∑
l
=
1
N
[
y
(
l
)
-
y
^
(
l
)
]
2
,
where N is number of samples, y
(
l
) is the real measured value, and y
^
(
l
) is the model predicted output value.The velocity and position of the particle no.i at iteration k
g are updated by (9) and (10). The inertia weighting factor w is responsible for the definition of the searching area, which provides a balance between the global and local explorations. Usually, a bigger w is advisable at the beginning of the searching procedure, rather than after some iterations. In this paper, the inertia weighting factor is decreased linearly during the optimization, which is defined as the following function:
(13)
w
(
k
g
)
=
w
max
⁡
-
k
g
G
(
w
max
⁡
-
w
min
⁡
)
,
where w
(
k
g
) is the inertia weighting factor at the k
gth iteration and w
max
⁡ and w
min
⁡ are the maximum and minimum value. If w
(
k
g
)
=
0, the velocity depends only on pbest and gbest, with no memory effect. Else if w
(
k
g
)
≠
0, the value of w decreases linearly from w
max
⁡ to w
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.


## 4.3. Algorithm Implementation

The flowchart of the identification model parameters with MPSO is shown in Figure7.Figure 7
Flowchart of the identification hysteresis model parameters with MPSO.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.


### 4.3.1. Part I: The Particle Swarm Initialization

Step 1.
Set the appropriate population size Size, the maximum number of iterationG, and the value of local acceleration coefficients c
1 and global acceleration coefficients c
2.Step 2.
Determine the maximum and minimum value of inertia weighting factorw
max
⁡ and w
min
⁡.Step 3.
Set the position boundary value[
x
min
⁡
,
x
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 ind-dimensions searching space. Randomly initialize the position vector and velocity vector of each particle according to the following equations:
(14)
x
(
i
,
j
)
=
x
(
i
,
j
)
min
⁡
+
(
x
(
i
,
j
)
max
⁡
-
x
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
v
(
i
,
j
)
=
v
(
i
,
j
)
min
⁡
+
(
v
(
i
,
j
)
max
⁡
-
v
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
p
p
p
p
p
p
p
i
=
[
1,2
,
…
,
Size
]
,
j
=
[
1,2
,
…
,
d
]
.


### 4.3.2. Part II: Fitness Evaluation

Step 5.
Acquire the output of PZT actuator and determine the number of samples of dataN. The higher the value N, 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).


### 4.3.3. Part III: Swarm Optimizing

Step 7.
Find the best positionpbest and the global position gbest 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 Ifk
g
=
G, output the results of parameter identification, and terminate the MPSO algorithm; otherwise, jump to Step 6.


## 4.3.1. Part I: The Particle Swarm Initialization

Step 1.
Set the appropriate population size Size, the maximum number of iterationG, and the value of local acceleration coefficients c
1 and global acceleration coefficients c
2.Step 2.
Determine the maximum and minimum value of inertia weighting factorw
max
⁡ and w
min
⁡.Step 3.
Set the position boundary value[
x
min
⁡
,
x
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 ind-dimensions searching space. Randomly initialize the position vector and velocity vector of each particle according to the following equations:
(14)
x
(
i
,
j
)
=
x
(
i
,
j
)
min
⁡
+
(
x
(
i
,
j
)
max
⁡
-
x
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
v
(
i
,
j
)
=
v
(
i
,
j
)
min
⁡
+
(
v
(
i
,
j
)
max
⁡
-
v
(
i
,
j
)
min
⁡
)
*
rand
(
)
,
p
p
p
p
p
p
p
i
=
[
1,2
,
…
,
Size
]
,
j
=
[
1,2
,
…
,
d
]
.


## 4.3.2. Part II: Fitness Evaluation

Step 5.
Acquire the output of PZT actuator and determine the number of samples of dataN. The higher the value N, 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).


## 4.3.3. Part III: Swarm Optimizing

Step 7.
Find the best positionpbest and the global position gbest 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 Ifk
g
=
G, output the results of parameter identification, and terminate the MPSO algorithm; otherwise, jump to Step 6.


## 5. Experimental Validation and Error Corrected Hysteresis Model

### 5.1. Experimental Design and MPSO Algorithm Parameters Selection

The experimental system is depicted in Section2. It considers that the input voltage of voltage amplifier is cosine voltage with the frequency by 0.1 Hz with the form of u
1
(
t
)
=
2
-
2
cos
⁡
⁡
(
0.2
π
t
) 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 N
=
10000. The criterion to stop the MPSO algorithm is that the maximum number of iterations is reached. The MPSO methods with different algorithmic parameters are compiled with MATLAB software. According to the experience of the previous studies, the position boundary values [
x
min
⁡
,
x
max
⁡
] of k
1, k
2, A
1, β
1, γ
1, and n are set to [
-
1,1
], [
-
10,10
], [
-
1,1
], [
-
10,10
], [
-
10,10
], and [
1,5
], respectively. Parameters of the MPSO algorithm are listed in Table 3.Table 3
List of MPSO algorithm parameters.
Parameters
Case1
Case2
Case3
Case4
θ
[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
]
d
6
x
min
⁡
[−1, −10, −1, −10, −10, 1]
x
max
⁡
[1, 10, 1, 10, 10, 5]
[
v
min
⁡, v
max
⁡
]
[−1, 1]
[
w
min
⁡, w
max
⁡
]
[0.10, 0.90]
Size
50
100
300
500
G
500
1500
5000
10000
c
1
1.3
c
2
1.7Equations (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.


### 5.2. 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:(15)
MAD
=
Max
⁡
1
≤
l
≤
N
|
y
^
(
l
)
-
y
(
l
)
|
,
MRD
=
MAD
y
Max
⁡
×
100
%
,
MAE
=
1
N
∑
l
=
1
N
|
y
^
(
l
)
-
y
(
l
)
|
,
NRMSD
=
1
y
Max
⁡
1
N
∑
l
=
1
N
[
y
^
(
l
)
-
y
(
l
)
]
2
,
where y
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 Table4. The hysteresis loops of voltage-displacement responses are shown in Figure 8 under the input of u
1
(
t
), within which, the simulation and experimental results are compared.Table 4
List of the identified parameters of hysteresis model.
k
1
k
2
A
1
β
1
γ
1
n
Case1
26.4355 * 10−6
2.2610
−15.4396 * 10−6
9.6240
−6.6160
1
Case2
28.3432 * 10−6
0.4260
−14.5375 * 10−6
5.3900
−5.3954
1
Case3
29.3383 * 10−6
0.2412
−15.3490 * 10−6
9.2620
−6.4164
1
Case4
29.3384 * 10−6
0.2410
−15.5320 * 10−6
9.2623
−6.0162
2The hysteresis loops of the PZT actuator obtained by experimental data and the response for the identified hysteresis model using MPSO algorithm under the cosine input voltageu
1
(
t
): (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4.
(a)
(b)
(c)
(d)Observing Figure8, the modelling error of Case 1 (Figure 8(a)) and Case 2 (Figure 8(b)) are significantly larger 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 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.


### 5.3. Error Correction of Hysteresis Model

Equation (8) can be manipulated and rewritten as follows:
(16)
h
˙
(
t
)
=
(
A
1
-
(
β
1
·
sign
⁡
(
u
˙
(
t
)
)
·
sign
⁡
(
h
(
t
)
)
+
γ
1
)
·
|
h
(
t
)
|
n
)
·
u
˙
(
t
)
.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:
(17)
y
(
t
)
=
k
1
u
(
t
)
+
k
2
h
(
t
)
,
h
˙
(
t
)
=
(
A
1
-
(
β
1
·
sign
⁡
(
V
˙
(
t
)
)
·
sign
⁡
(
h
(
t
)
)
+
γ
1
)
·
|
h
(
t
)
|
n
)
·
V
˙
(
t
)
+
Δ
Φ
,
Δ
Φ
=
ε
u
(
t
)
·
sign
⁡
(
u
˙
(
t
)
)
,
V
(
t
)
=
u
(
t
+
φ
)
,
where the ε is asymmetric coefficient. In the voltage rising process, the sign
⁡
(
u
˙
(
t
)
)
=
1, if ε
>
0; the higher the applied voltage u
(
t
), the more apparently the asymmetric factor Δ
Φ positive effects on h
(
t
); and the ε
<
0 causes the same trend of negative effects. In the voltage descending process, the sign
⁡
(
u
˙
(
t
)
)
=
-
1; these trends are exactly opposite. Moreover, the input bias φ can direct causes phase offset of output displacement y
(
t
). 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 A
1, β
1, γ
1, and n 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[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
] to the parameters identification results of Case 3 and Case 4. Secondly, obtain the optimal value of φ and ε by MPSO algorithm. Figures 9(a) and 9(b) show the comparison of experimental data and model predicted output of the two hysteresis models. The curve Exp, C model, and NC model represent the experiment curve, the corrected model curve, and noncorrected model curve, respectively. The output of the MPSO algorithm is φ
=
0.1047 and ε
=
1
.
35
×
1
0
-
8 for Case 3 and φ
=
-
0.0024 and ε
=
-
3
.
5
×
1
0
-
8 for Case 4. The trends of absolute error are shown in Figures 10(a) and 10(b), and quantitative error of the output displacements predicted is shown in Table 5.Table 5
Quantitative error of the output displacements predicted.
MAD
MRD
MAE
NRMSD
Non-corrected Case3
2.6782μm
2.2318%
1.2641μm
0.0131
Corrected Case3
0.6832μm
0.5693%
0.2340μm
0.0023
Non-corrected Case4
0.7726μm
0.6438%
0.3896μm
0.0037
Corrected Case4
0.5164μm
0.4303%
0.2200μm
0.0023The hysteresis loops of the PZT actuator obtained by experimental data and the model response for (17) using the parameters of Case 3 and Case 4: (a) the corrected Case 3 and (b) the corrected Case 4.
(a)
(b)The trends of absolute error. (a) The trends of absolute error for corrected Case 3 and noncorrected Case 3 and (b) the trends of absolute error for corrected Case 4 and noncorrected Case 4.
(a)
(b)Figures9(a) and 9(b) show that the fitting precision of the corrected hysteresis model is significantly better than the previous noncorrected model. According to Table 5, the maximum absolute modelling error of Figures 9(a) and 9(b) is 0.6832 um and 0.5164 um, that is, 0.5693% and 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 ofu
2
(
t
)
=
4
-
4
cos
⁡
⁡
(
0.2
π
t
), and use the MPSO to identify the parameters [
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
,
φ
,
ε
] for (17) with the Size
=
300 and G = 8000. The results of MPSO algorithm are k
1 = 39.3383*10−6, k
2 = 0.2412, A
1 = −15.3487*10−6, β
1 = 9.2620, γ
1 = −6.4160, n
=
2, φ
=
0.8612, and ε
=
1
.
05
×
1
0
-
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.Table 6
Quantitative error of the output displacements predicted.
MAD
MRD
MAE
NRMSD
Non-corrected model
1.3155μm
0.5481%
1.0396μm
4.35 × 10−5
Corrected model
0.8612μm
0.3588%
0.1835μm
8.29 × 10−6Figure 11
The hysteresis loops of the PZT actuator obtained by experimental data and the model response; it is using the MPSO to identify parameters withSize
=
300 and G
=
8000. The results of parameters identification are k
1 = 39.3383*10−6, k
2 = 0.2412, A
1 = −15.3487*10−6, β
1 = 9.2620, γ
1 = −6.4160, n
=
2, φ
=
0.8612, and ε
=
1
.
05
×
1
0
-
8.Figure 12
The time histories of the output displacements from the PZT actuator and the predictions of two models under the inputu
2
(
t
).Figure 13
The trends of absolute error of the two models under the inputu
2
(
t
).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 Figure14.Figure 14
The time histories of the output displacements from the PZT actuator and the response of the two models under the 0.1 Hz input voltage with decreasing amplitude.Observing Figures11, 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 the experimental results. Under the multiperiod input voltage, according the Table 6, relative to the noncorrected hysteresis model, the MAD, MAE, and NRMSD by the corrected hysteresis model are reduced by 34.5%, 82.3%, and 80.9%, respectively. Observing 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 rate-dependent hysteresis and constructing adaptive control system for the precise positioning control of the PZT actuator.


## 5.1. Experimental Design and MPSO Algorithm Parameters Selection

The experimental system is depicted in Section2. It considers that the input voltage of voltage amplifier is cosine voltage with the frequency by 0.1 Hz with the form of u
1
(
t
)
=
2
-
2
cos
⁡
⁡
(
0.2
π
t
) 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 N
=
10000. The criterion to stop the MPSO algorithm is that the maximum number of iterations is reached. The MPSO methods with different algorithmic parameters are compiled with MATLAB software. According to the experience of the previous studies, the position boundary values [
x
min
⁡
,
x
max
⁡
] of k
1, k
2, A
1, β
1, γ
1, and n are set to [
-
1,1
], [
-
10,10
], [
-
1,1
], [
-
10,10
], [
-
10,10
], and [
1,5
], respectively. Parameters of the MPSO algorithm are listed in Table 3.Table 3
List of MPSO algorithm parameters.
Parameters
Case1
Case2
Case3
Case4
θ
[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
]
d
6
x
min
⁡
[−1, −10, −1, −10, −10, 1]
x
max
⁡
[1, 10, 1, 10, 10, 5]
[
v
min
⁡, v
max
⁡
]
[−1, 1]
[
w
min
⁡, w
max
⁡
]
[0.10, 0.90]
Size
50
100
300
500
G
500
1500
5000
10000
c
1
1.3
c
2
1.7Equations (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.


## 5.2. 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:(15)
MAD
=
Max
⁡
1
≤
l
≤
N
|
y
^
(
l
)
-
y
(
l
)
|
,
MRD
=
MAD
y
Max
⁡
×
100
%
,
MAE
=
1
N
∑
l
=
1
N
|
y
^
(
l
)
-
y
(
l
)
|
,
NRMSD
=
1
y
Max
⁡
1
N
∑
l
=
1
N
[
y
^
(
l
)
-
y
(
l
)
]
2
,
where y
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 Table4. The hysteresis loops of voltage-displacement responses are shown in Figure 8 under the input of u
1
(
t
), within which, the simulation and experimental results are compared.Table 4
List of the identified parameters of hysteresis model.
k
1
k
2
A
1
β
1
γ
1
n
Case1
26.4355 * 10−6
2.2610
−15.4396 * 10−6
9.6240
−6.6160
1
Case2
28.3432 * 10−6
0.4260
−14.5375 * 10−6
5.3900
−5.3954
1
Case3
29.3383 * 10−6
0.2412
−15.3490 * 10−6
9.2620
−6.4164
1
Case4
29.3384 * 10−6
0.2410
−15.5320 * 10−6
9.2623
−6.0162
2The hysteresis loops of the PZT actuator obtained by experimental data and the response for the identified hysteresis model using MPSO algorithm under the cosine input voltageu
1
(
t
): (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4.
(a)
(b)
(c)
(d)Observing Figure8, the modelling error of Case 1 (Figure 8(a)) and Case 2 (Figure 8(b)) are significantly larger 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 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.


## 5.3. Error Correction of Hysteresis Model

Equation (8) can be manipulated and rewritten as follows:
(16)
h
˙
(
t
)
=
(
A
1
-
(
β
1
·
sign
⁡
(
u
˙
(
t
)
)
·
sign
⁡
(
h
(
t
)
)
+
γ
1
)
·
|
h
(
t
)
|
n
)
·
u
˙
(
t
)
.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:
(17)
y
(
t
)
=
k
1
u
(
t
)
+
k
2
h
(
t
)
,
h
˙
(
t
)
=
(
A
1
-
(
β
1
·
sign
⁡
(
V
˙
(
t
)
)
·
sign
⁡
(
h
(
t
)
)
+
γ
1
)
·
|
h
(
t
)
|
n
)
·
V
˙
(
t
)
+
Δ
Φ
,
Δ
Φ
=
ε
u
(
t
)
·
sign
⁡
(
u
˙
(
t
)
)
,
V
(
t
)
=
u
(
t
+
φ
)
,
where the ε is asymmetric coefficient. In the voltage rising process, the sign
⁡
(
u
˙
(
t
)
)
=
1, if ε
>
0; the higher the applied voltage u
(
t
), the more apparently the asymmetric factor Δ
Φ positive effects on h
(
t
); and the ε
<
0 causes the same trend of negative effects. In the voltage descending process, the sign
⁡
(
u
˙
(
t
)
)
=
-
1; these trends are exactly opposite. Moreover, the input bias φ can direct causes phase offset of output displacement y
(
t
). 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 A
1, β
1, γ
1, and n 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[
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
] to the parameters identification results of Case 3 and Case 4. Secondly, obtain the optimal value of φ and ε by MPSO algorithm. Figures 9(a) and 9(b) show the comparison of experimental data and model predicted output of the two hysteresis models. The curve Exp, C model, and NC model represent the experiment curve, the corrected model curve, and noncorrected model curve, respectively. The output of the MPSO algorithm is φ
=
0.1047 and ε
=
1
.
35
×
1
0
-
8 for Case 3 and φ
=
-
0.0024 and ε
=
-
3
.
5
×
1
0
-
8 for Case 4. The trends of absolute error are shown in Figures 10(a) and 10(b), and quantitative error of the output displacements predicted is shown in Table 5.Table 5
Quantitative error of the output displacements predicted.
MAD
MRD
MAE
NRMSD
Non-corrected Case3
2.6782μm
2.2318%
1.2641μm
0.0131
Corrected Case3
0.6832μm
0.5693%
0.2340μm
0.0023
Non-corrected Case4
0.7726μm
0.6438%
0.3896μm
0.0037
Corrected Case4
0.5164μm
0.4303%
0.2200μm
0.0023The hysteresis loops of the PZT actuator obtained by experimental data and the model response for (17) using the parameters of Case 3 and Case 4: (a) the corrected Case 3 and (b) the corrected Case 4.
(a)
(b)The trends of absolute error. (a) The trends of absolute error for corrected Case 3 and noncorrected Case 3 and (b) the trends of absolute error for corrected Case 4 and noncorrected Case 4.
(a)
(b)Figures9(a) and 9(b) show that the fitting precision of the corrected hysteresis model is significantly better than the previous noncorrected model. According to Table 5, the maximum absolute modelling error of Figures 9(a) and 9(b) is 0.6832 um and 0.5164 um, that is, 0.5693% and 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 ofu
2
(
t
)
=
4
-
4
cos
⁡
⁡
(
0.2
π
t
), and use the MPSO to identify the parameters [
k
1
,
k
2
,
A
1
,
β
1
,
γ
1
,
n
,
φ
,
ε
] for (17) with the Size
=
300 and G = 8000. The results of MPSO algorithm are k
1 = 39.3383*10−6, k
2 = 0.2412, A
1 = −15.3487*10−6, β
1 = 9.2620, γ
1 = −6.4160, n
=
2, φ
=
0.8612, and ε
=
1
.
05
×
1
0
-
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.Table 6
Quantitative error of the output displacements predicted.
MAD
MRD
MAE
NRMSD
Non-corrected model
1.3155μm
0.5481%
1.0396μm
4.35 × 10−5
Corrected model
0.8612μm
0.3588%
0.1835μm
8.29 × 10−6Figure 11
The hysteresis loops of the PZT actuator obtained by experimental data and the model response; it is using the MPSO to identify parameters withSize
=
300 and G
=
8000. The results of parameters identification are k
1 = 39.3383*10−6, k
2 = 0.2412, A
1 = −15.3487*10−6, β
1 = 9.2620, γ
1 = −6.4160, n
=
2, φ
=
0.8612, and ε
=
1
.
05
×
1
0
-
8.Figure 12
The time histories of the output displacements from the PZT actuator and the predictions of two models under the inputu
2
(
t
).Figure 13
The trends of absolute error of the two models under the inputu
2
(
t
).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 Figure14.Figure 14
The time histories of the output displacements from the PZT actuator and the response of the two models under the 0.1 Hz input voltage with decreasing amplitude.Observing Figures11, 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 the experimental results. Under the multiperiod input voltage, according the Table 6, relative to the noncorrected hysteresis model, the MAD, MAE, and NRMSD by the corrected hysteresis model are reduced by 34.5%, 82.3%, and 80.9%, respectively. Observing 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 rate-dependent hysteresis and constructing adaptive control system for the precise positioning control of the PZT actuator.


## 6. 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.



---

*Source: 932974-2014-02-13.xml*