# Modeling and Control for Giant Magnetostrictive Actuators with Rate-Dependent Hysteresis

Authors: Ping Liu; Zhen Zhang; Jianqin Mao
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The rate-dependent hysteresis in giant magnetostrictive materials is a major impediment to the application of such material in actuators. In this paper, a relevance vector machine (RVM) model is proposed for describing the hysteresis nonlinearity under varying input current. It is possible to construct a unique dynamic model in a given rate range for a rate-dependent hysteresis system using the sinusoidal scanning signals as the training set input signal. Subsequently, a proportional integral derivative (PID) control scheme combined with a feedforward compensation is implemented on a giant magnetostrictive actuator (GMA) for real-time precise trajectory tracking. Simulations and experiments both verify the effectiveness and the practicality of the proposed modeling and control methods.


---

## Body

## 1. Introduction

Magnetostrictive actuator has broad applications in the super-precision tracking system and microvibration control systems due to its prominent capabilities of high energy densities, large stroke, and fast response. Figure1 illustrates the components of a GMA. The stroke and output force is provided by the Terfenol-D rod in response to a varying magnetic field generated by the surrounding solenoid coils. A biased magnetic field generated by the permanent magnet and the prestress are introduced to produce bidirectional actuation and to improve performance of the Terfenol-D rod, respectively [1].Figure 1
Section view of GMA.However, the GMA’s wider application is restricted due to an inherent property of the magnetostrictive materials, the hysteresis nonlinearity, which can cause undesirable inaccuracy, oscillation, or even instability of the systems. The hysteresis is embodied as a nonlinear relationship between the input voltage and the output displacement and is dependent not only on the amplitude but also on the frequencies of the input signals; hence, the hysteresis system is rate-dependent [2]. As shown in Figure 2, the shape of the hysteresis loops changes significantly with the frequencies of the input signals. Consequently, it is a major challenge to capture the complicated rate-dependent hysteretic behavior precisely, and the modeling, identification, and control of the rate-dependent hysteresis have attracted much attention for its sophisticated nature and extensive applications.Figure 2
Rate-dependent hysteresis of GMA.The existing techniques for the rate-dependent hysteresis modeling can be classified into three categories. The first category is based on phenomenological models. Recent propositions include the modified Preisach model [3], the modified Prandtl-Ishlinskii model [4], and the generalized Prandtl-Ishlinskii model [5]. The second category contains the physical models. For instance, JinHyoung et al. developed the semilinear Duhem model to describe the rate-dependent hysteresis [6, 7]. The third category is the intelligent methods, such as artificial neural network (ANN) [8–10], fuzz tree (FT) [11], and all providing efficient ways to model the rate-dependent hysteresis. Control of smart structures with consideration of hysteresis nonlinearity has also received much attention. A very effective approach in dealing with the hysteresis in control system is to find an accurate model of the hysteresis and then to use its inverse compensation to cancel out the hysteresis [12, 13]. In addition, some feedback control scheme including adaptive control [14], robust control [15], and optimal control [16] were also adopted in the controller synthesis.In this paper, we define the ability for a model to correctly describe the data of some other frequencies not in the training set as the generalization ability of frequency. A universal model with generalization ability of frequency and high precision is always required in a practical application for convenience and efficiency.In this paper, we study the modeling and trajectory tracking control of a GMA with rate-dependent hysteresis nonlinearity under variable input signal frequencies. First of all, an RVM model, which is an inherent online machine learning technique using a more flexible and sparser function without additional regularization parameters, is employed as a universal model to capture the rate-dependent hysteresis. Then, the compensation with a PID feedback algorithm based on the proposed hysteresis model is applied to the GMA to cancel out hysteresis for real-time trajectory tracking control. Experiments are carried out to validate the feasibility and effectiveness of the proposed model and the control schemes.


## 2. Rate-Dependent Hysteresis Model

### 2.1. RVM Model

