A novel two-input asymmetric dynamic cross-coupled hysteresis modeling and identification for piezoelectric stack actuators

Abstract

Piezoelectric Stack Actuators (PEAs) have been widely used in high-precision positioning system because of its fast response. However, under different excitation voltages and external loads with dynamic frequencies, the width and inclination of PEAs hysteresis curves verified correspondingly, and exhibit asymmetry phenomenon. The main contribution of this paper is to present a two-input asymmetric dynamic cross-coupled hysteresis (TADCH) model to describe the complex nonlinear hysteresis of PEAs in dynamic excitation and loading engineering applications. Moreover, we also develop the corresponding identification method of TADCH model. Finally, several experiments are carried out to verify the accuracy of the proposed model.

Keywords

Prandtl–Ishlinskii hysteresis two-input dynamic cross-coupled asymmetric hysteresis

1. Introduction

Piezoelectric driving micro/nans-positioning systems, for instance, Nanocomposite [1], motion control [2], manipulators [3], and scanning probe microscopy [4], piezoelectric actuators (PEAs) always perform unsatisfactorily because of the complicated hysteresis phenomenon. Some models [5] based on physical physical microscopic mechanisms are proposed to eliminate or suppress nonlinear hysteretic behaviour. However, those methods cannot describe the piezoelectric mechanism well.

Experimental results demonstrate that the width of PEAs hysteresis loops widens with the increase of excitation frequency due to the existence of viscous effects. To solve this problem, rate dependent hysteresis models, such as the classical Preisach model (CPM) [6] and its modified model [7], and the improved Prandtl ishlinskii model (PI) [8], have been studied. The rate dependent, so-called dynamic models, are indeed more accurate than static ones [9]. Besides, the input-output curve of PEAs shows obviously rate and amplitude dependence due to the influence of external loads, which is called dynamic cross-coupled hysteresis. Several mathematical models, for describing the cross-coupled hysteresis nonlinearity of magnetostrictive actuators, were proposed. Among them, the most representative is the static coupled hysteresis model based on CPM.

Recently, based on CPM, lots of researches have been done on the formulas of static coupled hysteresis model [10]. Davino [11], developed a load-related model for magnetostrictive drive by employed the CPM model and the memory-less bivariate function. Also, for low-speed load applications, Zhang [12] proposed a weighted superposition model based on PI operator, which is called single side dead zone operator. However, there are few methods to describe the dynamic cross-coupled hysteresis characteristics of PEAs. The hysteresis representation, based on an intelligent algorithm, is also one of the main research branches at present. In reference [13], a double inputs dynamic model is established to describe the unknown cross-coupled effect of PEAs by using a multilayer neural network. However, this kind of model is not easy to integrate with the existing adaptive technologies, and it is difficult to realize the compensation of nonlinear system with unknown parameters.

Moreover, experiments also prove that the hysteresis loops of PEAs are not symmetric. Janaideh [14] proposed an improved PI model for characterizing hysteresis with saturated and asymmetric attributes in magnetostrictive actuators or shape memory alloys actuators by applying the two nonlinear envelope functions. By introducing an asymmetric backlash operator, an improved PI model [15] was presented. In literature [16], two asymmetric operators are used to describe the increase and the decrease of hysteretic curves separately. In addition, the Coleman–Hodgdon model [17] and the Bouc-Wen model [18] are often used to simulate the asymmetric hysteresis behaviour in PEAs.

The main motivation of this paper is the limited evaluation performance of the effect of dynamic load on PEAs hysteresis [19]. The purpose is to provide a two-input asymmetric dynamic cross-coupled hysteresis model (TADCH) to describe the complex hysteresis in dynamic excitation and loading applications and to develop the corresponding parameter identification processes.

2. Modelling of hysteresis nonlinearity

2.1. Prandtl–Ishlinskii model

