Hysteresis Modeling and Compensation of Piezoelectric Actuators Using Gaussian Process with High-Dimensional Input

: Rate-dependent hysteresis seriously deteriorates the positioning accuracy of the piezoelectric actuators, especially when tracking high-frequency signals. As a widely-used nonparametric Bayesian method, the Gaussian process (GP) has proven its effectiveness in nonlinear hysteresis modeling. In this paper, the dimension of the input to the GP model is extended to consider more dynamic features of the tracking signal so as to improve the rate-dependent hysteresis modeling accuracy. In contrast with the traditional training set containing only the position and speed information, the acceleration and jerk information, as well as their temporal distribution information, is also included in the input of the model. An inverse hysteresis compensator (IHC) is established in the same way, and open-loop and closed-loop controllers are developed by using the IHC. Experimental results on a PEA stage show that with the increase in the input dimension, the hysteresis modeling accuracy improves greatly and, thus, the controllers based on IHC can achieve a better tracking performance.


Introduction
Piezoelectric-driven nano-positioning stages are widely used in ultra-precision machining and measurement, such as fast tool servo (FTS) [1,2] and atomic force microscope (AFM) [3,4]. However, the rate-dependent hysteresis of the piezoelectric actuators (PEAs) seriously deteriorates the positioning accuracy, resulting in the morphological error of FTS and the image distortion of AFM [5]. In addition, as the input voltage frequency increases, the hysteresis loop becomes larger and rounder, which seriously limits the positioning performance, especially in high-speed positioning applications [6]. Therefore, in order to improve the positioning accuracy of the PEAs, it is important to model and compensate for the rate-dependent hysteresis nonlinearity.
In the past few decades, researchers have conducted extensive works to deal with the challenge of hysteretic nonlinearity modeling [7][8][9]. The relevant models can be generally divided into two types: the physics-based model and the phenomenological model [10,11]. The physics-based hysteresis model is derived from the basic physical principles of hysteresis materials via the relationship between physical quantities and empirical formulas [12]. Nevertheless, due to the complex physical causes of the actual nonlinear hysteresis system, it is difficult to establish the model based on physical principles. On the other hand, the phenomenological hysteresis model directly uses mathematical models to describe the nonlinear input and output relationship of hysteresis, without considering its inherent physical characteristics. Because of its relatively simple structure, the phenomenological model has become the most widely used one in hysteresis research, which mainly includes the Preisach model [13,14], Prandtl-Ishlinskii (P-I) model [15,16], Bouc-Wen model [17][18][19], Dahl model [20], Duhem model [21], etc. Unfortunately, most of these models can only describe the rate-independent hysteresis.
Owning to the fact that hysteretic nonlinearity is generally rate-dependent [6], more rate-dependent models were developed to consider the changing rate of the input voltage, including the rate-dependent Preisach model [22] and rate-dependent PI model [23][24][25][26]. Nonetheless, these models usually contain plentiful unknown parameters, which makes the identification more difficult. Therefore, in 2015, Yang et al. [27] introduced a velocity damping mechanism into the traditional PI model, which successfully constructed a modified rate-dependent PI (MPI) model. This model can accurately describe the dynamic hysteresis characteristics with a relatively simple model structure and relatively few unknown parameters. Recently, with the development of machine learning methods, intelligent models are becoming more popular due to their powerful nonlinear approximation capability [28,29]. Researchers have developed various machine learning methods to model the rate-dependent hysteresis, such as neural networks [30], fuzzy system [31], least-squares support vector machines [32], adaptive fuzzy internal model [33] and so on. These intelligent methods have shown great improvements in modeling accuracy, which provide a new way for modeling the hysteresis of the PEAs.
As a widely-used machine learning method, the Gaussian process (GP), which is based on Bayesian probabilistic inference, removes the need for selecting many parameters while retaining the power of nonlinear dynamic description. Therefore, the usage of GP could make the model flexible and accurate without specifying the functional form and the parameters. It is a proper choice for rate-dependent hysteresis modeling. In 2019, Tao et al. [34] firstly applied the GP to the modeling and feedforward compensation of the PEA hysteresis. Since the hysteresis behavior of PEA is rate-dependent, the voltage value and its changing rate were utilized as a two-dimensional input of the training set and the displacement of the PEA was used as the output, which outperforms the aforementioned MPI model and many other classical models. However, for mixed-frequency signals that are widely used in nano-positioning tracking, the model performance is still unsatisfactory. Therefore, it is necessary to develop an effective model for complex tracking signals.
In this work, besides the position and velocity information used in the traditional models, the acceleration and jerk information, as well as their temporal distribution information, is introduced into the input of the training set and the dimension of the input is increased. The effect of input dimension on prediction is studied and the inverse hysteresis model is obtained via exchanging the input and output signals of the experiment. The IHC is established through the inverse hysteresis model, and open-loop and closed-loop controllers are developed by using the inverse hysteresis compensator (IHC). Comparative experimental investigations on a commercial piezo-actuated stage validate the effectiveness of the proposed GP-based model with high-dimensional input and the controllers based on the IHC. The results show that the increase in the input dimension would improve the accuracy of model prediction and the controller, especially for complex tracking signals. Therefore, the tracking performance of the controllers with the proposed model is greatly improved compared to those controllers with the traditional 2-dimensional GP model and the MPI model.

