**1. Introduction**

With the growth of semiconductor and precision manufacturing industries, the positioning accuracy of micro/nano scale is highly required [1–4]. Because of the advantages of fast frequency response, high stiffness and displacement resolution, piezoelectric actuators (PEAs) based on inverse piezoelectric effect have gradually become one of the most widely used smart material actuators. In addition, PEAs are popularly applied as the actuators in micro/nano scale measurement, micro-electro-mechanical systems (MEMS), flexible electronics manufacturing, and biomedical engineering. However, strong nonlinear hysteresis of piezoelectric materials makes the output of PEAs difficult to predict, and the positioning accuracy and stability of the system are low [5,6]. Specifically, scholars provided plenty of mathematical models to describe the symmetric hysteresis. However, effective mathematical model to describe the asymmetric hysteresis is still lacking, as the asymmetric hysteresis is a more common hysteresis nonlinear phenomenon. Therefore, it is worth exploring the modeling of asymmetric hysteresis behavior and using it to compensate piezoelectric actuators.

To eliminate the influence of hysteresis on accuracy and stability of the system, scholars have proposed numerous hysteresis models and controllers. These hysteresis models can be divided into two categories [7]: mechanistic model and phenomenological model. The former are based on basic physical principles, and derived by energy-displacement or stress–strain methods. Among them, the famous models are Jiles–Atherton model [8,9] and Domain Wall model [10]. The latter uses directly mathematical means to characterize the hysteresis, and ignore the physical meaning behind the hysteresis phenomenon. Because of the high accuracy and flexibility, the phenomenological model is more popular in hysteresis modeling. The phenomenological models include Preisach model [11,12], polynomial model [13,14], Bouc–Wen model [15,16], Duhem model [17,18], neural network model [19,20], Prandtl–Ishlinskii (PI) model [21–23], and etc. Among them, because of simple expression and analytical inverse model, the PI model is the most widely used in hysteresis modeling and compensation. However, the PI model utilizes weighted superposition of the Play operators to describe hysteresis nonlinearity; the Play operator is a basic component of the PI model. Therefore, the symmetric structure of the Play operator determines that PI model can only describe symmetric hysteresis.

Actually, PEAs have slight or severe asymmetric characteristics, as shown in Figure 1. When the application scope of PEAs scales down to micro/nano-meter levels, the gap between symmetric and asymmetric hysteresis modeling approaches results in positioning error. To describe asymmetric hysteresis, scholars [24–26] have tried to modify the PI model, and designed the corresponding controller. Kuhnen [24] presented the dead-zone operator and the modified PI (Ku-PI) model to describe the asymmetric hysteresis. Janaideh et al. [25] provided the generalized PI (GPI) model based on a nonlinear play operator. Di fferent envelope functions make the GPI model have flexible ascending and descending edge. Hence, the GPI model can describe accurately the complex hysteresis phenomenon. By combining memoryless polynomial with PI model, Gu [26] proposed a new modified PI (Gu-PI) model to depict asymmetric hysteresis. These modified PI models and corresponding controller play a good role in the micro/nano positioning compensation. However, these models still have some limitations. The parameter identification of Ku-PI model is complex, which increases the di fficulty of compensator design; GPI model is flexible and accurate, but choosing envelope function still depends on the experience and recursive debugging; Gu-PI model is simple in structure and easy to identify, and has good performance with one-side Play operator when the initial loading curve is not considered. But we have tried to extend the application of Gu-PI model with two-side Play operator when the initial loading curve is considered, and its e ffect is poor. The reason is that the Play operator lacks accuracy in describing displacements near the zero voltage.

**Figure 1.** Asymmetric hysteresis loops of piezoelectric actuators (PEAs).

In this paper, a new polynomial-modified PI (PMPI) model is proposed to describe the asymmetric hysteresis of PEAs. Compared with PI model, the most important innovation of PMPI model is the introduction of Modified-Play (M-Play) operator and memoryless polynomial. M-Play operator replaces Play operator as the basic operator in the PMPI model. The memoryless polynomial enables PMPI model to describe the asymmetric hysteresis. The shape of asymmetric hysteresis is determined

