Transient and Harmonic Unipolar Hysteresis Model of Piezoelectric Actuators Using a System-Level Approach

: Thanks to their integrability and good electromechanical conversion abilities, piezoelectric actuators are a good choice for many actuation applications. However, these elements feature a frequency-dependent hysteresis response that may yield complex control implementation. The purpose of this paper is to provide the extension of a simple hysteresis model based on a system-level approach linking the strain derivative to the driving voltage derivative and taking into account the dynamic behavior of the hysteretic response of the actuator. The proposed enhancement consists of transient and harmonic regimes, allowing to extend the quasi-static model to dynamic behavior with any frequency. In particular, initial strain shift arising from stabilization and accommodation effects as well as frequency-dependent hysteresis shape are considered. The inclusion of the system dynamics in the model is obtained thanks to fractional derivatives and associated fractional transfer functions, allowing the consideration of the full actuator history as well as a ﬁne tuning of the system dynamics over a wide frequency band. Finally, a numerical example of linearized control through compensation loop is provided, demonstrating the interest in the proposed approach for providing a computationally-efﬁcient, simple yet efﬁcient way for ﬁnely predicting the actuator response and thus designing appropriate controllers.


Introduction
Thanks to their high integration potential and energy conversion abilities, piezoelectric materials have been of rising interest over the last decades for application in actuation systems [1,2], such as medical probes [3], inkjet devices [4], or car injectors to name just a few examples. However, when actuated at high electric field, such actuators exhibit nonlinear hysteretic effects that may compromise their correct operation. Therefore, properly designing actuators and their associated control system necessitates good understanding and model of such a hysteresis response.
Many hysteresis models are available in the literature [5], mainly for magnetic devices but with possible extension to piezoelectric systems, which can be possibly divided into two kinds of approaches. The first one consists in numerical models, where hysteresis is constructed from many elementary units, called "hysterons", that are combined together to shape the global actuator's hysteresis response. Such approaches include the well-known Preisach model [6,7] along with its derivatives such as Krasnosel'Skii-Pokrovskii (KP) technique [8,9], as well as Maxwell-Slip methods and associated dry friction models [10,11]. Usually based on phenomenological approach, such models are effective in finely relating the piezoactuator's behavior (including minor loops for instance), but are computationally and memory-intensive, therefore preventing their use in lightweight embedded systems. On the other hands, mathematical techniques, based on explicit or implicit formulations, offer hysteresis models with reduced parameter sets and possibly physical interpretation of their construction. This class includes Jiles-Atherton [12,13], Bouc-Wen [14,15], or Duhem [16,17] approaches for instance. However, proper parameter identification as well as stability and faithful representation of particular hysteretic behaviors may raise some issues when considering these techniques.
However, most of the previous models, at least in their initial formulations, where dedicated to quasi-static regime and did not included transient and harmonic behaviors of the actuators, which can however yield significant changes in the electromechanical response. Therefore, works aiming at enhancing such models by including time-dependent components have been proposed in order to provide more accurate representations of hysteretic behavior of piezoelectric or magnetic elements. This includes, for instance, modification of the Bouc-Wen [18], Jiles-Atherton [19], or dry friction [20] approaches, but also numerical techniques such as Preisach approaches [21]. In addition to dynamical effects such as creep or losses (potentially with nonfractional terms [22][23][24][25]), transient stage can also be related to quasi-static regime, for instance through, stabilization and accommodation that relate the partial repolarization due to irreversible polarization process that occurs in ferroelectric materials [26].
In that purpose, this paper proposes an enhancement of a simple mathematical approach inspired by the butterfly loop of piezoelectric materials previously introduced by the authors [27]. Compared to previously exposed models, this approach provides a simple and compact formulation of the hysteretic operator, allowing lightweight implementation with respect to hysteron-based models and part of mathematically-based ones, as well as efficient hysteretic behavior representation including for instance minor loops. This model is besides mainly dedicated to devices while the large majority of other approaches target material characteristics. However, this model, considering the coefficient relating the derivative of the strain with respect to the input voltage derivative, has so far been developed for quasi-static regime without any inclusion of dynamic effects directly within the model (hybrid approach was proposed in [28]). The reported study therefore proposes to directly include in the model relaxation effects as well as repolarization process for providing a generalized method allowing faithful response prediction as driving voltage frequency is varying and/or waveforms get more complex than a pure sine wave. The paper is decomposed as follows. Section 2, after recalling the original model basics, describes the proposed model refinement that takes into account both transient and harmonic effects. Then, Section 3 aims at validating this approach by comparing with experimental measurements performed on a stack actuator. As an example of application, Section 4 provides a numerical implementation of the present model in a compensation loop for precise position control, with the actuator modeled using a previously experimentally validated approach based on hybridization of the quasi-static model with Bouc-Wen approach and filtering [28]. Finally, Section 5 aims at recalling the main points of the paper and concludes about the model validity and potentials.