Rate-Dependent Hysteresis Nonlinear Modeling of PEA
For piezoelectric actuators, the mapping between the input voltage and the output displacement is nonlinear, rate-dependent and memorable because of hysteresis. Thus, the critical issue is to predict the output in the case of a particular input under the influence of hysteresis. The input and output of the GP training set can be expressed as where n is the length of the training input and output vectors and the training output y i ∈ R is the measured displacement of the PEA. Different from the traditional two-dimensional input where v i is the input voltage value and . v i is its relative changing rate, in this work, the dimensions of the input are extended to six, nine and twelve, which will be discussed in detail in Section 2.2. With these training sets, the latent function f (x) between the input and output could be mapped via training the GP model and the output displacement for a new input X * out of the training set can be predicted.

Principle of GP-Based Hysteresis Modeling
Essentially, GP is a series of random variables, in which any finite random variables are subject to joint Gaussian distribution over functions, p( f ). Therefore, if two or more points are picked in a function, observation of the outputs at these points will follow a joint Gaussian distribution. From a function-space view, a GP is formally specified by a mean function m(x) and a covariance function k(x, x ), where: (1) Therefore, the latent function of GP can be written as: In order to facilitate derivation, it is necessary to preprocess the data and subtract its mean to obtain the zero-mean distribution. The following inferences assume m(x) = 0. The covariance function k(x, x ), also called the kernel, plays a pivotal role in the Gaussian regression model.
Suppose that the data noise obtained by the sensor is ε, and assume it is independently and identically distributed and its variance is σ 2 n [34], then: and the joint prior distribution is: where K(X, X) is the covariance matrix between all the data sets in input X and it is expressed as: Other entries K(X, X * ), K(X * , X) and K(X * , X * ) have similar definitions. According to the joint Gaussian prior distribution obtained from the observed value, the joint posterior distribution can be calculated as: It is the critical predictive equation for GP regression with: Finally, y * will be taken as the prediction of the displacement of the PEA. As mentioned above, the covariance function k(x, x ) is the kernel function of the GP, which defines the similarity between samples. As one of the most popular covari-ance functions, the squared exponential (SE) function is chosen in this work, which is expressed as: where σ 2 f is the signal variance and l is the length scale. These parameters contained in the covariance function are called hyperparameters, which are defined as θ. Considering the noise term ε in Equation (4), the hyperparameters are: Hence, the logarithmic marginal likelihood is: The above logarithmic marginal likelihood can be maximized iteratively by the numerical optimization methods and the hyperparameters can be obtained using the training input information. For more information about the hyperparametric optimization, readers can refer to [34].

Extension of the Input Dimensions
As a rate-dependent behavior, the hysteresis nonlinearity is dependent on not only the voltage and its changing rate, but also the changing acceleration and even the jerk. Hence, in order to improve the modeling accuracy, it is necessary to increase the input dimension to include more information into the training set. Inspired by the work of Hu et al. to predict the contouring error dynamics of the multi-axis systems with a deep gated recurrent unit (GRU) neural network [35], the input dimension is extended to six, nine and twelve to improve the modeling performance.
Suppose that the voltage of PEA is: where T is the sampling period. By differentiating the voltage signal v, the changing rate . v, the changing acceleration .. v and the changing jerk ... v can be obtained as: . .
Different from the two-dimensional input of the training set, which contains only the voltage and its changing rate information at time t, the six-dimensional input, ninedimensional input and twelve-dimensional input at time t are constructed by considering not only more temporal derivative information but also the temporal distribution information. They are described as follows: It is worth noting that the original signal collected by the displacement sensor includes noise. The differential process is that the displacement difference between two adjacent sampling points divides by the sampling period. If there exists displacement noise, the noise will be enlarged greatly since it is divided by a rather small denominator, the sampling period. As the differential order increases, the noise will be consistently enlarged in the same way. Thus, a low-pass finite impulse response (FIR) filter is utilized to eliminate the noise of the training set.