The unilateral play operator F_r[u](t): C^m[0, T] → C^m[0, T], for any input u (t) ∈ C^m[0, T], is expressed as $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle F_{r}[u](0)=f_{r}(u(0),0)=w(0),\\ \displaystyle F_{r}[u](t)=f_{r}(u(t),F_{r}[u](t_{i}))=w(t),\quad t_{i}\lt t\lt t_{i+1},∼\text{and}∼0\leq i\leq N-1.\end{array} & & \displaystyle\end{eqnarray}$ (1) Where C^m is a function space with piecewise monotonicity and continuity In the function $f_{r}(u,w)=\max \{u-r,\min \{u,w\}\}$ , r is a nonnegative constant defined as the threshold. Set 0 = t₀ < t₁ < ⋯ < t_N = T is the sub-interval sequence of [0, T], so as to guarantee the input is monotonic in each [t_i, t_i₊₁]. In Eq. ((1)), F_r[u](⋅) is a function of threshold r that dominates the loop width. The relationship between input and output of operator ((1)) has central symmetry, and the width of hysteretic loop widens as the threshold r increases.

The PI model parameterized by the threshold r, is a superposition of F_r[u](⋅), and its formula is written as follows $\begin{eqnarray}\displaystyle H_{PI}(t)=au(t)+\displaystyle \int _{0}^{R}p(r)F_{r}[u,z(r)](t)dr. & & \displaystyle\end{eqnarray}$ (2) Herein, a is a positive constant. The density function p (r)(>0) satisfies $\int _{0}^{R}rp(r)dr<\infty$ that is evaluated by test data. When t = 0, z (r) means the initial value of H_PI(t) in (2) before u (0) is applied. However, the hysteresis curves of PI model are centrosymmetric. Besides, the PI model is employed to characterize static hysteresis owing to rate-independent property. Therefore, this is the main motivation to develop an effective representation method to characterize the asymmetric dynamic cross-coupled hysteresis of PEAs.

2.2. Two input Asymmetric dynamic cross-coupled hysteresis model

An asymmetric unilateral play operator $F_{r}^{ID}[.]$ : C^m [0, T] → C^m[0, T], based on the PI model, is proposed to characterize the asymmetric hysteresis in PEAs $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle F_{r}^{ID}[u](0)=f_{r}^{ID}(u(0),0)=w(0),\\[6.0pt] \displaystyle F_{r}^{ID}[u](t)=f_{r}^{ID}(u(t),F_{r}^{ID}[u](t_{i}))=w(t),\quad t_{i}<t<t_{i+1},∼0\leq i\leq N-1\end{array} & & \displaystyle\end{eqnarray}$ (3) with $\begin{eqnarray}\displaystyle f_{r}^{ID}(u,w)=\max \{I(u)-r,\min \{D(u),w\}\}. & & \displaystyle\end{eqnarray}$ (4) Where, we define I (⋅), D (⋅) as envelope functions, which determine the nonlinear and centrosymmetric properties of (3)–(4). I (⋅) and D (⋅): R¹ → R¹(I (u) < D (u)) denotes the ascending and decreasing case of hysteresis curve, respectively, which can be expressed as follows by employing the inverse hyperbolic functions $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle I(u)=a_{0}\ln (u+b_{0}\sqrt{u^{2}+1})+c_{0},\dot{u}\geq 0\\ \displaystyle D(u)=a_{1}\ln (u+b_{1}\sqrt{u^{2}+1})+c_{1},\dot{u}<0.\end{array}\right. & & \displaystyle\end{eqnarray}$ (5) Due to the introduction of I (⋅) and D (⋅), the asymmetric unilateral play operator possesses the property of noncentric symmetry, so it can be used to describe the properties of saturation and noncentric symmetry properties. a₀ > 0, a₁ > 0, and b₀, b₁ are constants estimated via experimental data.