by both weighted M-Play operators and memoryless polynomial. The main advantages of PMPI model are as following: (1) Whether the initial loading curve is considered or not, in PMPI model, the displacement near zero voltage can be described flexibly by M-Play operator. (2) The inverse of PMPI model can be derived directly from PI model. The feedforward compensation of hysteresis in real-time application can be carried out. To validate the proposed PMPI model and I-M compensator, simulation and experiment are conducted on a piezoelectric micromotion platform.

This paper is organized as follows: Section 2 introduces PMPI model and examines the congruency property. Section 3 develops an I-M compensator and analyzes its stability. Section 4 verifies the PMPI model and I-M compensator on a piezoelectric micromotion platform. Section 5 concludes this paper.

#### **2. Polynomial-Modified Prandtl–Ishlinskii Model**

Before introducing the proposed PMPI model, it is necessary to review the PI model in brief.

## *2.1. Prandtl–Ishlinskii Model*

The PI model utilizes weighted superposition of the Play operators to describe hysteresis nonlinearity, so the Play operator is the basic component of PI model. Without special description in the following paper, the Play operator refers to two-side Play operator. For any piecewise monotonic input signal *v*(*t*) ∈ *Cm*(0 *tN*), the output *w*(*t*) = *Fr*(*v*)(*t*) of the Play operator with threshold *r* is defined as

$$\begin{aligned} w(0) &= F\_r(v)(0) = f\_r(v(0), 0) \\ w(t) &= F\_r(v)(t) = f\_r(v(t), w(t\_{i-1})) \end{aligned} \tag{1}$$

for *ti*-1 < *t*< *ti*, 1 < *i* < *N*, with

$$f\_r(v, w) = \max(v - r, \min(v + r, w))\tag{2}$$

where 0<*t*1<*t*2<···<sup>&</sup>lt;*tN* is a division of the time domain (0 *tN*) to ensure that the input signal *v*(*t*) is monotonic within each subinterval (*ti*-1 *ti*). The Play operator have partial memory and rate-independent properties. That is, the output of the Play operator depends not only on current input, but also on the input history. However, the input rate does not change the output shape. As an illustration, Figure 2 shows the input–output relationship of play operator. It is easy to notice that there is a symmetry center in the input–output trajectory.

**Figure 2.** The input–output relationship of the Play operator.

The PI model utilizes the Play operator *Fr*(*v*)(*t*) to describe the hysteresis relationship between input *v* and output Γ*PI*:

$$\begin{cases} \Gamma\_{Pl}(\upsilon)(k) = p\_0 \upsilon(k) + \sum\_{i=1}^n p(r\_i) F\_{r\_i}(\upsilon)(k) \\ \text{s.t. } p(r\_i) \ge 0 \end{cases} \tag{3}$$

where *n* is the number of Play operators, *p*0 is a positive constant. *ri* is the threshold of the Play operator, *p*(*ri*) represents the weighted coe fficient for the threshold *ri* and approaches 0 as *ri* becomes larger. If *p*(*ri*) < 0, the PI model will fail to correctly describe hysteresis minor loops.

Observed from Equation (3), the PI model is composed of the weighted superposition of Play operator *Fr*(*v*)(*t*) and linear input signal. The input–output curve is parallelogram, which is in central symmetry, so the PI model can only be used to characterize symmetric hysteresis. As an illustration, Figure 3 shows the hysteresis loop generated by the PI model with *p*0 = 2.1, *p*(*ri*) = 2.6*e*<sup>−</sup>0.25*ri*, *ri* = 0:0.5:9.5, and input *v*(*t*) = 5 + 5sin(2 π*t* − π/2). It is obvious that the hysteresis curve generated by PI model is symmetrical. But PEAs often exhibit more or less asymmetric characteristics as shown in Figure 1. In addition, the max displacement of zero voltage in the PI model is uniquely determined by the initial loading curve, and the relationship can be expressed by Equation (4). In fact, the max displacement of the zero voltage is not necessarily related to the initial loading curve, which reflects that the Play operator lacks accuracy and flexibility in describing displacement near zero voltage. One of the motivations of the paper is to modify the PI model to flexibly and accurately represent the asymmetric hysteresis of PEAs.

 *n*