The relevance vector machine (RVM) introduced by Tipping [17, 18] is a probabilistic model similar to the support vector machine (SVM), but where the training takes place in a Bayesian framework, and predictive distributions of the outputs instead of point estimates are obtained. RVM evades the complexity by producing models that have both structure and parameterization processes that are appropriate to the information content of the data.Given a data set of input-target pairs{xn,yn}n=1N considering scalar-valued target functions only, we follow the standard probabilistic formulation and assume that the targets are samples from the model with additive noise as follows:
(1)y=∑i=1N‍wiK(x,xi)+ɛ=Φw+ɛ,
where ɛ is independent sample from some noise process which is further assumed to be mean-zero Gaussian with variance σ2, Φ is the N×(N+1) design matrix with Φ=[φ1,φ2,…,φN]T, wherein φi=[1,K(xi,x1),…,K(xi,xN)]T, K(x,xi) is kernel function, and w=(w0,…,wN)T is the weight vector. In this work, the radial basis function (RBF) kernel function is used. Due to the assumption of independence of the y, the likelihood of the complete data set can be written as
(2)p(y∣w,σ2)=(2πσ2)-N/2exp{-12σ2∥y-Φw∥2}.With as many parameters in the model as training examples, we would expect maximum likelihood estimation ofw and σ2 from (2) to lead to severe overfitting. To avoid this, a common approach is to impose some additional constraints on the parameters, for example, through adopting a Bayesian perspective, and constrain the parameters by defining an explicit prior probability distribution over them.We encode a preference for smoother (less complex) functions by making the popular choice of a zero-mean Gaussian prior distribution overw:
(3)p(w∣α)=∏i=0N‍N(wi∣0,αi-1)=(2π)-N/2∏i=0N‍αi1/2exp(-(αiwi2)2)
with the hyperparameters α=[α0,…,αN]T. Importantly, there is an individual hyperparameter associated independently with every weight, moderating the strength of the prior thereon. According to Bayes’ rule, the posterior distribution over the weights is thus given by
(4)p(w∣y,α,σ2)=p(w∣y,σ2)p(w∣α)p(y∣α,σ2)=(2π)-(N+1)/2|Σ|-1/2×exp{-(w-μ)TΣ-1(w-μ)2},
where the posterior covariance and mean are, respectively;
(5)μ=σ-2ΣΦT,Σ=(σ-2ΦTΦ+A)-1
with A=diag(α0,α1,…,αN). From the marginal likelihood, the probability distribution of is given by
(6)p(y∣α,σ2)=∫‍p(y∣w,σ2)p(w∣α)dw=(2π)-N/2|C|-1/2exp{-12yTC-1y}
with C=σ2I+ΦA-1ΦT. Values of α and σ2 cannot be obtained in closed form, and here we summarize formulae for their iterative re-estimation. For α, differentiation of (6), equating to zero and rearranging, gives
(7)αinew=γiμi2,
where μi is the ith posterior mean weight from (5), and we have defined the quantities γi by
(8)γi=1-αiΣii
with Σii the ith diagonal element of the posterior weight covariance from (5) computed with the current α and σ2 values. For the noise variance σ2, differentiation leads to re-estimation
(9)(σ2)new=∥y-Φμ∥2N-Σiγi.Note that the in the denominator refers to the number of data examples and not the number of basis functions.The learning algorithm thus proceeds by repeated application of (7) and (9), concurrent with updating of the posterior statistics Σ and μ from (5) until some suitable convergence criteria have been satisfied.In practice, during re-estimation, we generally find that many of theα(i) tend to infinity. From (4), this implies that p(wi∣y,α,σ2) becomes highly peaked at zero; that is, we are a posteriori certain that those wi are zero. The corresponding basis functions can thus be pruned, and sparsity is realized.


### 2.2. Modeling Performance