Model Principles
Based on the original model exposed in [27], this section aims at investigating and including transient and harmonic effects in order to relate the actuator behavior for more general excitation cases (waveform, frequency. . . ).

System-Level Model Basics
The principles of the original hysteresis model for piezoelectric actuators lie in considering the coefficient α linking the strain derivative dS to the voltage derivative dV as dS = αdV (1) Obviously, if α is a constant, linear relationship holds. However, introducing voltage dependency in α and inducing a butterfly shape of this coefficient as a function of the voltage leads to hysteresis in the voltage-strain plane. Specifically, this butterfly form is achieved by applying an effective voltage that equals the driving voltage V minus a memory voltage V shi f t , set to the driving voltage value when the later reaches a local extremum, or equivalently when its derivative cancels. The coefficient α is then obtained through a constant value α 0 , relating the linear part of the actuator, and a hysteretic part achieved through a predefined mathematical function h: Note that the choice of the mathematical function allows obtaining different hysteresis shapes.
Specifically, considering the absolute effective voltage V − V shi f t in the function yields a symmetric hysteretic behavior, which will be assumed all along this study. Moreover, excluding the linear part α 0 from the hysteretic one represented by h leads to the required property h(0) = 0 for properly relating the small signal behavior. Considering the simple case of linear function such as α(V) = α 0 + γ V − V shi f t as an example yields plots depicted in Figure 1, where the driving voltage is a unipolar shifted sine wave. As the driving voltage and memory voltage are zero (position (1)), the increase in the driving voltage to position (2) also yields a linear increase of the coefficient α, so that the strain S evolves in a quadratic way. Once in position (2), corresponding to a driving voltage maximum position, the memory voltage is set to this maximal value, so that h V − V shi f t = 0 and α goes back to α 0 (position (3)). Then, the driving voltage decreases, leading to an increase of V − V shi f t (as V shi f t > V) and thus of α, until position (4) where the driving voltage reaches a minimum value (actually 0 V in this example) to which is set V shi f t , as in position (1) where α = α 0 .
(1) (1), (4) (3)  Therefore, this approach provides a quite simple and computationally efficient modeling method, as it only requires a few parameters (typically 1 for the low-field behavior-α 0 -and 1-2 for the function h) and only one memory (V shi f t ) while relating quite well the quasi-static behavior of actuators under various cases [27]. It also provides versatility and adaptivity as the function h can be any function (so that the strain variation is not necessarily quadratic), such as exponential or even Gaussian to get closer to material hysteresis loop [27]. One can note that in the particular case of linear function, the method is quite close to a differential form of Rayleigh hysteresis model [29], but this differential approach permits avoiding issues that arise in unipolar driving when considering the regular Rayleigh approach [30,31].
It has to be recalled here that the model targeting applicative aspects, it mainly aims at relating actuator behaviors under lower excitation magnitude compared to the case of material characterization. In other words, it is considered that compared to the full hysteretic/butterfly shape of the material, only minor material loops are under consideration. Therefore, it is assumed that no significant domain switching occurs for instance. In addition, working quantities in that high-level system view consist in strain and voltage instead of (potentially local) strain and electric field, and global average response of the whole device is under consideration. However, even under these moderate excitations, the actuator still exhibits significant hysteresis behavior.