The mathematical expression of two-input asymmetric dynamic cross-coupled hysteresis model, based on (3)–(5), is proposed for characterizing complex hysteresis of PEAs $\begin{eqnarray}\displaystyle H_{o}(t)=h[u](t)+\int _{0}^{R}p(r,G[l](t))F_{r[\dot{u}](t)}^{ID}[u,z](t)dr+\int _{0}^{R^{\ast }}p^{\ast }({\lambda},G^{\ast }[u](t))F_{{\lambda}[\dot{l}](t)}^{ID^{\ast }}[l,z^{\ast }](t)d{\lambda} & & \displaystyle\end{eqnarray}$ (6) with the dynamic threshold functions $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle r[\dot{u}](t)=r+d|\dot{u}(t)|\\ \displaystyle {\lambda}[\dot{l}](t)={\lambda}+d^{\ast }|\dot{l}(t)|.\end{array}\right. & & \displaystyle\end{eqnarray}$ (7) To ensure the output saturation of the proposed model when the excitation signals are too large, the function h[u](t) is defined in the same form as the envelope function (5) $\begin{eqnarray}\displaystyle h[u](t)=\left\{\begin{array}{@{}l@{}}\displaystyle a_{2}\ln (u+b_{2}\sqrt{u^{2}+1})+c_{2},\quad \dot{u}\geq 0\\ \displaystyle a_{3}\ln (u+b_{3}\sqrt{u^{2}+1})+c_{3},\quad \dot{u}<0.\end{array}\right. & & \displaystyle\end{eqnarray}$ (8) Where a₂, b₂, a₃, and b₃ are constants. Herein, u (t), l (t), and H_o(t) denotes the input voltage, the load, and the displacement of PEAs, or the magnetic field, the stress, and the strain in magnetostrictive actuators, respectively. p (r, G[l](t)), p ∗ (𝜆, G ∗ [u](t)) denote the cross-coupled density functions. Herein, $F_{r[\dot{u}](t)}^{ID}[\cdot ]$ and $F_{{\lambda}[\dot{l}](t)}^{ID^{\ast }}[.]$ are asymmetric unilateral play operators depending on dynamic thresholds $r[\dot{u}](t)$ and ${\lambda}[\dot{l}](t)$ .

3. Identification of TADCH model

3.1. Discretization of TADCH model

For inputs u (t) and l (t), we set h is the sample time, and then the discretized input sequences are u (k), l (k)_{k= 1,2, …,}K = T∕h, with h = t (k) − t (k −1). Based on Eq. (6), the discretization form of TADCH model is presented as $\begin{eqnarray}\displaystyle H_{O}(k)\text{} & =\text{} & \displaystyle h[u](k)+\mathop{\sum }_{i=1}^{n}p_{1}(r_{i})F_{r_{i}[\bar{u}](k)}^{ID}[u,z_{i}](k)+G[l](k)\mathop{\sum }_{i=1}^{n}p_{2}(r_{i})F_{r_{i}[\bar{u}](k)}^{ID}[u,z_{i}](k)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\mathop{\sum }_{i=1}^{n}p_{1}^{\ast }({\lambda}_{i})F_{{\lambda}_{i}[\bar{l}](k)}^{ID^{\ast }}[l,z_{i}^{\ast }](k)+G^{\ast }[u](k)\mathop{\sum }_{i=1}^{n}p_{2}^{\ast }({\lambda}_{i})F_{{\lambda}_{i}[\bar{l}](k)}^{ID^{\ast }}[l,z_{i}^{\ast }](k)\end{eqnarray}$ (9) the weighting functions in Eq. (9) are expressed as $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}p_{k}=A_{k}e^{-B_{k}r}\\ p_{k}^{\ast }=A_{k}^{\ast }e^{-B_{k}^{\ast }{\lambda}}\end{array}\right.,\quad k=1,2. & & \displaystyle\end{eqnarray}$ (10) Where A_i is a positive real number, and B_i denotes an arbitrary constant.

Herein, the discretized $r_{i}[\bar{u}](\cdot ),{\lambda}_{i}[\bar{l}](\cdot )$ in Eq. (9) are defined as follows $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle r_{i}[\bar{u}](k)=d_{0}i+d\left|{\displaystyle \frac{u(k)-u(k-1)}{t(k)-t(k-1)}}\right|\\[10.0pt] \displaystyle {\lambda}_{i}[\bar{l}](k)=d_{0}^{\ast }i+d^{\ast }\left|\frac{l(k)-l(k-1)}{t(k)-t(k-1)}\right|\end{array}\right.,\quad i=1,\ldots ,n.k=1,\ldots ,K. & & \displaystyle\end{eqnarray}$ (11) Based on Eqs (3)–(5), and (10), we can obtain the dynamic asymmetric unilateral play operator.