In this paper, the multivalued mapping is transformed into a single-valued one by employing the current and previous inputs and previous outputs as exogenous inputs, the nonlinear regression model of the form(10)y^(k+1)=f(x(k),x(k-1),…,x(k-m);y(k),y(k-1),…,y(k-n))
which is called feedforward model, NARX model, or series-parallel model [19], where x(k) and y(k) are the input and the true output of the system at discrete time instant k, y^(k+1) is the estimated output of the model at time instant k+1, and m and n are the dynamic orders.The effects of the order (m and n) selection on the modeling results are investigated based on multisimulation. It is found that the larger the model order, the lower the training error. However, the testing error is not monotonic with respect to the order selection. We choose the order (m and n) to make a compromise between the training and testing errors. Then, we choose m=3, n=2.We first apply the sinusoidal scanning signals data from 1 Hz to 100 Hz to obtain the nonlinear system based on RVM, then we test the rate-dependent model using respectively the single-signal data from 1 Hz to 100 Hz and some sets of compound-signal data. Figure3 gives the sinusoidal scanning signals data generated by the MATLAB order and the modeling result of the model data.The input signal and the hysteresis loop of the sinusoidal scanning signal.
(a)
(b)Our experimental data are obtained at a sampling frequency of 10 KHz from GMA manufactured by the Beijing University of Aeronautics and Astronautics [20]. At every frequency, we take 500 pairs of data, (xi,yi), i=1,…,500 from all experimental data and from the model simulation, respectively. The modeling results can be evaluated through the root mean square error (RMSE) and the relative error (RE), which are defined, respectively, as
(11)RMSE=∑i=1N|yi-yi′|2N,RE=∑i=1N|yi-yi′|2∑i=1N|yi|2,
where N is the number of data, yi is the experimental output, and yi′ is the output calculated from the model.Figure4 and Table 1 show the comparison between the rate-dependent hysteresis loop measured by experiments and those simulated based on the RVM modeling method at several single-signal and compound-signal inputs.Table 1
Modeling results.
Frequency (Hz)
RMSE (um)
RE
1
0.4833
0.0287
5
0.3948
0.0232
10
0.4127
0.0243
20
0.5569
0.0333
30
0.5048
0.0304
40
0.4833
0.0296
50
0.5493
0.0347
60
0.7150
0.0429
70
0.9346
0.0501
80
0.7939
0.0519
90
1.0523
0.0568
100
0.9319
0.0552
10/20/40
0.5986
0.0543
30/60/90
0.7654
0.0602Results of RVM modeling.
(a)
(b)
(c)
(d)Remarks
(1) The number of the modeling data must be enough big, then the frequency range of the sinusoidal scanning signal can cover the needed frequency.
(2) The experimental data are compared with the modeling results in Figure4 and Table 1 to reveal the outstanding generalization performance of the universal rate-dependent RVM hysteresis model.


### 2.3. Trajectory Tracking Control

The experimental equipment was constructed to achieve 1-DOF trajectory tracking control as shown in Figure5. The GMA was manufactured by Zhang et al. [20], measuring Φ50 mm × 200 mm with Φ7 mm × 80 mm Terfenol-D rod, and was driven by GF-20 amplifier operating in current mode, which was controlled by a computer with dSPACE (DS1103) control board. The displacement was measured by an eddy current sensor with a 0.1 resolution. The input signal was generated by the computer with dSPACE and transported to the GMA through amplifier; then, the output displacement generated by the GMA was transported by the eddy current sensor to the computer, and then we got the experiment data (xi,yi).Figure 5
The diagram of experiment equipment.The basic idea for controller synthesis for hysteresis system is to design a right inverse model to cancel the hysteresis nonlinearity. However, there is no any exact analytical inverse solution to the RVM model, and there are deviations with the variable input frequency in the direct inverse compensation control. In order to improve the accuracy of tracking control, a PID feedback controller combined with compensation in the feed-forward loop is used for real-time tracking control. Figure6 shows the block diagram of the tracking control system. The hysteresis nonlinearity could be eliminated by the inverse compensation in the feed-forward loop, and the PID feedback controller is adopted to deal with remaining nonlinear uncertainties generated by the inverse compensation with the change on the input frequency, which can be seen as the disturbance.Figure 6
The diagram of the combined control scheme.The nonlinear system (10) is said to be invertible at [x(k-1),…,x(k-m);y(k),…,y(k-n)]T if there is a subset A of Rm+n+1, such that for [x(k-1),…,x(k-m);y(k),…,y(k-n)]T∈A,
(12)f(x1(k),…,x(k-m);y(k),…,y(k-n))≠f(x2(k),…,x(k-m);y(k),…,y(k-n)),
for any distinct inputs x1(k), x2(k) [21]. So we construct the inverse compensation of the model based on
(13)x^(k+1)=f(y(k),y(k-1),…,y(k-n);x(k),x(k-1),…,x(k-m)).Although the inverse compensator is not the analytical inverse solution of the original model, the compensation effect is also good, which can be seen in Figure7. And the small part uncompensated can be treated by the PID controller.Figure 7
The compensation effect of the inverse model.The combined feedback controller, suitable for both single signal and compound signal, is fixed with the different input signal. Figure8 and Table 2 show the trajectory tracking control results of the combined control scheme. For convenience, the time intervals shown in the figures are shorter with the increasing frequency. Note that the sampling frequency of tracking control is the same as that used for modeling.Table 2
Tracking control results.
Amplitude (um)
Frequency (Hz)
RMSE (um)
RE
20
1
0.3226
0.0228
20
5
0.3441
0.0243
20
10
0.3838
0.0271
20
20
0.4913
0.0347
20
30
0.6217
0.0440
20
40
0.7874
0.0557
20
50
0.9699
0.0686
20
60
1.1860
0.0839
20
70
1.3793
0.0975
20
80
1.6374
0.1158
20
90
1.7541
0.1196
20
100
1.8659
0.1245
20
10/20/40
0.4836
0.0592
20
30/60/90
0.9719
0.1192
20
15/35/55/75
0.7664
0.1084Results of combined tracking control.
(a)
(b)
(c)
(d)Remarks
(1) The reference signal should be constrained in the range of the output of the GMA betweenumin and umax; otherwise, the overtopped reference signal may not be attained. To solve the problem, one method is to predict the range of GMA output during the modeling process.
(2) The RVM model we derived is a numerical model, and an analytical form is not available, so the PID controller’s parameters cannot be directly computed. Thus, we adjust the parameters of PID controller through offline simulation.
(3) The tracking control results in Figure8 and Table 2 reveal that the combined controller is effective at both single-input and compound-input frequencies, and the RE all less than 13%, the allowable range of the project.


