A Digitized Representation of the Modified Prandtl–Ishlinskii Hysteresis Model for Modeling and Compensating Piezoelectric Actuator Hysteresis

Abstract

Piezoelectric actuators are widely used in micromanipulation and miniature robots due to their rapid response and high repeatability. The piezoelectric actuators often have undesired hysteresis. The Prandtl–Ishlinskii (PI) hysteresis model is one of the most popular models for modeling and compensating the hysteresis behaviour. This paper presents an alternative digitized representation of the modified Prandtl–Ishlinskii with the dead-zone operators (MPI) hysteresis model to describe the asymmetric hysteresis behavior of piezoelectric actuators. Using a binary number with n digits to represent the classical Prandtl–Ishlinskii hysteresis model with n elementary operators, the inverse model can be easily constructed. A similar representation of the dead-zone operators is also described. With the proposed digitized representation, the model is more intuitive and the inversion calculation is avoided. An experiment with a piezoelectric stacked linear actuator is conducted to validate the proposed digitized MPI hysteresis model and it is shown that it has almost the same performance as compared to the classical representation.

1. Introduction

Piezoelectric actuators are widely used in micromanipulation applications due to their rapid response and high repeatability, for example, microsurgical robots [1,2], autofocus optical systems [3,4,5], precise fabrication [6,7,8] and other applications [9,10,11]. One of the biggest challenges while using the piezoelectric actuators in dynamic applications is to model and compensate for the undesired complex hysteresis.

Various methods have been proposed for modeling and compensating the hysteresis behavior. The existing methods can be classified into physics-based models and phenomenology-based models [12]. Physics-based models [13] are often derived on the basis of physical principles of certain material or system properties [14,15,16]. Physics-based models often require a deep understanding of the causes of hysteresis and are often specific to the related properties.

The phenomenology-based model can also be roughly divided into three main groups [17]: differential-based models, such as the Bouc-wen [18] model, the Duhem [19] model, the Dahl, and the LuGre [20] model; neural-network models, such as the back propagation neural network based model [21], the gated recurrent unit based model [22], the neural network adaptive control method [23] and etc; operator-based models, such as the Preisach [24] model, the Krasnosel’skii–Pokrovskii(KP) [25] model, the Maxwell-slip [26] model, the Prandtl–Ishlinskii(PI) [27] model, and their variations.

The PI hysteresis model is one of the most popular models because of its simple implementation and analytical inversion [28]. Ang [29] concluded the relationship between the rate and the operators and modified the PI hysteresis model to model the rate-dependent hysteresis behavior. Tan [30] extended the PI hysteresis model to compensate for the ill-conditioned hysteresis behavior with a negative gradient. Mohammad [31] proposed a nonlinear play operator for modeling the asymmetric hysteresis of the Shape Memory Alloy (SMA) operators. Kuhnen [32] proposed to use the dead-zone operators for modeling the memory-free asymmetric hysteresis behavior.

In the above literature, the design procedure generally consists of parameters identification of the description model and the construction of the inverse model as a desired compensator. The mathematical complexity of the identification and inversion problem depends on the selected modeling approaches [33]. The amount of calculation may increase significantly with an increase in the number of elementary operators.

In this paper, we present an alternative digitized representation with which the inverse model can be easily constructed to compensate for the hysteresis behavior. The inversion calculation is avoided. The proposed representation is also more intuitive due to the geometric meanings of the parameters of the operators.

The paper is organized as follows: Section 2 describes the classical representation from which the proposed representation can be derived. Section 3 gives the experimental results with discussion. Section 4 covers the conclusion.

2. Hysteresis Mathematical Model

This section describes the classical representation and an alternative digitized representation of the classical PI hysteresis model. A similar representation of the dead-zone operators is also described, the construction of the inverse model with the proposed representation is also presented.

2.1. Classical PI Hysteresis Model

The classical PI hysteresis model has a series of elementary operators, the system output can be obtained by a weighted sum of n elementary operators with different thresholds. As shown in Figure 1, the working principles of the classical representation of the PI hysteresis model can be represented using the following equations:

(1)$y\left(0\right)=y_0$

(2)$y\left(t\right)=\sum _{i=1}^n{w}_i·max\{u\left(t\right)-r_i,min\{u\left(t\right)+r_i,y_{ei}(t-T)\}\}$

where $w_i$ is the weight corresponding to the ith elementary operator, $r_i$ is the threshold of the ith elementary operator, n is the number of the elementary operators and $y_{ei}$ is the output of the ith elementary operator.

When the parameter $r_i$ and $w_i$ are identified, the hysteresis curve model in the $(u,y)$ plane is shifted so that it is in the positive section of the plane. An example with three elementary operators is shown in Figure 2, $b{w}_i=2r_i$ is the input value $x_i$ of turning point where the slope of the hysteresis curve is changed. The $k_{th}$ output value $y_k$ can be formulated as follows [34]:

(3)$y_k=\sum _{i=1}^k(b{w}_{k+1}-b{w}_i)·w_i$

The Figure 2 can also be regarded as a series of 4 segments. Let $(\mathsf{\Delta }{x}_i,\mathsf{\Delta }{y}_i)$ by the input and output value of each of the n segments, the value of the $\mathsf{\Delta }{x}_i$ is the change of $b{w}_i$ while the value of the $\mathsf{\Delta }{y}_i$ is the change of $y_k$. The slope $gra{d}_k$ of the kth segment can be calculated using the following equation:

(4)$gra{d}_1=0;$

(5)$gra{d}_k=\sum _{i=1}^k{w}_i;$

As $b{w}_i=2r_i$, the value of $\mathsf{\Delta }{x}_i$ can be calculated using the following equations:

(6)$\mathsf{\Delta }{x}_1=2·r1;$

(7)$\mathsf{\Delta }{x}_i=2·(r_{i+1}-r_i);$

(8)$\mathsf{\Delta }{x}_n=2·(u_A-r_{n-1});$

The value of $\mathsf{\Delta }{y}_k$ can then be obtained using the following equations:

(9)$\mathsf{\Delta }{y}_1=0;$

(10)$\mathsf{\Delta }{y}_k=(\sum _{i=1}^{i}w(i))·\mathsf{\Delta }{x}_k$

The paired property $(\mathsf{\Delta }{x}_i,\mathsf{\Delta }{y}_i)$ represents the change of input and output value of each segment of the hysteresis curve and is thus more intuitive and the parameters are easier to be identified. The value of the paired property can be derived from the parameters of the classical representation. With this paired property, an alternative digitized representation of the classical PI hysteresis model is described in the next section.

2.2. Digitized Classical PI Hysteresis Model

The state of the classical PI hysteresis model with n operators is represented with a binary number $B=b_n\dots b_2{b}_1$ with n digits. Each digit of the binary number is a elementary operator with a binary state that can be either 0 or 1. Each elementary operator $b_i$ also has a paired property $(\mathsf{\Delta }{x}_i,\mathsf{\Delta }{y}_i)$, where $\mathsf{\Delta }{x}_i$ represents a change in the system input and $\mathsf{\Delta }{y}_i$ represents a change in the system output. An example with two operators is shown in Figure 3. With an increase in the number of operators, the system output and input can be approximated to be a linear relationship within each operator. The paired properties H of a PI hysteresis model with n elementary operators can be expressed in one equation:

(11)$H=\left[\begin{array}{ccc}\mathsf{\Delta }{x}_1& \cdots & \mathsf{\Delta }{x}_n\\ \mathsf{\Delta }{y}_1& \cdots & \mathsf{\Delta }{y}_n\end{array}\right]$

The model also has a paired property $(X_B,Y_B)$ that corresponds to its state. The state of the system always changes to either one of its two neighbors: upper state or lower state. Let $b_1$ be the least significant digit and $b_n$ be the most significant digit, the change of the state is always from the least significant digit. Let the least significant digit with a binary state of 0 at the current state be the A digit, let the least significant digit with a binary state of 1 at the current state be the Z digit.

When the state of the system changes to its upper neighbor, the binary state of the $b_A$ digit is changed from 0 to 1, the paired property of the system is increased by $(\mathsf{\Delta }{x}_A,\mathsf{\Delta }{y}_A)$. Denote $X_{B_t}$, $Y_{B_t}$ be the system input and output of the current state and $X_{B_{t+1}}$, $Y_{B_{t+1}}$ be the system input and output after a change, the change to its upper neighbor can be shown in the following equations:

(12)$X_{B_{t+1}}=X_{B_t}+\mathsf{\Delta }{x}_A$

(13)$Y_{B_{t+1}}=Y_{B_t}+\mathsf{\Delta }{y}_A$

Similarly, when the state of the system changes to its lower neighbor, the binary state of the $b_Z$ digit is changed from 1 to 0, the paired property of the system is decreased by $(\mathsf{\Delta }{x}_Z,\mathsf{\Delta }{y}_Z)$. Denote $X_{B_t}$, $Y_{B_t}$ be the system input and output of the current state, and $X_{B_{t+1}}$, $Y_{B_{t+1}}$ be the system input and output after a change, the change to its lower neighbor can be shown in the following equations:

(14)$X_{B_{t+1}}=X_{B_t}-\mathsf{\Delta }{x}_Z$

(15)$Y_{B_{t+1}}=Y_{B_t}-\mathsf{\Delta }{y}_Z$

The traditional PI hysteresis model with n elementary operators has $2^n$ possible states. A simple example with 3 operators can be shown in Figure 4. The modeling and control accuracy may be improved with an increase in the number of operators.