In Eq. (9), it demonstrates that the numerical realization of TADCH model is difficult, which mainly in the estimations of I (⋅) and D (⋅), h[u](t), the weight functions p_k, $p_{k}^{\ast }$ , and the dynamic threshold functions. Moreover, the analytical form of G[l](⋅), G^∗[u](⋅) are also required. So, the identification method of TADCH model divides into three steps owing to the cross-coupled hysteresis effects.

3.2. The first step identification of TADCH model

Several groups of voltage u (t) with different frequency are applied to PEAs respectively, and without load applied at this time. Therefore, G[l](⋅) in Eq. (9) is 0. In the first two terms of Eq. (9), the estimating parameters X ₁ = (a₀, b₀, c₀, a₁, b₁, c₁, a₂, b₂, c₂, a₃, b₃, c₃, A₁, B₁, d₀, d), based on nonlinear least square method, can be evaluated by minimizing an error square function 𝛹 ( X ₁) $\begin{eqnarray}\displaystyle {\Psi}(X_{1})=\mathop{\sum }_{j=1}^{J}\mathop{\sum }_{k=1}^{K}W_{j}(H_{m}(k)-H_{o}(k))^{2} & & \displaystyle\end{eqnarray}$ (12) with $\begin{eqnarray}\displaystyle W_{j}={\displaystyle \frac{w_{j}}{w_{1}}},\quad j=1,\ldots ,J & & \displaystyle\end{eqnarray}$ (13) and meets the constraints $\begin{eqnarray}\displaystyle (a_{0},c_{0},A_{1},d_{0},d)>0. & & \displaystyle\end{eqnarray}$ (14) Where H_m(⋅) is the experimental data, and H_o(⋅) means the output of TADCH model. J indicates the number of hysteresis curves. K denotes the number of data chosen on each loop. w_j is the ratio of the hysteresis width to maximum output. W_j are weighting coefficients. Similarly, the third term in Eq. (9) can be estimated, under the condition of various frequencies of loads exerted, and without voltage applied, set the parameters vector as $\boldsymbol{X}_{2}=∼(a_{0}^{\ast },b_{0}^{\ast },c_{0}^{\ast },a_{1}^{\ast },b_{1}^{\ast },c_{1}^{\ast },A_{1}^{\ast },B_{1}^{\ast },d_{0}^{\ast },d^{\ast })$ .

3.3. The second step identification of TADCH model

We adopt the “first step” excitation voltage conditions in this step. Meanwhile, a series of constant load l_s is applied to PEAs, for s = 1, …, S. In the case of applied loads or not, the displacement variation can be calculated by the following formula $\begin{eqnarray}\displaystyle {\Delta}H_{m}^{s}(k)=H_{m}^{s}(k)-H_{m}(k),\quad s=1,\ldots ,S & & \displaystyle\end{eqnarray}$ (15) Herein, $H_{m}^{s}(k)$ denotes measured data at specific excitations. ${\Delta}H_{m}^{s}$ presents the measurement of displacement variation.

From Eq. (9), the expression of displacement variation is written as $\begin{eqnarray}\displaystyle {\Delta}H^{s}(k)=G(l_{s})\mathop{\sum }_{i=1}^{n}p_{2}(r_{i})F_{r_{i}[\bar{u}](k)}^{ID}[u,z_{i}](k),\quad s=1,\ldots ,S. & & \displaystyle\end{eqnarray}$ (16) Where, G (l_s) in Eq. (16) represents the variation magnitude of output deformation under l_s. The analytical expression G[l](⋅), based on the magnitude in (15) and the least square curve fitting method, can be obtained conveniently. Similarly, the analytical formula G^∗[u](⋅) is derived in (9) under the condition that u (t) is a slow time-varying function or constant.

3.4. The third step identification of TADCH model

In the third step, $H_{m}^{s}(k)$ is estimated by adopting the signal of the same voltage in the “first step”, and a constant load is applied. Based on Eq. (16), the displacement reduction ${\Delta}H_{m}^{s}(k)$ can be calculated. And then, the parameters X ₃ = (A₂, B₂) of the third term (9) in ΔH^s(k), according to the identified X₁ and G[l](⋅), can be obtained from “the first step”. Herein, the weighting constants are defined as W_j = 1, j = 1, 2, 3. Using the same method, we can get the identification of parameters $\boldsymbol{X}_{4}=∼(A_{2}^{\ast },B_{2}^{\ast })$ of the fifth term in function (9).

