On Implementation of the Preisach Model Identification and Inversion for Hysteresis Compensation

A challenge for precise positioning in nanopositioning using smart materials is hysteresis, limiting positioning accuracy. The Preisach model, based on the delayed relay operator for hysteresis modelling, is introduced. The model is identified from experimental data with an input function ensuring information for all input levels. This paper presents implementational issues with respect to hysteresis compensation using the Preisach model, showing the procedure to follow, avoiding pitfalls in both identification and inversion. Issues due to the discrete nature of the Preisach model are discussed, and a specific linear interpolation method is tested experimentally, showing effective avoidance of excitation of vibrational dynamics in the smart material. Experimental results of hysteresis compensation are presented, showing an approximate error of 5% between the reference and measured displacement. Consequences of an insufficient discretization level and a high frequency reference signal are illustrated, showing significant deterioration of the hysteresis compensation performance.


INTRODUCTION
Materials exhibiting both sensing and actuation capabilities, arising from a coupling of mechanical properties with applied electromagnetic fields, are commonly referred to as smart materials Moheimani and Goodwin (2001).These materials include, for instance, piezoelectric materials and shape memory alloys, which are materials that show nonlinear hysteretic behavior Tan and Baras (2005).Hysteresis causes the output to lag behind the input, and is the main form of nonlinearity in piezoelectric materials Devasia et al. (2007).With an accurate description of hysteresis, the effect can be compensated for Eielsen et al. (2012); Iyer and Tan (2009); Leang and Devasia (2006); Croft et al. (2001).
Hysteresis models can be roughly classified into two groups, physical-based and phenomenological-based Tan and Baras (2004); Esbrook et al. (2013).The Jiles-Atherton model of ferromagnetic hysteresis Jiles and Atherton (1986), is an example of a model based on a physical description that can describe hysteresis Rosenbaum et al. (2010).However the derivation of such physical models can be an arduous task, and often result in high order models not suited for practical applications Ismail et al. (2009).Phenomenological models, on the other hand, aim only to approximate the physics, thus giving simpler models.The Preisach operator model Iyer and Tan (2009); Hughes and Wen (1997); Iyer and Shirley (2004); Tan et al. (2001); Tan and Baras (2005); Leang and Devasia (2006); Zhao and Tan (2006) for hysteresis modelling was first regarded as a physical model for hysteresis based on some hypothesis concerning the physical mechanisms of magnetization Mayergoyz (2003).However in later years the model has been recognized as a phenomenologicalbased model mainly due to a more general description.
The idea of Preisach models, including the Preisach operator, the Krasnosel'skii-Pokrovskii (KP) operator Riccardi et al. (2012); Tan and Bennani (2008) and the Prandtl-Ishlinskii (PI) operator Al Janaideh et al. (2011); Macki et al. (1993) is to calculate an output as a sum of weighted basis operators.The Preisach operator is described as a delayed relay element, depicted in Fig. 1, the KP operator as a delayed relay element with final slope and the PI operator as play and stop elements.Consequently, due to their similar structure, these operators are referred to as Preisach-type operators Tan and Bennani (2008); Iyer and Tan (2009); Rosenbaum et al. (2010).This paper will focus specifically on the Preisach operator, for more information about the PI and KP operators the reader is advised to consult Tan and Bennani (2008) as a starting point.
Numerous papers have presented hysteresis compensation by employing the Preisach model, however, there are few references to the difficulties and challenges with the implementation of such a hysteresis compensation scheme.This paper demonstrates important aspects to keep in mind when implementing an identification and inversion scheme based on the Preisach operator.The Preisach operator has been chosen due to its popularity, and its common use in the literature.Several implementational issues are discussed and explained with respect to the Preisach model, showing the procedure to follow for hysteresis compensation, avoiding pitfalls in both identification and inversion.In Tan et al. (2001) interpolation is proposed to avoid discrete inputs from the inversion algorithm, however no method for choosing interpolation points was proposed.This paper experimentally tests a specific linear interpolation method, showing effective avoidance of vibrational dynamics in the smart material.
In order to compensate for hysteresis with the Preisach model, the Preisach density function, which corresponds to a weight for each delayed relay element, must be identified.One of the first methods proposed for identification was to twice differentiate the Everett function, obtained by applying first order reversal inputs to the material Mayergoyz (2003), i.e. changing the input signal to the opposite direction.However this method involves differentiation of a measurement, and will consequently be highly sensitive to measurement noise.The most common method for identification of the Preisach density function is based on the constrained least squares method Iyer and Tan (2009); Iyer and Shirley (2004); Galinaitis et al. (2001); Tan (2002); Eielsen (2012).In contrast to linear systems, hysteresis is a rate-independent effect, i.e. that in order to maximize the information through identification, only the input magnitude needs to be varied, and not the input frequencies Iyer and Shirley (2004).
Compensation of hysteresis is commonly achieved by feed forward of an inverse hysteresis input Devasia et al. (2007), illustrated in Fig. 2.This has, for instance, been done for the Coleman-Hodgdon model Eielsen et al. (2012) and the Preisach model with the construction of a right inverse of the hysteresis model Tan and Baras (2005); Hughes and Wen (1997); Croft et al. (2001).However analytical solutions of Preisachtype operators generally do not exist, with the exception of the PI model described with play operators, which inverse turns out to be a stop-type PI operator Kuhnen (2003); Tan and Bennani (2008).Consequently, the inversion of the Preisach operator often has to be carried out iteratively Iyer et al. (2005).In this paper this is done based on a closest match algorithm Tan et al. (2001).
It is well known that negative phase introduces a phase lag between the input and output, which must not be confused with the lagging behind property of the hysteresis nonlinearity Devasia et al. (2007).Due to this, if a time delay and/or a low pass filter is present, the system will experience a phase lag in addition to the real hysteresis.However such a model is only valid for frequencies close in value to the frequencies used to identify it Stakvik (2014).A time delay of T d = 4.58 × 10 −4 s has, for instance, been reported in a commercial available atomic force microscopy (AFM), where the source of this time delay is speculated to be the displacement sensor in the AFM Ragazzon et al. (2014).
The remainder of this paper is organized as follows.The Preisach model is described in Section 2. Section 3 first states the experimental setup, before the identification procedure is explained together with some implementational aspects.In Section 4 experimental results of hysteresis compensation are presented, studying both the effect of low discretization and high frequency on the reference signal.Finally, in Section 5 some concluding remarks are drawn.