Dynamic Effects-Electromechanical Transient Regime and Harmonic Relaxation
However, polarization mechanisms in ferroelectrics, at the origin of piezoelectricity, show timedependent effects that are not related by the original model described above. Particular relaxation-like phenomena may yield effects such as stabilization and accommodation, creep, or change in the hysteresis loop with the frequency for instance. This study aims at enriching the previously exposed approach by including dynamical effects. In particular, two phenomena will be under focus. First, transient regime, resulting from depolarization for instance and yielding bias in strain even at rest (zero or constant applied voltage) or with time constants not directly related to the driving voltage frequency, will be under consideration. Such an effect can be seen as a change in the initial condition of the actuator. Therefore, during the first cycles upon the voltage application, a drift in the strain response may be observed, which is not related by the initial model exposed in Section 2.1. This effect can physically be partly explained by irreversible polarization and domain rearrangement. Additionally to this transient regime, frequency-dependent behavior can be observed [23], with the hysteresis shape and thus strain response changing as the voltage changes in frequency (but not in magnitude). The strategy adopted in the present study consists in treating such behaviors as particular relaxation phenomena. Furthermore, as these effects may be related to the electrical polarization state of the material, it will be assumed that these relaxations are applied to the driving voltage only ( Figure 2). These transient and harmonic effects are besides related to the whole history of the actuator. Therefore, these relaxation-related effects cannot be taken into account as a classical, Debye-like relaxation model. Instead, non-integer (i.e., fractional) derivatives are good tools for relating these particular features. Indeed, such operators take into account the whole history, with an unlimited memory, of signals and are therefore well adapted to ferroelectrics. On a physical aspect, they are between energy storage terms (potential or kinetic-i.e., zero and second derivation orders of typical quantities such as strain or voltage for instance) and pure losses (first derivation order of typical quantities), therefore allowing taking into account complex phenomena that occur at the material or system level. An example of such an approach is the well-known Cole-Cole relaxation model in dielectrics [32,33], that is, however, classically only considered in the frequency domain (and without hysteresis as well). It could also be noted that, compared to the present approach, the Cole-Cole representation assesses the material behavior at low-field excitation (dielectric spectroscopy for instance), and cannot represent the hysteretic behavior at a system-level as the presently proposed model.
Mathematically, fractional derivatives arise from pseudo-differential operators, and correspond to the fact that, to obtain a classical first-order derivation, a n-fractional derivative operator D n (with n the possibly non-integer derivation order) has to be applied n times. For example, for n = 1/3, the fractional derivation operator should be applied three times to obtain the first-order derivative of a function f : where s is recalled as the Laplace variable and n is the fractional order (n ∈ IR). in(s) and out(s), respectively, denote the input and output signals of the system represented by the transfer function H(s).  [34,35]. Therefore, it can be seen from Figure 3 that the frequency response can be finely tuned both in terms of magnitude response and phase, allowing changing the attenuation slope and bandwidth in a precise way. Phase tuning through the fractional order n also permits shaping differently the elliptical input-output characteristics and its frequency behavior as shown in Figure 4. Moreover, a remarkable property when considering unipolar input (Figure 4b) is the ability of relating an almost constant shift after a transient stage, while the entire derivative order n = 1 yields offsets that significantly change as the frequency increases. Finally, the possibility of finely controlling the system dynamics can also be observed from the step response depicted in Figure 5, showing that low fractional order expands the settling time and thus well relates long-term phenomena.   Therefore, in order to take into account transient and harmonic effects into the static model and according to Figure 2, the driving voltage V is changed into an equivalent static model input voltage V * that is obtained after two filtering processes: one representing the transient effect whose transfer function is denoted G t and the other harmonic phenomena with an associated transfer function G d . This, therefore, yields in the Laplace domain: with G t (s) and G d (s), respectively, defined as with τ i and λ i (i = t, d) the time constant and fractional order associated to each effect, respectively. Obviously, as the transient effect is expected to have a response spanning a wider time period compared to harmonic losses and conversion, one can expect that τ λ d d .s λ d < τ λ t t .s λ t for relatively small value of s and with λ t < λ d .
Alternatively, instead of computing the voltage response in Laplace domain to go back in time domain (with time variable t) for applying the hysteresis model, it is possible to obtain the equivalent voltage by convolution of the impulse responses g t (t) and g d (t), respectively associated with G t (s) and G d (s): where * is the convolution operator, and g t (t) and g d (t) are defined as [36] g t (t) = 1 with E α,β (z) the Mittag-Leffler function:

Set-Up and Model Identification
This section details an evaluation of the model performance considering measurements performed on an experimental apparatus. As a system-oriented approach, the objective of this experimental validation is to asses to what extent the model provides a flexible tool to be implemented in a wide variety of control systems for instance.
In order to experimentally investigate the validity of the proposed model, a piezoelectric stack actuator (PICMA R P-888.51 from PI TM Ceramic GmbH, Lederhose, Germany, with an internal capacitance of 6 µF and 0-100 V unipolar driving voltage) has been considered. The actuator, that is left in stress-free conditions, is driven by a dSPACE TM , Inc, Paderborn, Germany, DS2003/DS2102 acquisition/control system through an amplifying module (PI TM Ceramic GmbH, Lederhose, Germany, E-503). Strain is monitored using a 120.3 Ω strain gauge (BQ120-3CA, Donghua Testing Technology, Jingjiang, China) with a sensitivity of 1.085 mV·µ −1 connected to a signal conditioning module (Donghua Testing Technology, Jingjiang, China, DH 3840). Finally, an oscilloscope (Tektronix TM , OR, États-Unis, TDS 1002) is used for additional monitoring of actuator driving voltage.
Experimental set-up schematic is depicted in Figure 6. Experiments consisted in varying applied voltage magnitude (50 and 100 V) and frequency (1, 5, 10, 20, and 40 Hz). Exponential function has been considered for the quasi-static model: as it relates the hysteresis behavior of the considered actuator with better accuracy than the linear model [27]. In Equation (10), γ is the maximal variation of the strain derivative-voltage derivative coefficient, and V c refers to the voltage coefficient allowing the control of the saturation occurrence of the function h. Identification procedure first consisted in characterizing static model parameters (α 0 , γ, and V c ) by considering steady state at low frequency (1 Hz). More precisely, the procedure consisted in identifying these parameters through the experimentally evaluated value of α, as depicted in Figure 7. It can be noted that the parameter for α 0 has been found to 6 µ ·V −1 , which is quite far from the manufacturer's datasheet (8.33 µ ·V −1 ). This is explained by the fact that transient effect still arises at this frequency. Hence, next procedures will need to include a refinement for α 0 . Then, still at the same frequency, transient response (strain shift) was used for finding the values of transient regime fractional derivative order (λ t ) and time constant (τ t ), as well as refining α 0 (Figure 8). Finally, high frequency excitation (5-40 Hz) permitted extracting the values of the harmonic regime fractional derivative order (λ d ) and time constant (τ d ). Model parameters are listed in Table 1. It is interesting to note that the best value found for the harmonic regime fractional derivation order, λ d , is greater than 1, meaning that the associated behavior is between losses and kinetic energy, while one may expect it to be between potential energy and losses (i.e., 0 ≤ λ d ≤ 1). However, the fractional derivation order associated to transient regime follows this expectation quite well, with a quite small value well denoting a slowly varying process with frequency. It could also be noted that the use of the strain derivative-voltage derivative coefficient α in the characterization procedure may yield some inaccuracy in the identification procedure as a numerical differentiation is performed (in Figures 7  and 8 a 10-point smoothing is applied, which, however, does not yield a significant distortion, as the sampling frequency of 1 ms is 1000 times the driving signal frequency). However, the fact that the theoretical results match these experimental data over a wide voltage range implies a quite high confidence in the found parameter values.