Experimental Setup
The experiments are conducted on a commercial one-dimension piezoelectric-driven nano-positioning stage (P66X30, Harbin Core Tomorrow Science and Technology Co., Ltd., Harbin, China). Its displacement is realized by the deformation of a flexure hinge guiding mechanism with a travel range from 0~26.52 µm. The stage is driven by a piezo-ceramicmade PEA composed of lead zirconate titanate and an integrated high-resolution strain gauge sensor is used to measure the displacement. A voltage amplifier with an amplification factor of 15 is employed to drive the PEA and a signal conditioner is used to capture the position signal from the strain gauge sensor. The control algorithms are programmed in the environment of MATLAB/Simulink and downloaded to a rapid prototyping system (dSPACE-DS1103) equipped with several 16-bit Analog-to-Digital Converters (ADCs) and the 16-bit Digital to Analog Convertors (DACs) to realize real-time control. The sampling frequency of the experiments is selected as 20 kHz. The whole experimental setup is illustrated in Figure 1a, while its signal flow block diagram is shown in Figure 1b.

Modeling Results with Different Input Dimension
Since the bandwidth of the proportional-integral (PI) controller of the experimental stage is approximately 1700 Hz, 50-1000 Hz is selected as the interesting frequency range. Therefore, to obtain the training set, a chirp signal with a frequency from 25 to 1025 Hz and an amplitude of up to 30 V is selected as the voltage input to drive the PEA, as shown in Figure 2. And the displacement of the PEA stage is captured by the strain gauge sensor, which is used as the output of the training set. Additionally, the MPI rate-dependent model [27], a classical phenomenological model (please refer to the Appendix A for more details), is built to describe the hysteresis nonlinearity via the same dataset for comparison. To evaluate the modeling accuracy, the widely-used normalized root mean squared error (NRMSE) and the relative maximum error (RME) are chosen as the evaluation indexes [7,8,27,34], which are defined respectively as: where y i andŷ i are the true displacements and the model predictions of the PEA, respectively. To evaluate the modeling accuracy with different input dimensions, sinusoidal signals with a frequency of 100 Hz, 200 Hz, 300 Hz, 500 Hz and 1000 Hz and an amplitude of up to 30 V, a two-frequency mixed-signal x(t) = 7.5(sin(120 × 2πt) + sin(180 × 2πt)) + 22.5, a four-frequency mixed-signal x(t) = 6 sin(100 × 2πt) + 3 sin(150 × 2πt) + 3 sin(200 × 2πt) + 6 sin(250 × 2πt) + 22.5 and a triangular wave of 50 Hz and an amplitude of up to 20 V, are chosen as the driven signals, respectively. It is noted that under high frequency inputs, the actual control voltage will be greatly enlarged due to the lightly-damped resonance of the nanopositioning stage, especially for the closed-loop controllers. To avoid the control saturation, this study is dedicated to the range from 0-10 µm as a reference for any range study. The hyperparameters are obtained by a typical GP algorithm implemented in our MATLAB environment, converged after some 200 iterations on the measured data. No specific parameter optimization is applied at this stage of the investigation that is focused on the adequate performance of the GP method. The results obtained with the MPI model and the GP-based model with different input dimensions are listed in Table 1, and a comparison results with different methods is shown in Figure 3.
From Table 1, it can be observed that for the sinusoidal signal inputs, the GP-based models with different input dimensions all outperform the MPI model. With the increase in the input dimension, the prediction accuracy is slightly enhanced. With the increase in the signal frequency, the predictive capability of the MPI model deteriorates rapidly. However, for the GP-based model, the prediction accuracy retains satisfactory even when the input frequency is very high. For the 1000 Hz sinusoidal signal, the NRMSE for the 12-dimensional input GP-based model is just 1.2855%, much smaller than that of the MPI model. For the mixed-frequency signals, the NRMSEs for the GP-based models with 2-dimensional and 6-dimensional inputs are similar to those for the MPI model. With the increase in the input dimension, the prediction accuracy can be greatly improved. The NRMSE for the GP-based model with the 12-dimensional input is only 1.0759% for the four-frequency mixed-signal, which is just 37.26% of that for the MPI model and 35.62% of that for the GP-based model with the 2-dimensional input. For the triangular wave, the GP-based model yields the similar prediction accuracy as compared with the MPI model. The above testing results verify the effectiveness of GP-based rate-dependent hysteresis modeling with high-dimensional input and show that the modeling accuracy is greatly improved with the enhancement of the input dimension. The proposed highdimensional GP-based model outperforms the classical phenomenological MPI model and the traditional 2-dimensional GP-based model for all testing signals. Compared with Refs. [7][8][9], this model method is capable of realizing more accurate prediction results with fewer hyperparameters. It has been proven a great success in modeling complex signals, including mixed-frequency signal and triangular wave. Nonetheless, the effectiveness of modeling these complex signals were not validated in these references.