*p*(*ri*)*ri* (4)

*i*=1

*y*(0) =

**Figure 3.** Hysteresis loops generated by the Prandtl–Ishlinskii (PI) model.

#### *2.2. Polynomial-Modified Prandtl–Ishlinskii Model*

Utilizing PI model to describe the asymmetric hysteresis of PEAs will produce considerable errors, which cannot meet the needs of precise positioning. To characterize accurately the asymmetric hysteresis, many modified PI models have been proposed, such as Ku-PI model, Gu-PI model, GPI model, and etc. Although these modified PI models can describe the asymmetric hysteresis to some extent, they still have some limitations. To improve the flexibility and accuracy of the model in describing asymmetric hysteresis, we propose Modified-Play (M-Play) operator and polynomial-modified PI model (PMPI).

The M-Play operator is derived from the Play operator, and it is formed by multiplying the threshold value of Play operator on the descending edge by a delay coe fficient η > −1. The M-Play operator can be written as

$$\begin{cases} w(0) = F\_{r, \eta\_i}(v)(0) = f\_{r, \eta\_i}(v(0), 0) \\ w(t) = F\_{r, \eta\_i}(v)(t) = f\_{r, \eta\_i}(v(0), w(t\_i)) \end{cases} \tag{5}$$

where

$$f\_{r, \eta\_i}(v, w) = \max(v - r, \min(v + \eta\_i r, w))\tag{6}$$

The coe fficient η alters the threshold of the descending edge. The larger the η value, the more obvious the delay in the descending state. With proper values of η for every individual M-Play operator, the flexibility and accuracy of the model can be significantly enhanced.

Figure 4 shows the response of M-Play operators with di fferent delay coe fficient η. Obviously, one-side Play operator and two-side Play operator are special cases of M-Play operator. When η = 0, M-Play operator is equivalent to one-side Play operator, and when η = 1, M-Play operator is equivalent to two-side Play operator.

**Figure 4.** The response of M-Play operators with di fferent delay coe fficient η*.*

Remarks: From Figure 3, although the delay coe fficient η is introduced, it can be seen from the input and output trajectories of M-Play operator still exist in the symmetry center. The weighted superposition of M-Play operators alone cannot characterize the asymmetric hysteresis, and the delay coe fficient η is only used to improve the description of the displacement near zero voltage.

In this paper, the proposed PMPI model is formulated as:

⎪⎪⎪⎪⎪⎪⎪⎪⎨

$$\begin{cases} H(\boldsymbol{\upsilon})(k) = \Gamma p\_{l,\boldsymbol{\eta}}(\boldsymbol{\upsilon})(k) + P(\boldsymbol{\upsilon})(k) \\ \Gamma p\_l(\boldsymbol{\upsilon})(k) = p\_0 \boldsymbol{\upsilon}(k) + \sum\_{i=1}^n p\_i F\_{r\_i \boldsymbol{\eta}\_i}(\boldsymbol{\upsilon})(k) \\ P(\boldsymbol{\upsilon})(k) = a\_1 \boldsymbol{\upsilon}(k)^3 + a\_2 \boldsymbol{\upsilon}(k)^2 + a\_3 \\ \text{s.t. } p\_i \ge 0, \ \boldsymbol{\eta}\_i \ge -1 \end{cases} \tag{7}$$

where *p*0 and *p*(*ri*) are defined the same as the ones in PI model (3), and *Fr,*η*<sup>i</sup>* are the M-Play operator. A third-degree memoryless polynomial *P*(*v*)(*t*) is used to characterize asymmetric hysteresis of PEAs.

$$y(0) = \sum\_{i=1}^{n} p(r\_i)\eta\_i r\_i + a\_3 \tag{8}$$