## 2.1. RVM Model

The relevance vector machine (RVM) introduced by Tipping [17, 18] is a probabilistic model similar to the support vector machine (SVM), but where the training takes place in a Bayesian framework, and predictive distributions of the outputs instead of point estimates are obtained. RVM evades the complexity by producing models that have both structure and parameterization processes that are appropriate to the information content of the data.Given a data set of input-target pairs{xn,yn}n=1N considering scalar-valued target functions only, we follow the standard probabilistic formulation and assume that the targets are samples from the model with additive noise as follows:
(1)y=∑i=1N‍wiK(x,xi)+ɛ=Φw+ɛ,
where ɛ is independent sample from some noise process which is further assumed to be mean-zero Gaussian with variance σ2, Φ is the N×(N+1) design matrix with Φ=[φ1,φ2,…,φN]T, wherein φi=[1,K(xi,x1),…,K(xi,xN)]T, K(x,xi) is kernel function, and w=(w0,…,wN)T is the weight vector. In this work, the radial basis function (RBF) kernel function is used. Due to the assumption of independence of the y, the likelihood of the complete data set can be written as
(2)p(y∣w,σ2)=(2πσ2)-N/2exp{-12σ2∥y-Φw∥2}.With as many parameters in the model as training examples, we would expect maximum likelihood estimation ofw and σ2 from (2) to lead to severe overfitting. To avoid this, a common approach is to impose some additional constraints on the parameters, for example, through adopting a Bayesian perspective, and constrain the parameters by defining an explicit prior probability distribution over them.We encode a preference for smoother (less complex) functions by making the popular choice of a zero-mean Gaussian prior distribution overw:
(3)p(w∣α)=∏i=0N‍N(wi∣0,αi-1)=(2π)-N/2∏i=0N‍αi1/2exp(-(αiwi2)2)
with the hyperparameters α=[α0,…,αN]T. Importantly, there is an individual hyperparameter associated independently with every weight, moderating the strength of the prior thereon. According to Bayes’ rule, the posterior distribution over the weights is thus given by
(4)p(w∣y,α,σ2)=p(w∣y,σ2)p(w∣α)p(y∣α,σ2)=(2π)-(N+1)/2|Σ|-1/2×exp{-(w-μ)TΣ-1(w-μ)2},
where the posterior covariance and mean are, respectively;
(5)μ=σ-2ΣΦT,Σ=(σ-2ΦTΦ+A)-1
with A=diag(α0,α1,…,αN). From the marginal likelihood, the probability distribution of is given by
(6)p(y∣α,σ2)=∫‍p(y∣w,σ2)p(w∣α)dw=(2π)-N/2|C|-1/2exp{-12yTC-1y}
with C=σ2I+ΦA-1ΦT. Values of α and σ2 cannot be obtained in closed form, and here we summarize formulae for their iterative re-estimation. For α, differentiation of (6), equating to zero and rearranging, gives
(7)αinew=γiμi2,
where μi is the ith posterior mean weight from (5), and we have defined the quantities γi by
(8)γi=1-αiΣii
with Σii the ith diagonal element of the posterior weight covariance from (5) computed with the current α and σ2 values. For the noise variance σ2, differentiation leads to re-estimation
(9)(σ2)new=∥y-Φμ∥2N-Σiγi.Note that the in the denominator refers to the number of data examples and not the number of basis functions.The learning algorithm thus proceeds by repeated application of (7) and (9), concurrent with updating of the posterior statistics Σ and μ from (5) until some suitable convergence criteria have been satisfied.In practice, during re-estimation, we generally find that many of theα(i) tend to infinity. From (4), this implies that p(wi∣y,α,σ2) becomes highly peaked at zero; that is, we are a posteriori certain that those wi are zero. The corresponding basis functions can thus be pruned, and sparsity is realized.