MODELING 2.1 The Preisach Operator
The Preisach operator is built up by the delayed relay operator, shown in Fig. 1, where the output y to an input u, for a pair of thresholds (β, α), with β ≤ α is given as for u(t) < β (1) where t is the time, y(t) is the relay output at time t, defined from the previous output y(t − 1) = ξ where ξ is the state of the relay, and R β,α the output of the relay corresponding to the (β, α) pair.
Further, the output of the Preisach operator is calculated as a weighted superposition of delayed relays.For an applied input u β,α ∈ {+1, −1}, the output of the Preisach operator is given by Mayergoyz (2003) where µ(β, α), called the Preisach density function, is a nonnegative weight function, representing the weights of each hysteron in the Preisach plane where α m and β m refer to the highest and lowest values for α and β respectively.The Preisach plane S can be geometrically divided into two subregions, S + , which corresponds to relays with output value +1, and S − , corresponding to output values of −1.
To understand the changing states of the Preisach plane, assume that at some initial time t 0 , the input u(t 0 ) = u 0 < β m .Hence, all delayed relay elements have output value −1.Further, assume that the input is increased monotonically until some maximum at time t 1 , where the input value is u 1 .The corresponding Preisach plane is illustrated in Fig. 3a, where the boundary between S + and S − is the horizontal line α = u 1 .Then all elements below the memory curve are turned on, that is, they have a α value lower than u 1 , likewise, all hysterons with α > u 1 are turned off.
Next, the input is decreased monotonically until some input u 2 at time t 2 , then all hysterons with a β value higher than u 2 will be switched off.This results in the staircase memory curve shown in Fig. 3b.
Based on the regions S + and S − the output in (2) can be rewritten to where the negative contribution is subtracted from the positive contribution.Based on this equation the rateindependent property of the Preisach operator is observed.The output of (3) only depends on the regions S + and S − , which again is defined from the past sequence of local maxima and minima of the input.Consequently, since the output only depends on the levels of the input, and not the frequency, the Preisach operator is rate-independent.

Discretization
In order to implement the discrete Preisach model numerically, the Preisach plane must be discretized.
A common method of doing this is to partition the Preisach plane into subregions, shown in Fig. 4, where a discretization level of n h = 4 results in n h + 1 input levels.The weight µ(β, α) is then described as a center weight, shown as red dots in the same figure.This implies that if the red dot is a part of the S − region, the contribution from this subregion will be negative, correspondingly, if the dot is part of S + the contribution will be positive.The input levels n = 1 and n = h + 1 refer to β m and α m respectively, which define the range of the Preisach model as shown in Fig. 3.