The proposed PMPI consists of two parts: several weighted M-Play operators (denoted as MPI model) and memoryless polynomial *P*(*v*)(*t*). The introduction of M-Play operator enables the PMPI model to describe accurately the displacement near zero voltage. It can be seen from Equation (8) that, in PMPI model, the max displacement of zero voltage is adjusted flexibly by the delay coe fficient η. The combination of M-Play operators and third-order memoryless polynomial can characterize the asymmetric hysteresis of PEAs. To illustrate the advantages of PMPI model, Figure 5 shows the hysteresis loops generated by PMPI model (7) with *P*(*v*)(*t*) = −0.05*v*<sup>3</sup> + 1.2*v*<sup>2</sup> + 0.2 and η = 0:0.025:0.475. It is worth noting that weight function *p*(*r*), positive constant *p*0, and input *v*(*t*) are the same as of the ones used in the hysteresis loops in Figure 3. In contrast, it is the M-Play operator and memoryless polynomial that enabled the PMPI model to characterize flexibly and accurately the serious asymmetric hysteresis.

**Figure 5.** Hysteresis loops generated by polynomial-modified PI (PMPI) model.

#### *2.3. Congruency Property*

The congruency property is one of the basic properties of hysteresis model. The congruency property means that minor hysteresis loops with different input history are congruen<sup>t</sup> in the same input range. In this section, we prove the congruency property of PMPI model and establish the minor loop mathematical model.

Because of the introduction of M-Play operator and memoryless polynomial *P*(*v*)(*t*), the congruency property of PMPI model needs to be proven. Memoryless polynomial *P*(*v*)(*t*) is a bijective function and has no partial memory, its output depends only on the current input. Therefore, the congruency property of PMPI model depends on the MPI model.

To illustrate the congruency property of MPI model, Figure 6 shows the hysteresis loops of M-Play operator with different input history in the same input range. Figure 6a,c shows two input signals with different input history in the same range, and Figure 6b,d shows the corresponding hysteresis curves. Combining these four graphs, it can be clearly seen that, although the input history of the two signals is different, the shapes of the corresponding minor loop formed by M-Play operators are same in the same input range (*vm vM*). This example illustrates the congruency property of the M-Play operator, and further demonstrates that MPI model has congruency property.

Sprekels's monograph [27] has proved that the shape of minor loop is uniquely determined when the vertical height *h* of minor loop are constant. Therefore, we quantitatively represent the shape of minor loop with its vertical height *h*. Figure 7 shows the geometric relations of the minor loop of M-Play operator. It can be seen from the figure that the relationship between vertical height *hr* and interval length *x* can be formulated as:

$$h\_r(\mathbf{x}) = \max(\mathbf{x} - (1 + \eta)r, 0) \tag{9}$$

where *x* is the interval length of input range (*vm vM*). When *x* < (1 + η)*<sup>r</sup>*, the M-Play operator cannot form minor loop, that is, *hr* = 0. Based on Equation (8), the vertical height *hm* of the minor loop of MPI model can be expressed as:

$$h\_{\mathfrak{m}}(\mathbf{x}) = p\_0 \mathbf{x} + \sum\_{i=1}^{n} p\_i h\_{r\_i}(\mathbf{x}) \tag{10}$$

As mentioned above, the output of memoryless polynomial function *P*(*v*)(*t*) only depends on the current input, and there is no partial memory. The PMPI model is based on the MPI model plus *P*(*v*)(*t*), so the vertical height *hpm* of the minor loop of PMPI model can be expressed as:

$$h\_{\rm pwl}(\mathbf{x}) = p\mathbf{u}\mathbf{x} + \sum\_{i=1}^{n} p\_i h\_{\Gamma\_i}(\mathbf{x}) + \left[P(\mathbf{v}\mathbf{u}) - P(\mathbf{v}\_{\rm m})\right] \tag{11}$$

**Figure 6.** Hysteresis loops of M-Play operator with different input history in the same input range. (**a**) Input history 1; (**b**) hysteresis curve generated by input history 1; (**c**) input history 2; (**d**) hysteresis curve generated by input history 2

**Figure 7.** Geometric relations of the minor loop of M-Play operators.