Harmonic Excitation at Different Amplitudes and Frequencies
First experimental measurement set consisted in applying a 50 or 100 V driving voltage magnitude at varying frequencies. Maximal frequency is limited by the maximal output current of the E-503 amplifier [37], yielding a limit approximately equal to 55 Hz for the considered actuator and maximal driving voltage. Therefore, experimental driving frequency has been limited to 40 Hz to ensure that no distortion occurs.
Experimental results, along with model predictions (where the fractional derivatives were obtained using FOMCON toolbox [34,35]), are depicted in Figure 9. Therefore, results show very good agreement between experimentally measured responses and predicted ones. In particular, the model shows good abilities in relating the change in the hysteresis loop, as also illustrated by Figure 10a that depicts the hysteresis loop area, which appears to follow, in the considered frequency range, a rather linear increase. Such an accurate description is enabled by the fractional derivative approach along with its combination with the static model which permits shaping the hysteresis loops as well as its dependence with the driving voltage magnitude. In addition to this harmonic regime consideration, the transient stage, which yields a bias in the strain response in spite of zero minimal voltage due to accommodation and stabilization effects for instance, is also well related by the proposed approach. This can also be demonstrated by the zero-voltage strain depicted in Figure 10b. Again, this ability of relating transient regime are permitted by the fractional derivative included in the quasi-static hysteresis model; considering a first-order integer derivative would lead to a transient regime that would exhibit too important variations in the strain shift as the frequency is changing.

Varying Voltage Magnitude
In addition to the previously exposed harmonic response, the second experimental investigation aimed at assessing the effect of a non-sinusoidal excitation, by introducing a vanishing envelope to the driving voltage (Figure 11). Considering a center frequency of 1 Hz, results are shown in Figure 12, both in terms of voltage-strain curve and strain time-domain waveforms. Again, good agreement between experimental measurements and model prediction can be observed, with a slight overestimation of the strain when considering the theoretical analysis. Still, the model quite accurately relates the hysteretic effect of the actuator, including minor loops. It can be noted that while the sole quasi-static model would yield minor loops sticking to the main one, the introduction of time-dependent fractional derivatives allows predicting detachment of such minor loops from the main path, as well as smooth strain rate reversal.
Nevertheless, one of the major assets of the enhancement provided by this study is the ability to relate well the strain offset that appears. In particular, time-domain waveforms in Figure 12b show that the proposed approach permits assessing the change in such a strain shift as the voltage magnitude decreases. Moreover, it can be noted that, even when the driving voltage keeps zero value, strain continues to decrease (time longer than 5.5 s), which is also very well predicted by the model as shown in the enlargement in Figure 12b.