Experimental setup
The aim of this work is to identify and compensate for hysteresis by running tests on a piezoelectric actuator in a nanopositioning laboratory.The experimental setup consists of a dSPACE DS1103 hardware-inthe-loop system, an ADE 6810 capacitive gauge, an ADE 6501 capacitive probe from ADE Technologies, a Piezodrive PDL200 voltage amplifier, the custommade long-range serial-kinematic nano-positioner from EasyLab (see Fig. 5), and two SIM 965 programmable filters.The capacitive measurement has a sensitivity of 1/5 V/µm and the voltage amplifier has a gain of 20 V/V.The programmable filters were used as reconstruction and anti-aliasing filters.With the DS1103 system, a sampling frequency of 100 kHz is used in all experiments.

Constrained Linear Least Squares Identification
In this paper, all identification schemes were conducted off-line with a constrained linear least squares method.
Recall that the Preisach density function, µ(β, α) is nonnegative, which gives the constraints for the least  squares estimate μ(β, α) ≥ 0, where μ(β, α) is the estimate of µ(β, α).In order to estimate μ(β, α), a time series of measured input and outputs were used to create an over-determined system of linear equations based on a discretized form of ( 2) (4) where the vector consisting of y values is the measured hysteretic output, η 0 is a constant contribution, R i (t j ) is the output calculation for Preisach element i at time j, corresponding to R β,α [u, ξ](t i ) in (2), n t is the number of samples and n q is the number of Preisach elements, given as where n h is the level of discretization.The parameters to be identified, μ and η0 , which is the estimate of η 0 , is then found using the pseudo inverse as where η0 is a measurement bias component, A + is the pseudo inverse of the matrix A in (4), μ is the vector of estimated Preisach weights and y is the measured output of the hysteresis.This identification scheme can be performed by the procedure outlined in the next subsection.
When identifying the Preisach density function, the input used for identification is of vital importance.Choosing an input with an equal or larger number of    input reversals compared to the number of discretization levels, n h , results in at least one reversal for each Preisach operator.Consequently, such a number of input reversals will guarantee that every delayed relay element is switched on and off at least once, providing excitation of all Preisach elements.Such an input is shown in Fig. 6a, where a sinusoidal signal with increasing amplitude provides a sufficient number of input reversals for identifying a Preisach model with n h = 15.An example of an identified Preisach distribution, obtained from the experimental setup, can be seen in Fig. 6b.
Since the hysteresis nonlinearity is rate-independent, the frequency of the input signal does not affect the hysteresis response.However due to the dynamics in the smart material and low pass filters in the experimental setup, the frequency of the identification signal must be restricted.This is the case in piezoelectric actuators, where a creep effect appears for low frequencies, while vibration dynamics are excited at high frequencies.In addition, low pass filters have an increasing negative phase for increasing frequencies.For the experimental setup applied in this work, the creep effect occurs below 5 Hz, while the vibration dynamics have a resonance frequency of about 780 Hz Stakvik (2014).To avoid vibrations, the fundamental frequency of the input should be at least a factor of ten below the resonance frequency, preferably even less Devasia et al. (2007).Consequently, the identification signal should be chosen to avoid both high and low frequency issues.A more detailed discussion of these issues are presented in Croft et al. (2001), where the effect of both creep and vibration are modelled prior to the hysteresis identification.