4. Experiment and analysis

4.1. Test platform

To test the prediction and compensation performance of TADCH model, the platform and corresponding data acquisition system are setup (as shown in Fig. 1), which include the control system and data acquisition board (DAB). Herein, the DAB possesses a multi-channel sampling/holding board for recording multiple analogue signals simultaneously. Furthermore, a non-contact capacitance sensor, with a resolution of 5 nm (2 Hz), 50 nm (8 kHz), and a maximum range of 0.5 mm, is adopted for measuring the output deformation of the test actuator. Herein, the selected actuator is the multilateral executing agency (wtyd08055) of the 26th Research Institute of China Electronic Technology Corporation. In Z direction, the selected load generator can produce a load signal with the range of 0n to 200N and the frequency range of 0 to 100 Hz. The load sensor adopts the one with high rigidity (1 kn/M), the accuracy of 0.6N (0.03%), and the measurement range of 0–2000 N.

Fig. 1.

Test platform.

4.2. Experiments and analyses

4.2.1. Prediction performance verification

Based on the established platform, we aim to verify the validation of TADCH model by comparison with PI model. Given the proposed three-step parameter identification strategy, the parameters X ₁, X ₂, X ₃, X ₄, and the analytical expressions G[l](⋅) and G^∗[u](⋅) can be determined. Figure 2 shows two sets of comparisons, under the excitation voltages u (t) = 50sin(2π ft), frequency of 1 Hz and 100 Hz, and the applied load 0 N and l (t) = 75sin(2πft) frequency of 1 Hz, respectively.

As shown in Fig. 2(a), it is found that TAHCD model is slightly better than PI in describing asymmetric characteristics under the condition of small excitation voltage frequency, and without external load applied. Furthermore, there is little difference between TAHCD model and PI model to predict the amplitude, slope, and width of loops. In Fig. 2(b), with the increase of dynamic excitation frequency and amplitude, the accuracy of PI model decreases significantly. The prediction error of TAHCD model is smaller than the PI model, especially with the introduction of parameters b_i, i = 0, …,3, the functions G[l](⋅), and G^∗[u](⋅).

Fig. 2.

Comparison results. (a) u (t) = 50sin(2πft), frequency of 1 Hz, and applied 0 N; (b). u (t) = 50sin(2πft), frequency of 100 Hz, and applied s (t) = 75sin(2πft) frequency of 1 Hz.

Fig. 3.

Feedforward compensation control. (a). Linearity analysis; (b). Tracking error.

4.2.2. Compensation performance verification

Based on the similar process within the literature 19, the discretization of the inverse TADCH model can be derived correspondingly. In this section, a typical tracking experiment, by employing the TADCH inverse model, is developed for verifying the validity of the feedforward compensation strategy. The desired trajectory is defined as the amplitude attenuation triangle wave signal with a frequency of 100 Hz. Without losing its generality, a sinusoidal load with the frequency of 1 Hz is selected following the prediction experiments 4.2.1. As shown in Fig. 3(a)–(b) present the comparison results between the expected trajectory and the measured ones. In Fig. 3(a), the input-output hysteresis curves are approximately linear, and the ratio of hysteresis width to the maximum output is less than 1%. Besides, as can be seen from Fig. 3(b), the maximum absolute tracking error can be controlled within 120 nm. Experiments show that the introduction of inverse compensation controller based on TADCH model essentially suppresses the effect of hysteresis on PEAs.

5. Conclusions

By introducing the double dynamic threshold functions, the asymmetric anti-hyperbolic envelope functions, and the continuous cross-coupled density functions, we propose a two-input dynamic asymmetric cross-coupled hysteresis model. In addition, based on the nonlinear least square method and several sets of input-output data under specific excitations, a parameter identification strategy for TADCH model is also developed. Experimental results demonstrate that the TADCH model, compared with the PI model, has an excellent performance in describing the asymmetric cross-coupled hysteresis, under the condition of high frequency and high amplitude double excitations. At the same time, the TADCH model and its inverse based feedforward compensation control error can be limited in the hundreds of milliseconds for actual engineering requirements.

Footnotes

Acknowledgements