Numerical Example: Application to Linearized Control
Based on the previous model that includes static, transient, and harmonic regimes, it is proposed in this section to illustrate the interest in the exposed approach for providing a mean for efficiently and accurately controlling a piezoelectric actuator with respect to a given strain command. This illustrative example will be numerically evaluated by considering an experimentally validated model of the same actuator used in the previous section and exposed in [28], using a Bouc-Wen model for relating the transient regime and a fourth-order transfer function for the harmonic behavior ( Figure 13). Note that compared to the work in [28], the harmonic transfer function has been converted to Laplace domain through bilinear transform as: S(s) S stat (s) = 92144s 3 − 68696 × 10 3 s 2 + 114816 × 10 6 s + 432 × 10 9 s 148196s 3 + 154912 × 10 3 s 2 + 123824 × 10 6 s + 432 × 10 9 (11) with S the actuator output strain and S stat the output of the quasi-static model (made of linear, system-level, and Bouc-Wen blocks). Also, to accommodate with the slight discrepancies between the hybrid model in [28] and actual experimental device, the values of λ t , τ t , and τ d have been marginally changed to 0.2, 32 × 10 −5 s and 2.2 × 10 −3 s, respectively (note that for the transient-related fractional transfer function, this changes the coefficient of the frequency-dependent term from 0.19 s 0.15 to 0.2 s 0.2 ). All the other parameters are kept the same as per Table 1. The principles of the proposed control scheme lie in a compensation approach as exposed in [38]. More precisely, based on the model related in this study (Figure 2), the expression of the strain S in the Laplace domain yields: with H(.) the hysteresis operator (defined by the function h and voltage memory V shi f t , but not including the linear part given by α 0 ). Therefore, considering that S com is the desired output strain, the corresponding required voltage V req can be found from Equation (12) as The corresponding implementation for the control is depicted in Figure 14 (the fractional derivative blocks being also part of the FOMCON toolbox [34,35]), along with the actuator model of [28] presented in Figure 13 and a direct control path consisting of the inversion of the sole linear part. Note that compared to Equation (12), model reduction was performed by including the compensation loop before the inverse transient and harmonic functions, and thus removing the associated direct functions before the hysteresis block.  Figure 15 shows the system response under sine excitation at various frequencies. Therefore, it can be seen that the compensation approach including the proposed model allows significantly limiting the hysteretic response of the actuator and provides a mean of allowing a good and accurate control from a desired strain. Additionally, it can be observed that the proposed approach also addresses well the strain drift, as the strain bias at zero voltage appearing in the direct control case is significantly suppressed when using the compensation technique. Slightly lower and larger strains can be noted however at low and high frequencies respectively, which can however be attributed to slight deviation of the actuator model compared to the real one [28]. Moreover, the case of a sine input strain with exponential vanishing magnitude is depicted in Figure 16. In addition to the significant hysteresis reduction including minor loops, such results also show the ability of ensuring that the strain goes back to zero as the command vanishes towards this value, thanks to the inclusion of the fractional transient term. Such characteristics can also be observed on the driving voltage obtained in the compensation approach and depicted in Figure 17, where the strain shift in open loop is efficiently addressed by computed voltages below 0 V.

Conclusions
This study reported the enhancement of a system-level hysteresis model for piezoelectric actuators based on voltage-dependent strain derivative/voltage derivative coefficient, and using shift and reversal of a mathematical function, therefore providing a simple yet efficient way for assessing the actuator behavior. While the original model works in quasi-static regime only, this paper extended its applicability to transient and harmonic regimes. To this end, fractional transfer functions applied to the driving voltage were considered. Thanks to their ability of taking into account the whole history of the actuator and their capability of providing finely tunable system dynamics, fractional derivatives have been shown to allow taking into account transient strain shift effects as well as frequency-dependence of the hysteresis shape in harmonic regime, without compromising the model simplicity. Finally, a numerical example of the application of this model to linearized actuator control based on compensation approach was provided. Therefore, it was demonstrated that the model provides a computationally-efficient yet accurate way for describing the behavior of piezoelectric actuators (but with possible extension to any other kind of actuators, e.g., electromagnetic) including transient and harmonic aspects, with potential use in embedded controllers for instance. To further enhance the model applicability, other effects such as loading or temperature dependent behaviors can be considered as potential further works. In particular, regarding the latter aspect (temperature dependence), it is possible to take advantage of the proposed model. As an example, the effect of temperature on polarization (e.g., working near Curie temperature) can be included in the hysteretic operator (strain derivative/voltage derivative coefficient α) while the potential additional losses may be better included in the fractional terms.

Conflicts of Interest:
The authors declare no conflict of interest.