For example, the input *v*(*t*) is given to verify the correctness of Equation (11). In this case, the sequence of input maxima and minima is defined as {0→7→4→10→4→7→0}. Thus, the input voltage has the same input range (4 7). The parameters of the PMPI model are chosen as *r* = 0:0.5:9.5, *p*0 = 3.93, *pi* = 2.5*e*<sup>−</sup>0.44*ri*, η = 0:0.05:0.95, *P*(*v*)(*t*) = −0.02*v*<sup>3</sup> − 0.48*v*2. For the input voltage shown in Figure 8a, the associated hysteresis loops of PMPI model are given in Figure 8b. It can be seen from the partial enlarged drawing of Figure 8b that the shapes of the minor loops are identical in the same input range and the height of the minor loops are both 19.02 μm. According to the Equation (11) derived above, the vertical height of the minor loops also are 19.02 μm. This example verifies that the PMPI model also has congruency property and the minor loop mathematical model is correct.

(**a**)

**Figure 8.** Simulation of the congruency property. (**a**) Input voltage. (**b**) Hysteresis loops.

#### **3. The Design and Analysis of the Inverse Model Compensator**

Feedforward compensation based on inverse hysteresis model is efficient and practical for reducing hysteresis effect in open-loop systems. The main idea of this method is to construct inverse hysteresis model, which is cascaded to the controlled object for feedforward control. After the cascade, the piezoelectric system can be approximated as a linear system, as shown in Figure 9. In this section, based on PMPI model, we will design the inverse model (I-M) compensator to compensate hysteresis and analyze its stability.

**Figure 9.** Schematic illustration of feedforward compensation based on inverse hysteresis model.

#### *3.1. Inverse model Compensator Design*

The MPI model introduces delay coefficient η on the basis of PI model, but this has no influence on the inverse MPI model, that is, the parameter expression of inverse MPI model is the same as that of inverse PI model. This is because that the parameter expression of inverse PI model is derived by the initial loading curve, the delay coefficient η only has influence on the descending edge but not the initial loading curve. The initial loading curve is still expressed as:

$$p(r) = p\_0 r + \int\_0^r p(\zeta)(r-\zeta)d\zeta. \tag{12}$$

The MPI model has an analytical inverse model, but the PMPI model has one more memoryless polynomial *P*[*u*](*k*). Hence, there is no analytical inverse model. The inverse PMPI model can only be obtained by iterating the memoryless part. In this section we utilize the iterative structure to design a compensator as shown in Figure 10. In practical applications, most compensator systems are discrete, when the sampling frequency is sufficiently high, *ui* ≈ *ui*−1. So the I-M compensator can be expressed as:

$$u\_d(k) = H^{-1}[y\_d](k) = \Gamma\_{MPI}^{-1}(y\_d[k] - P[u\_d](k-1))\tag{13}$$

where *yd*(*k*) is the desired displacement, and *ud*(*k*) is the desired control voltage. Γ−1 *MPI*[·](*k*) stands for the inverse MPI model and can be calculated according to the following formula:

$$
\Gamma\_{MPl}^{-1}[y\_d](k) = \hat{p}\_0 y\_d(k) + \sum\_{i=1}^n \hat{p}\_i F\_{\mathbb{H}\_i \eta\_i}[y\_d](k) \tag{14}
$$

with