The authors would like to acknowledge the financial supports from the NSFC (Natural Science Foundation of China, Nos. 11802118, 51705241), the NSFJP (National Science Foundation of Jiangsu Province, No. BK20170808).

Conflict of interest

All authors (the name of author) declare that they have no conflict of interest or financial conflicts to disclose.

References

Rafiee

X. Q.

Mareishi

, Nonlinear response of piezoelectric nanocomposite plates: Large deflection, post-buckling and large amplitude vibration, International Journal of Applied Mechanics7 (2015).

Tan

Dong

, State Estimation of micro-positioning stage with piezoactuators, IEEE Transactions on Industrial Electronics99 (2017), 1–1.

Wei

and Xiao

T.R.

, Modeling of a three degrees of freedom piezo-actuated mechanism, Smart Materials & Structures26 (2017), 015006.

Lau

Liu

and Chen

, Effect of annealing temperature on the morphology and piezoresponse characterization of poly (vinylidene fluoride- trifluoroethylene) Films via Scanning Probe Microscopy, Advances in Condensed Matter Physics2013 (2013), 1–5.

Chen

J.P.

and Montgomery

T.S.

, A macroscopic theory for the existence of the hysteresis and butterfly loops in ferroelectricity, Ferroelectrics23 (1980), 199–207.

Mayergoyz

I.D.

, Mathematical Models of Hysteresis and their Applications, Elsevier, New York, 2003.

Xiao

and Li

, Modeling and high dynamic compensating the rate-dependent hysteresis of piezoelectric actuators via a novel modified inverse preisach model, IEEE Transactions on Control Systems Technology21 (2013), 1549–1557.

Janaideh

M.A.

and Krejčí

, Inverse rate-dependent Prandtl–Ishlinskii model for feedforward compensation of hysteresis in a piezomicropositioning actuator, IEEE/ASME Transactions on Mechatronics18 (2013), 1498–1507.

and Li

, Dahl Model-based hysteresis compensation and precise positioning control of an XY parallel micromanipulator with piezoelectric actuation, Journal of Dynamic Systems Measurement & Control132 (2010), 558–564.

10.

Dong

Y.Y.

and Wang

H.J.

, Identification and experimental assessment of two-input preisach model for coupling hysteresis in piezoelectric stack actuators, Sensors & Actuators A Physical220 (2014), 92–100.

11.

Davino

Giustiniani

and Visone

, Experimental properties of an efficient stress-dependent magnetostriction model, Journal of Applied Physics105 (2009), 368370.

12.

Zhang

Chen

Q.W.

and Mao

J.Q.

, A generalized stress-dependent Prandtl-Ishlinskii model and its adaptive inverse compensation with model reference for GMA, Control Conference IEEE1 (2011), 535–540.

13.

Y.H.

Mao

J.Q.

and Zhang

, On Generalized dynamic Preisach operator with application to hysteresis nonlinear systems, IEEE Transactions on Control Systems Technology19 (2011), 1527–1533.

14.

Janaideh

M.A.

Rakheja

and Su

C.Y.

, A generalized Prandtl-Ishlinskii model for characterizing the hysteresis and saturation nonlinearities of smart actuators, Smart Materials & Structures18 (2009), 0450011–0450019.

15.

Aguirre

Janssens

and Brussel

H.V.

, Asymmetric-hysteresis compensation in piezoelectric actuators, Mechanical Systems & Signal Processing30 (2012), 218–231.

16.

Jiang

and Qiu

, A Modified Prandtl-Ishlinskii model for modeling asymmetric hysteresis of piezoelectric actuators, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control57 (2010), 1200–1210.

17.

Merry

R.J.E.

Uyanik

Koops

K. R.

, Modeling and compensation of asymmetric hysteresis in a piezo actuated metrological AFM, American Control Conference IEEE (2009).

18.

Qin

H.C.

N.B.

and Chen

, An asymmetric hysteresis model and parameter identification method for piezoelectric actuator, Mathematical Problems in Engineering6 (2014), 1–14.

19.

Zhang

Z.J.

and Dong

Y.Y.

, Asymmetrically dynamic coupling hysteresis in piezoelectric actuators: Modeling identification and experimental assessments, International Journal of Applied Mechanics11 (2019).