## 2.2. Modeling Performance

In this paper, the multivalued mapping is transformed into a single-valued one by employing the current and previous inputs and previous outputs as exogenous inputs, the nonlinear regression model of the form(10)y^(k+1)=f(x(k),x(k-1),…,x(k-m);y(k),y(k-1),…,y(k-n))
which is called feedforward model, NARX model, or series-parallel model [19], where x(k) and y(k) are the input and the true output of the system at discrete time instant k, y^(k+1) is the estimated output of the model at time instant k+1, and m and n are the dynamic orders.The effects of the order (m and n) selection on the modeling results are investigated based on multisimulation. It is found that the larger the model order, the lower the training error. However, the testing error is not monotonic with respect to the order selection. We choose the order (m and n) to make a compromise between the training and testing errors. Then, we choose m=3, n=2.We first apply the sinusoidal scanning signals data from 1 Hz to 100 Hz to obtain the nonlinear system based on RVM, then we test the rate-dependent model using respectively the single-signal data from 1 Hz to 100 Hz and some sets of compound-signal data. Figure3 gives the sinusoidal scanning signals data generated by the MATLAB order and the modeling result of the model data.The input signal and the hysteresis loop of the sinusoidal scanning signal.
(a)
(b)Our experimental data are obtained at a sampling frequency of 10 KHz from GMA manufactured by the Beijing University of Aeronautics and Astronautics [20]. At every frequency, we take 500 pairs of data, (xi,yi), i=1,…,500 from all experimental data and from the model simulation, respectively. The modeling results can be evaluated through the root mean square error (RMSE) and the relative error (RE), which are defined, respectively, as
(11)RMSE=∑i=1N|yi-yi′|2N,RE=∑i=1N|yi-yi′|2∑i=1N|yi|2,
where N is the number of data, yi is the experimental output, and yi′ is the output calculated from the model.Figure4 and Table 1 show the comparison between the rate-dependent hysteresis loop measured by experiments and those simulated based on the RVM modeling method at several single-signal and compound-signal inputs.Table 1
Modeling results.
Frequency (Hz)
RMSE (um)
RE
1
0.4833
0.0287
5
0.3948
0.0232
10
0.4127
0.0243
20
0.5569
0.0333
30
0.5048
0.0304
40
0.4833
0.0296
50
0.5493
0.0347
60
0.7150
0.0429
70
0.9346
0.0501
80
0.7939
0.0519
90
1.0523
0.0568
100
0.9319
0.0552
10/20/40
0.5986
0.0543
30/60/90
0.7654
0.0602Results of RVM modeling.
(a)
(b)
(c)
(d)Remarks
(1) The number of the modeling data must be enough big, then the frequency range of the sinusoidal scanning signal can cover the needed frequency.
(2) The experimental data are compared with the modeling results in Figure4 and Table 1 to reveal the outstanding generalization performance of the universal rate-dependent RVM hysteresis model.


## 2.3. Trajectory Tracking Control