$$\begin{aligned} \hat{p}\_i &= -\frac{\frac{1}{p\_0} \text{ s.t. } p\_0 > 0\\ \hat{p}\_i &= -\frac{\frac{i}{p\_0} \text{ s.t. } p\_i}{(p\_0 + \sum\_{j=1}^i p\_j)(p\_0 + \sum\_{j=1}^{i-1} p\_j)}\\ \hat{r}\_i &= p\_0 r\_i + \sum\_{j=1}^{i-1} p\_j (r\_i - r\_j) \end{aligned} \tag{15}$$

**Figure 10.** Schematic diagram of inverse model (I-M) compensator.

## *3.2. Stability Analysis*

It can be seen from Figure 10 and Equation (11) that I-M compensator obtains the desired control voltage *ud*(*k*) by iterative solution. The divergence problem may occur during the iterative solution process. Therefore, it is necessary to analyze the stability of the I-M compensator. In this paper, the small gain theorem is used to analyze the stability. The small gain theorem states that when the

maximum closed-loop gain of a closed-loop system satisfies |K|<1, then the system is stable in any case. The stable condition can be described as following:

$$\begin{split} \max \{ K\_{inv}(k) | K\_{inv}(k) = \left| \frac{d \Gamma\_{MP}^{-1}[y\_d](k)}{dy\_d} \cdot \frac{dP[u\_d](k)}{du\_d} \right| \} < 1 \\ \Rightarrow p\_0 > \max \{ 3a\_1 u\_d(k)^2 + 2a\_2 u\_d(k) \} \end{split} \tag{16}$$

where *Kinv*[*u*](*k*) denotes the absolute gain of the I-M compensator, and *p*0 *a*1 and *a*2 are the parameters of the PMPI model.

Proof: To satisfy Equation (16), the maximum value of the product of absolute gain |*dP*[*ud*](*k*)/*dud*(*k*)| and |*d*(Γ− 1*MPI*[*yd*](*k*)/*d*(*yd*(*k*)))| should be less than 1. From Equation (6), we can ge<sup>t</sup>

$$0 \le d(F\_{\mathbb{P}\_i, \mathbb{P}\_l}[y\_d](k) / d(y\_d(k))) \le 1 \tag{17}$$

Combining Equations (15), the following relationship can be derived

$$\begin{aligned} \sum\_{i=1}^{n} \hat{p}\_i &= \sum\_{i=1}^{n} \left\{ \frac{1}{p\_0 + \sum\_{j=1}^{i} p\_j} - \frac{1}{p\_0 + \sum\_{j=1}^{i-1} p\_j} \right\} \\ &= \frac{1}{p\_0 + \sum\_{i=1}^{n} p\_i} - \frac{1}{p\_0} = \frac{1}{p\_0 + \sum\_{i=1}^{n} p\_i} - \hat{p}\_0 \Rightarrow -\hat{p}\_0 < \sum\_{i=1}^{n} \hat{p}\_i < 0 \end{aligned} \tag{18}$$

Therefore, the upper limit of the gain *d*(Γ− 1 *MPI*[*yd*](*k*)/*d*(*yd*(*k*))) is 1/*p*0. Since the memoryless polynomial *P*[*u*](*k*) is di fferentiable, its gain can be expressed by the following formula:

$$\left|dP[\mu\_d](k)/d\mu\_d\right| = \left|3a\_1\mu\_d(k)^2 + 2a\_2\mu\_d(k)\right|\tag{19}$$

Combining Equations (18) and (19), it can be concluded that when the identified parameters of PMPI model meet the condition (16), I-M compensator is globally stable.

The I-M compensator using iterative structure has superior performance in accuracy and response speed, but it may appear that the parameter identified by PMPI model does not satisfy Equation (16). Once this happens, the proportional gain *keud*(*k*) addressed in Equation (20) can be introduced to adjust the ratio between the MPI and memoryless polynomial part. In addition, the proportional gain *keud*(*k*) can greatly improve the convergence speed of I-M compensator. Compared with the dichotomy of pure iterative, the I-M compensator utilizes the characteristic of the analytical inverse of the MPI model, and its iterative process is approximately open-loop. Hence, the I-M compensator has fewer iterative steps, faster convergence speed, and higher accuracy.

$$\begin{aligned} H[u\_d](k) &= \left(\Gamma\_{MPl}[u\_d](k) + k\_t u\_d(k)\right) + \left(P[u\_d](k) - k\_t u\_d(k)\right) \\ &\Rightarrow p\_0 + k\_t > \max\left|3a\_1 u(k)^2 + 2a\_2 u(k) - k\_t\right| \end{aligned} \tag{20}$$

#### **4. Experimental Verification and Discussion**

In this section, an experimental platform is established. The experiment is conducted to verify the effectiveness of the PMPI model and I-M compensator in hysteresis modeling and compensation.