Controller Design Based on Hysteresis Compensation
As discussed above, the high-dimensional GP model can effectively predict the ratedependent hysteresis. It can be used to correct the control of the PEAs via the IHC designed with the inverse GP model. Based on the IHC, open-loop and closed-loop controllers are constructed and tested to validate the effectiveness of the proposed model in this section.

The Open-Loop Controller Based on IHC
Owing to the fact that the inversion of hysteresis effect is by nature hysteresis loops, the direct inverse hysteresis compensation concept is widely used in the literature to help eliminate the hysteresis effect [36], which is shown in Figure 4. The compensation process can be written as: y out (t) = H H −1 [y r ] (t) (22) where H[·] denotes the rate-dependent hysteresis model, H −1 [·] denotes its inverse model, y r (t) is the reference trajectory and y out (t) is the real output trajectory. The input voltage of the PEA, v f f (t), is generated from the IHC and the IHC can be directly modeled by GP with different input dimensions via interchanging the voltage and displacement as the input and output of the aforementioned GP model since the inverse of the hysteresis model is still a hysteresis model. The hyperparameters selected are the same as for Equation (11). For more details about the GP-based IHC modeling, readers can refer to [34]. To model the IHC so that v f f (t) can be obtained, the training data is chosen as , where x i is a two, six, nine or twelve-dimensional vector containing the PEA output displacement y i and its temporal derivative and temporal distribution information. For example, the 12-dimensional input x at time t is expressed as:  For comparison, the inverse MPI model is built using the same training data as that used for the two-dimensional GP-based inverse hysteresis model. Sinusoidal sig-nals with different frequencies from 100 Hz to 1000 Hz, a two-frequency mixed-signal x(t) = 7.5(sin(120 × 2πt) + sin(180 × 2πt)) + 22.5, a four-frequency mixed-signal x(t) = 6 sin(100 × 2πt) + 3 sin(150 × 2πt) + 3 sin(200 × 2πt) + 6 sin(250 × 2πt) + 22.5 and a triangular wave of 50 Hz are chosen as the reference trajectories, respectively. The testing results of the open-loop controllers with different compensators are listed in Table 2. From Table 2, it can be observed that the GP-based compensators can provide better tracking accuracy for the sinusoidal signals as compared with the MPI compensator, especially when the tracking frequency becomes high. The 12-dimensional and 9-dimensional GP-based compensators outperform the 6-dimensional and 2-dimensional GP-based compensators. For the mixed-frequency signals and triangular wave, the 12-dimensional GP-based compensator and 9-dimensional GP-based compensator show better tracking performance as well. For the 1000 Hz sinusoidal signal, the NRMSE for the 12-dimensional GP-based compensator is only 34.90% of that for the MPI compensator and 70.21% of that for the 2-dimensional GP-based compensator. For the 120 and 180 Hz mixed-frequency signal, the NRMSE for the 12-dimensional GP-based compensator is just 65.98% of that for the MPI compensator and 49.78% of that for the 2-dimensional GP-based compensator. The evaluation index values summarized in Table 2 strongly validate the effectiveness of the proposed high-dimensional GP-based compensator.
Comparisons of tracking accuracy for some representative signals with different compensators are shown in Figure 5. It can be observed that, even with complex signals, the proposed compensator still retains an excellent tracking accuracy as more information, including more temporal derivation as well as the temporal distribution, is introduced into the training sets compared to the normal 2-dimensional GP-based compensator and the inverse MPI based compensator. The superiority of this method is demonstrated clearly via comparing these tracking curves.