The experimental equipment was constructed to achieve 1-DOF trajectory tracking control as shown in Figure5. The GMA was manufactured by Zhang et al. [20], measuring Φ50 mm × 200 mm with Φ7 mm × 80 mm Terfenol-D rod, and was driven by GF-20 amplifier operating in current mode, which was controlled by a computer with dSPACE (DS1103) control board. The displacement was measured by an eddy current sensor with a 0.1 resolution. The input signal was generated by the computer with dSPACE and transported to the GMA through amplifier; then, the output displacement generated by the GMA was transported by the eddy current sensor to the computer, and then we got the experiment data (xi,yi).Figure 5
The diagram of experiment equipment.The basic idea for controller synthesis for hysteresis system is to design a right inverse model to cancel the hysteresis nonlinearity. However, there is no any exact analytical inverse solution to the RVM model, and there are deviations with the variable input frequency in the direct inverse compensation control. In order to improve the accuracy of tracking control, a PID feedback controller combined with compensation in the feed-forward loop is used for real-time tracking control. Figure6 shows the block diagram of the tracking control system. The hysteresis nonlinearity could be eliminated by the inverse compensation in the feed-forward loop, and the PID feedback controller is adopted to deal with remaining nonlinear uncertainties generated by the inverse compensation with the change on the input frequency, which can be seen as the disturbance.Figure 6
The diagram of the combined control scheme.The nonlinear system (10) is said to be invertible at [x(k-1),…,x(k-m);y(k),…,y(k-n)]T if there is a subset A of Rm+n+1, such that for [x(k-1),…,x(k-m);y(k),…,y(k-n)]T∈A,
(12)f(x1(k),…,x(k-m);y(k),…,y(k-n))≠f(x2(k),…,x(k-m);y(k),…,y(k-n)),
for any distinct inputs x1(k), x2(k) [21]. So we construct the inverse compensation of the model based on
(13)x^(k+1)=f(y(k),y(k-1),…,y(k-n);x(k),x(k-1),…,x(k-m)).Although the inverse compensator is not the analytical inverse solution of the original model, the compensation effect is also good, which can be seen in Figure7. And the small part uncompensated can be treated by the PID controller.Figure 7
The compensation effect of the inverse model.The combined feedback controller, suitable for both single signal and compound signal, is fixed with the different input signal. Figure8 and Table 2 show the trajectory tracking control results of the combined control scheme. For convenience, the time intervals shown in the figures are shorter with the increasing frequency. Note that the sampling frequency of tracking control is the same as that used for modeling.Table 2
Tracking control results.
Amplitude (um)
Frequency (Hz)
RMSE (um)
RE
20
1
0.3226
0.0228
20
5
0.3441
0.0243
20
10
0.3838
0.0271
20
20
0.4913
0.0347
20
30
0.6217
0.0440
20
40
0.7874
0.0557
20
50
0.9699
0.0686
20
60
1.1860
0.0839
20
70
1.3793
0.0975
20
80
1.6374
0.1158
20
90
1.7541
0.1196
20
100
1.8659
0.1245
20
10/20/40
0.4836
0.0592
20
30/60/90
0.9719
0.1192
20
15/35/55/75
0.7664
0.1084Results of combined tracking control.
(a)
(b)
(c)
(d)Remarks
(1) The reference signal should be constrained in the range of the output of the GMA betweenumin and umax; otherwise, the overtopped reference signal may not be attained. To solve the problem, one method is to predict the range of GMA output during the modeling process.
(2) The RVM model we derived is a numerical model, and an analytical form is not available, so the PID controller’s parameters cannot be directly computed. Thus, we adjust the parameters of PID controller through offline simulation.
(3) The tracking control results in Figure8 and Table 2 reveal that the combined controller is effective at both single-input and compound-input frequencies, and the RE all less than 13%, the allowable range of the project.


## 3. Conclusion

The main contribution of this paper is to propose a systematic approach for applications of a rate-dependent hysteresis nonlinear model to precise tracking control of the GMA. We model the rate-dependent hysteresis based on the RVM in a wide frequency range. Based on the model, an inverse compensation in the feed-forward loop combined with a PID feedback controller has been developed. Experimental results demonstrate that the RVM model can capture the rate-dependent hysteresis effects, and our identification and control schemes are effective.



---

*Source: 427213-2013-09-16.xml*