When the system input is changed to x with the current state $B_T$, the system output can be calculated with the following procedure:

1. If the system input value, x, is greater than or equal to $X_{B_T}$ and less than $X_{B_T}+\mathsf{\Delta }{x}_A$, the value of the system output y is calculated using (16).

2. If the system input value, x, is greater than $X_{B_T}+\mathsf{\Delta }{x}_A$, the state of the system is changed to its upper neighbor by setting the $b_A$ digit from 0 to 1, the $(X_{B_T},Y_{B_T})$ is updated using (12) and (13). The procedure is then repeated.

3. If the system input value, x, is less than $X_{B_T}$ and greater than $X_{B_T}-\mathsf{\Delta }{x}_Z$, the value of the system output y is calculated using (17).

4. If the system input value, x, is is less than $X_{B_T}-\mathsf{\Delta }{x}_Z$, the state of the system is changed to its lower neighbor by setting the $b_Z$ digit, the $(X_{B_T},Y_{B_T})$ is updated using (14) and (15). The procedure is then repeated.

(16)$y=Y_{B_T}+\frac{\mathsf{\Delta }{y}_A}{\mathsf{\Delta }{x}_A}\times (x-\mathsf{\Delta }{X}_{B_T})$

(17)$y=Y_{B_T}+\frac{\mathsf{\Delta }{y}_Z}{\mathsf{\Delta }{x}_Z}\times (x-\mathsf{\Delta }{X}_{B_T})$

2.3. Digitized Dead-Zone Operators

With the dead-zone operators, the modified PI hysteresis model can model the memory-free asymmetric hysteresis behavior. The state of the dead-zone operators can be represented with a number C. The n paired properties can be used to help describe the values of the n dead-zone operators as shown in (18). An example with two operators is shown in Figure 5.

(18)$S=\left[\begin{array}{ccc}\mathsf{\Delta }{x}_1& \cdots & \mathsf{\Delta }{x}_m\\ \mathsf{\Delta }{y}_1& \cdots & \mathsf{\Delta }{y}_m\end{array}\right]$

The model also has a paired property $(X_C,Y_C)$ that corresponds to its state, the relationship between the paired property and the state of the model is described in the following equations:

(19)$X_C=X_0+\sum _{i=1}^C\mathsf{\Delta }{X}_i;C=1,\dots ,m$

(20)$Y_C=Y_0+\sum _{i=1}^C\mathsf{\Delta }{Y}_i;C=1,\dots ,m$

When the system input is changed to x with the current state C, the system output can be calculated with the following strategy:

• 1.. If the system input value, x, is greater than or equal to $X_C$ and less than $X_C+\mathsf{\Delta }{x}_C$, the value of the system output y is calculated using (21).

• 2.. If the system input value, x, is greater than $X_C+\mathsf{\Delta }{x}_C$ and C is less than m, the state of the system is changed to its upper neighbor by increasing C, the $(X_C,Y_C)$ is updated by adding $(\mathsf{\Delta }{x}_C,\mathsf{\Delta }{y}_C)$. The procedure is then repeated.

• 3.. If the system input value, x, is less than $X_C$ and C is greater than 1, the state of the system is changed to its lower neighbor by decreasing C, the $(X_C,Y_C)$ is updated by subtracting $(\mathsf{\Delta }{x}_C,\mathsf{\Delta }{y}_C)$. The procedure is then repeated.

(21)$y=Y_C+\frac{\mathsf{\Delta }{y}_C}{\mathsf{\Delta }{x}_C}\times (x-\mathsf{\Delta }{X}_C)$

2.4. Inverse Model

The key idea of the parameters identification of the inverse model is to find the reflection of the resultant hysteresis curves about the ${45}^{\circ }$ line. A linear response is obtained by cascading the inverse hysteresis operator ${\mathsf{\Gamma }}^{-1}$ as a feedfoward controller with the actual hysteresis. The proposed model for the inverse feedfoward controller is illustrated in Figure 6.

Unlike previous work, the parameters of the inverse model using the proposed representation can be found by simply exchanging $\mathsf{\Delta }x$ with $\mathsf{\Delta }y$ within the paired property, as shown in the following equations:

(22)$H^{-1}=\left[\begin{array}{ccc}\mathsf{\Delta }{y}_1& \cdots & \mathsf{\Delta }{y}_n\\ \mathsf{\Delta }{x}_1& \cdots & \mathsf{\Delta }{x}_n\end{array}\right]$

(23)$S^{-1}=\left[\begin{array}{ccc}\mathsf{\Delta }{y}_1& \cdots & \mathsf{\Delta }{y}_m\\ \mathsf{\Delta }{x}_1& \cdots & \mathsf{\Delta }{x}_m\end{array}\right]$