Closed-Loop Controller
Although the GP-based IHC open-loop controllers perform well for the reference trajectories, especially when the dimensions of the input are high, it is difficult for them to cope with the creep of PEA and other kinds of external disturbances. Therefore, in order to improve the stability of the control system and the ability of the disturbance rejection, a feedforward/feedback closed-loop controller is designed based on the IHC, whose block diagram is shown in Figure 6. Here, P(s) denotes the PEA system, and a feedback loop is used to eliminate the influence of the modeling error and various kinds of external disturbances. A widely used PI controller is utilized as the tracking controller C(s). It can be written as: where e(t) is the tracking error, K p is the proportional gain and K i is the integral gain of the PI controller. It is known that the tracking performance of the PI controller becomes better with an increase in the control gain. However, overlarge control gains may cause a low relative stability margin. Hence, the specific parameters of K p and K i are finally maximized with the trial-and-error method in the step response experiments before the unstable vibration occurs. Hence, the input voltage v(t) of PEA is composed of two parts: where v f f (t) is the voltage generated from the IHC and v f b (t) is the voltage generated from the PI controller C(s).
To evaluate the performance of the feedforward/feedback closed-loop controller, the reference trajectories are selected the same as those used in the open-loop tracking tests. The pure PI controller without the feedforward compensator is also tested for comparisons. In the experiments, K p = 0.3 and K i = 20, 000 for all controllers. The testing results are summarized in Table 3. From the table, it can be observed that the GP-based closed-loop controllers outperform the pure PI controller and the MPI-based closed-loop controller. The MPI-based controller performs better than the pure PI controller when tracking low-frequency sinusoidal trajectories, mixed-frequency trajectories and triangular trajectories since hysteresis is compensated to some extent. Nonetheless, when the tracking frequency becomes high, the tracking accuracy of this controller decreases rapidly and is even worse than that of the pure PI controller since it suffers from the rapid deterioration of its modeling accuracy. For the GP-based method, however, the better tracking accuracy can still be maintained at high frequencies due to its smaller modeling errors, especially for the models with 12dimensional and 9-dimensional inputs. It is worth mentioning that for the mixed-frequency signal and triangular wave, with the increase in the input dimension of the model, the tracking errors decrease greatly. For the four-frequency mixed reference trajectory, the NRMSE for the 12-dimensional GP-based method is just 35.93% and 72.70% of those for the pure PI method and the 2-dimensional GP-based method, respectively. For the triangular reference trajectory, the NRMSE for the 12-dimensional GP-based method is just 32.92% and 50.16% of those for the pure PI method and 2-dimensional GP-based method, respectively. In order to better illustrate the performance of the proposed high-dimensional GP-based closed-loop controller, comparisons of tracking performance with different methods and input dimensions are plotted in Figure 7. From this figure, it can be observed that the results with the proposed controller outperform those obtained with the traditional 2-dimensional GP-based controller and the pure PI controller, which demonstrates its priority. In conclusion, the experimental results shown in the table and the figure validate the effectiveness of the high-dimensional GP-based method.

Conclusions
GP regression is a promising method to model rate-dependent hysteretic nonlinearity. However, when tracking high-frequency (i.e., higher than 500 Hz) or complex (i.e., mixedfrequency and triangular) signals, the accuracy of the traditional GP-based model with 2-dimensional input will deteriorate significantly. In this paper, more temporal derivative information, as well as the temporal distribution information, is introduced into the training set by increasing the input dimension of the GP model to improve its modeling accuracy. Experimental results show that with the increase in the input dimension, the prediction accuracy of the model improves greatly, especially for the complex signals, such as the mixed-frequency signal and the triangular wave. For the four-frequency mixed reference trajectory, the NRMSE for the 12-dimensional GP-based method is just 35.93% and 72.70% of those for the pure PI method and the 2-dimensional GP-based method, respectively. As these two types of trajectories are widely used in PEA applications, such as AFM and FTS, the improvement of modeling accuracy is of great importance. GP-based IHCs with different input dimensions are also constructed, and the open-loop and closed-loop controllers based on these IHCs are designed. Testing results show that the 12-dimensional GP-based open-loop and closed-loop controllers both exhibit the best performance among the similar controllers and the tracking performance is much better than that of the classical MPI based controllers, which validates the effectiveness of modeling the rate-dependent hysteresis of PEAs with high-dimensional input GPs.