Aspects of Implementation
The estimated bulk contribution, η0 , in eq. ( 6) can be viewed as a constant bias component outside the range of the Preisach model.To identify both μ and η0 , a repeating two step procedure has been applied.Firstly, the Preisach weight μ is identified with an initial η0 , typically zero, satisfying the constraints.Secondly, the bulk contribution η0 is estimated based on the first estimate of μ.This procedure is repeated until the estimate of η0 converges to some value.In this paper the MATLAB functions lsqnonneg and lsqnonlin are employed to implement this identification procedure.If the bulk contribution η0 is not identified, the Preisach weight function μ has to contain this constant contribution, which can cause a significant degradation of performance of the identified model.
The discrete Preisach model is defined by two parameters, the discretization level n h and the range (β m , α m ).During an identification process it is important that the range of the input signal and the range of the Preisach model correspond.In the rest of this paper it is assumed that all Preisach elements have an initial state of −1, i.e. in the S − region, and that the input resembles the increasing signal shown in Fig. 6a.This ensures that the Preisach elements corresponding to small (β, α) values are excited first, due to the initial negative input.Below, three different input ranges are employed in identification of a Preisach model with discretization level n h = 15, β m = −90 and α m = 90.
• In Fig. 7a a too small input range is applied for the Preisach model.As a result of this, the end diagonal elements in the Preisach distribution is zero.Moreover, such an input range does not excite any border element of the distribution more than maximally once (depending on if all Preisach weights are defined as positive or negative initially), making the model store some of the static bulk contribution η in these border elements.
• By applying a more suitable input range, i.e. [-91,91], all elements of the Preisach distribution are excited, illustrated in Fig. 7b.This also removes  the stored static contribution seen in the border elements in the previous case.To excite all elements in the discrete Preisach model, the maximum amplitude of the input signal has to exceed α m for the states to switch from -1 to 1.This is achieved with the input range in this case.
• When employing a too large input range, as shown in Fig. 7c, the end diagonal elements have to contain the weights for all inputs exceeding (β m , α m ).This will make the model a poor fit for all input values larger than the Preisach model range.
An identified Preisach model will, when provided with an input signal, produce a hysteretic output.This should be done to validate that the identified model captures the measured hysteresis behavior.In Fig. 8a this verification has been done on a sinusoidal signal of 20 Hz, using a Preisach model with discretization level n h = 20.The hysteresis loop between the input voltage and output measurement can be seen in Fig. 8b.These plots also illustrate the discrete nature of the Preisach model, where the output changes whenever the values of the delayed relay elements switches.The steps of the model output vary in step size, resulting from different Preisach weights in the identified model.
The discrete behavior of the Preisach operator, discussed above, creates an error between the measured and modelled output.To reduce this error, the discretization level of the model can be increased, which will improve the fit between the measured output and the model output.Concerns related to a higher discretization level is that the model order increases, which in turn makes both identification and inversion more time consuming.

Inversion Algorithm and Interpolation
As previously mentioned, to compensate for hysteresis the Preisach model must be inverted.However the delayed relay operator does not have an analytical inverse, and therefore requires a numerical approximation.In this paper a closest match algorithm is implemented, exploiting the discrete nature of the model by finding the input value which produces an input as close as possible to the reference signal.Due to the closest match property, the inverted signal will always have the optimal value based on the discrete model.The implementation of the closest match algorithm is outlined in more detail in Tan et al. (2001).
The inverted input of the discrete Preisach model naturally has a discrete behavior, which introduces high frequencies in the inverted signal in the discrete steps.These frequencies can excite the vibration dynamics in stiff structures with little damping, such as devices using piezoelectric actuators, causing significant loss of precision in hysteresis compensation.Consequently, to avoid these vibrations, the input should be smoothed, which can be performed by, for instance, linear interpolation as proposed in Tan et al. (2001).In this paper, a specific method for choosing interpolation points is proposed, illustrated in Fig. 9, where each point to interpolate is chosen at a discrete step of the inverted signal.The discrete input in the figure refers to the inverted signal without interpolation.The resulting inverted input signal significantly reduces the magnitude of the higher frequency components (resulting from the discrete steps), subsequently attenuating the excitation of the vibration dynamics.The performance of a few other interpolation methods are studied in Stakvik (2014).
By employing interpolation for smoothing, the hysteresis compensation scheme cannot be applied for realtime implementation with closed loop control.The interpolation procedure requires the next step of the discrete input, however since the reference signal is calculated in real-time of a controller the interpolation will be affected by a varying time-delay.This time-delay is due to the uncertainty of the future reference signals for the inversion algorithm.This issue can be circumvented by applying a direct feed-forward compensation from the reference trajectory, and in this way avoiding the varying time-delay.

Experimental Results of Hysteresis Compensation
If a suitable input range is applied in the identification of the model, experimental tests of hysteresis compensation can be performed.This section will exemplify some practical issues concerning the choice of discretization level and reference frequency for hysteresis compensation.The following three scenarios will illustrate the consequence of a too small discretization level, a too high reference frequency, while the last scenario exemplifies the best performance of hysteresis compensation achieved by applying the Preisach model.For all experiments an input with 100 reversals is used for identifying the models.The reference signal is a sinusoidal with 10 µm amplitude, with the frequencies of the input signal defined for each scenario.

Low Discretization
Applying a low discretization level on the Preisach model is desirable due to low computational effort in identification and inversion.Fig. 10 shows hysteresis compensation based on a model with a discretization level n h = 10 and reference input as a sinusoidal signal of 5 Hz.The error between the reference and the measurement is significant, in addition, the measurement lags behind the reference, i.e. it does not compensate for the hysteresis sufficiently.By comparing the measurement in Fig. 10 with the interpolation method in Fig. 9, the same decreasing ramp is observed in the maximum of the reference.This implies that for low n h the interpolation method causes a significant error in the extremums of the reference.Additionally, the low discretization level introduces an error in the extremum of the reference that can not be explained by the interpolation method.This error originates from a poor model description of the hysteresis, where the output is modelled incorrectly for these values, despite the fact that the identification was conducted with a larger input range than the maximum amplitude of the reference signal.If the discretization level is increased the model description should be able to capture these variations more precisely.