Due to its recursive nature, the calculation of the parameters identification of the inverse model does not change much with an increase in the number of the operators, while that of the modified PI hysteresis model is proportional to the number of the elementary operators. With the proposed digitized representation, the inversion calculation is avoided.

3. Experimental Results

In this section, the asymmetric hysteresis behavior of a piezoelectric stacked linear actuator, MPO-050100 (Nanofaktur, Stuttgart, Germany), is modeled using the classical representation as well as the proposed representation. Two experiments have been performed. The first experiment is to identify the parameters and to model the hysteresis behavior of the piezoelectric actuator. The second experiment is to test the performance of the proposed representation deploying the inverse model as a feedfoward compensator.

3.1. Experimental Setup

As seen from Figure 7, the input for the driver is a 10 V peak to peak, 1 Hz sinusoidal wave, is produced by a 16-bit D/A card. The input is passed to an amplifier, EMO-050100 (Nanofaktur, Stuttgart, Germany). The piezoelectric stacked linear actuator will deform, and the output displacement will be measured by a laser sensor, ZW-7010 (Omron, Kyoto, Japan) and is then converted to an analog signal which is received by a 16-bit A/D card. The experimental setup is shown in Figure 8.

3.2. Experimental Results

The hysteresis behavior of the piezoelectric stacked linear actuator under periodic control inputs is recorded. The parameters of the digitized representation and the classical representation of the modified PI model with the dead-zone operators are then identified with $n=25$, the identified parameters are used to model the hysteresis behavior of the piezoelectric stacked linear actuator. With the least square fitted models, the modeling results using the classical and the digitized representation with dead-zone operators are superimposed to the measured hysteresis behavior as shown in Figure 9. It is shown that the modeling results are almost the same using classical representation and the digitized representation. The root mean square error (RMSE) between the modeled hysteresis behavior and the measured hysteresis behavior are 0.50 $\mathsf{\mu }$m using either of the two representation with a stroke of 66.14 $\mathsf{\mu }$m. The modeling error is around $0.76\%$ of its stroke. This proves that the proposed representation has almost the same performance compared to the classical representation while modeling the asymmetric actuator hysteresis.

The inverse model is constructed and is then applied to compensate for the hysteresis behavior of the system at a changing amplitude sinusoidal wave with $x_d(t)=A\times sin(2\times \pi \times t)+68.15$ ($\mathsf{\mu }$m). The compensation results are shown in Figure 10. The RMSE between the measured compensated position and the desired position is $0.36$ $\mathsf{\mu }$m, with a maximum amplitude of 60 $\mathsf{\mu }$m. The error is around $0.6\%$.

3.3. Discussion

Experimental results show that the proposed digitized representation performs almost the same as the classical representation while modeling the hysteresis behavior of the piezoelectric stacked linear actuator. This is because the digitized representation is derived from the classical representation and the two representation are equal mathematically. With $n=25$, the RMSE for modeling the piezoelectric actuator hysteresis is around $0.76\%$. The compensation experimental results with the RMSE of around $0.6\%$ is even slightly better than the modeling results. This also proves that the proposed representation can compensate for the asymmetric actuator hysteresis. The description and compensation experimental results prove that the proposed digitized representation and its inverse model can model and compensate for the asymmetric asymmetric actuator hysteresis.

With the proposed digitized representation, the construction of the inverse model is less mathematically complicated and this improvement may be significant with a large value of n. The parameters of the proposed representation is also more intuitive with their geometric meanings; this may help with the parameters identification. The proposed representation is more intuitive and may be easier to be modified further.

4. Conclusions

The proposed digitized representation is an alternative to the classical representation. The inverse model can be easily constructed with the proposed representation and the inversion calculation is avoided. The proposed representation is also more intuitive. The proposed representation and its inverse model are validated with description and compensation experimental results.

Author Contributions

Conceptualization, C.Z., Y.N.A. and W.T.A.; Software, C.F.; Validation, C.Z.; Writing—original draft, C.Z.; Supervision, W.T.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Working principles of the PI hysteresis model.

Figure 2. Example of (shifted) hysteresis obtained with three elementary operators.

Figure 3. The paired properties of the digitized representation with two operators.

Figure 4. A digitized representation of the classical PI hysteresis model with 3 operators.

Figure 5. Paired properties of the dead-zone operators for extending the classical PI hysteresis model for asymmetric hysteresis behavior.

Figure 6. Block diagram of the inverse model as a feedforward controller.

Figure 7. Experimental architecture.

Figure 8. Experimental setup.

Figure 9. Modeling results using the classical representation and the digitized representation with dead-zone operators at 1 Hz, sinusoidal wave.

Figure 10. Compensation results at a changing amplitude sinusoidal wave.