High Frequency
Due to the poor performance of the low discretization level above, the discretization is increased to n h = 100.The interpolation procedure was introduced to avoid large frequency components in the input signal due to steps in the original signal.However if the frequency of the reference signal is increased, the frequency components of the system also increase, causing vibrations in the piezoelectric actuator.These vibrations are illustrated in Fig. 11 where the reference signal has a frequency of 50 Hz.The measurement reveals that the vibrations have a frequency of about 750 Hz, corresponding to the resonant dynamics of the piezoelectric actuator.Even though there are significant vibrations in the measurement, the amplitude of the error is considerably smaller than for the previous case.This is mainly due to the increase of the discretization level.
In general, vibrations as in Fig. 11, are present when the reference frequency is lower.However for a reference signal with sufficiently low fundamental frequency, the magnitude of the frequency components, of the inverted input around the resonance frequency, is sufficiently low for the vibrations to be dominated by other sources of noise.To apply reference signals with high fundamental frequencies, some kind of vibrational compensation scheme must be performed.

High Discretization and Low Frequency
By combining the high discretization level of n h = 100 with low reference frequency of 5 Hz, a best performance of hysteresis compensation using the Preisach   model is illustrated in Fig. 12. From Fig. 12b the reference-measurement relationship shows that the hysteresis is compensated, and from Fig. 12a, the measurement follows the reference without lagging behind.
The error between reference and measurement is approximately 5% at the maximum amplitude, which can be due to varying environmental conditions, e.g.temperature variations of the piezoelectric material, between the time of identification and compensation.
Another reason for the error can be due to the fact that there are some inaccuracies in the identified model close to the extremums of the reference.This hypothesis is supported by the fact that the measurement has a larger amplitude than the reference.For hysteresis compensation using the Coleman-Hodgdon model, the error has been shown to be less than 1% Eielsen et al. (2012).When comparing with these results, compensation of hysteresis with the Preisach model performs poorly, however, an adaptive identification procedure could increase the compensation performance.

CONCLUDING REMARKS
This paper presents implementational issues of hysteresis compensation by employing the Preisach model, based on the delayed relay operator.Suggestions of how to avoid these issues are presented to obtain a satisfactory performance of the hysteresis compensation.
In particular, aspects of the identification and inversion procedure are shown, and the effect of different model parameters and input ranges are illustrated.In the discrete model, the maximum input range should exceed the range of the Preisach model with a small amount.A proposal of how to estimate both the Preisach weight function and the bulk modulus was given, and results from an identification were shown.
Further, the inversion procedure was explained, fo-cusing firstly on the need for interpolation of the inverted signal to avoid unnecessary high frequencies.A specific method for choosing interpolation points was proposed and illustrated.Issues with low discretization levels and high frequency reference signals were illustrated with experimental results, showing significant errors and vibrations, respectively.At last, a best performance of hysteresis compensation by applying the Preisach model was presented.

Figure 3 :
Figure 3: Illustration of the staircase memory curve of the Preisach Plane.

Figure 4 :
Figure 4: Partition of Preisach plane with center weighting masses.
Estimated distribution of the Preisach weights with n h = 15.

Figure 6 :
Figure 6: Fig.(a) shows an example of an input signal with sufficient reversals for identification of a Prisach model with discretization level n h = 15.Fig.(b) shows the corresponding Preisach distribution based on the input signal in Fig. (a).

Figure 7 :
Figure 7: Showing the detrimental effects on the Prisach weight function using wrong input amplitudes for identification of the discrete model.The model is defined for an input range of [-90,90].

Figure 8 :
Figure 8: Model output when identifying using a 20 Hz sine wave, with 90 V amplitude, with n h = 20 as discretization level.

Figure 9 :
Figure 9: Illustration of the proposed linear interpolation applied in the inverse method.

Figure 12 :Figure 10 :
Figure 12: Best performance of Preisach model hysteresis compensation with discretization level n h = 100.

Figure 11 :
Figure 11: Hysteresis compensation with a too high reference frequency for Preisach model with discretization level n h = 100, where the vibrations are caused by actuator dynamics.