Hysteresis inverse compensation-based model reference adaptive control for a piezoelectric micro-positioning platform

The classic BW model was proposed by Bouc in 1967, and Wen extended the model to describe a large class of hysteretic non-linear systems [34]. The model describe by using the first-order non-linear differential equation

$$\begin{eqnarray} \dot{H}(t) = & \alpha\dot{\upsilon}(t) - \beta \left|\dot{\upsilon}(t) \right|H(t){\left| H(t) \right|^{n - 1}} &- \gamma\dot{\upsilon}(t){\left| H(t) \right|^n}\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where H(t) is the hysteresis non-linear term at time t, $\dot{H}(t)$ is the differential of H(t) with respect to t, υ(t) is the input of the BW model, α, β and γ are the model parameters, which controls the shape and size of the hysteresis loop, respectively. The parameter n control the sharpness of the hysteresis loop. Reasonable selection of the above parameters can simulate different types of hysteretic non-linearity of the PMP platform. Usually let n = 1, and the BW model is expressed as

$$\begin{equation} \begin{cases} {y(t)}{ = du(t)-H(t)} \\ {\dot{H}(t)}{ = \alpha \dot {\upsilon}(t) - \beta \left| {\dot{\upsilon}(t)} \right|H(t) - \gamma \dot {\upsilon}(t)\left| H(t) \right|} \end{cases} \end{equation} \tag{ 2 }$$

where d is the actuator factor, and y(t) is the output of the BW model. Owing to describing the hysteresis of the PMP platform, the parameters α, β, γ and d need to be identified. The intelligent optimization algorithms are often used to identify the BW model [35]. This study propose a CMTPDE algorithm to identify the four unknown parameters.

The classical differential evolution (CDE) algorithm has the advantages of simple structure, fast convergence speed, and few adjustable parameters. However, the CDE algorithm has some limitations such as slow convergence and easy to fall into local minimum points. The CMTPDE algorithm is a kind of the intelligent optimization algorithms, which has ability to improve the global search capability. By using the simultaneous evolution of two populations and sharing of optimal individual resources through information exchange, the speed of optimization is improved, and the local optimum is avoided.

To facilitate the parameters identification of the BW model, equation (2) is rewritten as:

$$\begin{eqnarray} &&y(t + 1) = G(t,y(t),\upsilon(t),\Phi ) \end{eqnarray} \tag{ 3 }$$

where Φ = [d, α, β, γ] is the unknown parameter vector of BW model for the PMP platform. The flow diagram of the CMTPDE algorithm for identifying the parameters of the BW model is shown in figure 3. The steps of the proposed algorithm are as follows:

• (a)  
  Define two populations, and both sizes are D. The maximum iteration of the CMTPDE algorithm is $G_{\textrm{max}}$. The iterations that two populations need to achieve before the exchange of information are both K. The value ranges of the parameters which need to be identified are $[P_{\min},P_{\textrm{max}}]$.
• (b)  
  The initial populations are expressed as

  $$\begin{eqnarray} &&{P^j}(0) = P_{\min }^j + \textrm{rand}(D,{N_p})(P_{\textrm{max}} - P_{\min}) \end{eqnarray} \tag{ 4 }$$

  where N_p is the dimension of the parameters which need to be identified. $p_i^j(g)$ is the ith individual of the jth population in the gth generation, j = 1, 2 and i = 1, 2, ..., D. rand(D, N_p ) is a matrix with D rows and N_p columns, where the matrix element value is random from 0 to 1.
• (c)  
  Define the fitness function

  $$\begin{eqnarray} &&f = \frac{1}{2}\sum\limits_{t = 1}^T {{{\left( {y(t) - {y_p}(t)} \right)}^2}} \end{eqnarray} \tag{ 5 }$$

  where y(t) is the output of BW model at time t, and y_p (t) is the actual output of the PMP platform. T is the number of samples. The optimal values of the objective functions for the two populations are f₁ and f₂, and their individuals are Best₁ and Best₂, respectively.
• (d)  
  Update the two populations respectively by doing mutation, crossover and selection. Define the variation factors of the two population as F^j . Select three individuals x₁, x₂ and x₃ randomly, where $i\neq x_1 \neq x_2 \neq x_3$. The formula of mutation is as following

  $$\begin{eqnarray} &&{h_{i}^j}(g + 1) = {p_{x_1}^j}(g) + F^j({p_{x_2}^j}(g) - {p_{x_3}^j}(g)). \end{eqnarray} \tag{ 6 }$$

  Then define the crossover probability as CR. The formula of crossover is as following

  $$\begin{eqnarray} {v_{i}^j}(g + 1) = \begin{cases} {{h_{i}^j}(g + 1)}&{\mathrm{rand}\ {l_{i}^j} \le CR}\\ {{p_{i}^j}(g)}&{\mathrm{rand}\ {l_{i}^j} \gt CR} \end{cases} \end{eqnarray} \tag{ 7 }$$

  where random number ${l_{i}^j}\in[0,1]$. Next the formula of selection is given as

  $$\begin{equation} {p_i^j}(g + 1) = \begin{cases} {{v_i^j}(g + 1)} & {f(v_{i}^j(g))\lt f(p_{i}^j(g))}\\ {{p_{i}^j}(g)} & {f(v_{i}^j(g))\geq f(p_{i}^j(g))}. \end{cases} \end{equation} \tag{ 8 }$$
• (e)  
  Calculate the objective function and storage the population whose the objective function is smaller. For the first population, back to the third step until Kth iteration. Otherwise, execute the seventh step. For the second population, execute the sixth step.
• (f)  
  For the second population, define P_cm as the probability of crossover and mutation. If ${P_{cm}} \lt 0.5$, the crossover operator is adopted, else the mutation operation is adopted. Then back to the third step until Kth iteration. Otherwise, execute the seventh step.
• (g)  
  Determine the global optimum fitness value f and the global optimum individual Best

  $$\begin{equation} f = \begin{cases} f_1&f_1\lt f_2\\ f_2&f_1\geq f_2 \end{cases} \end{equation} \tag{ 9 }$$

  $$\begin{equation} Best = \begin{cases} Best_1&{Best_1\lt Best_2}\\ Best_2&{Best_1\geq Best_2} \end{cases} \end{equation} \tag{ 10 }$$
• (h)  
  G = G + 1. If $G \lt {G_{\textrm{max} }}$, back to the third step. Otherwise, output the parameter vector of the BW model Φ = Best and stop.

Zoom In Zoom Out Reset image size

To identify the parameters of the BW model, a sinusoidal input signal defined as $y(t) = 18 \textrm{sin}(2\pi t + 1.5\pi)+18$ is adopted to drive the PMP platform. The identified results are d = 1.0988, α =−0.5257, β = 1.0185 and γ =−0.2077.

In order to use the RLS algorithm to identify dynamic linear part, it is necessary to collect data through the PMP platform to obtain a set of signals with varying frequencies. In this paper, a sinusoidal swept input signal with a frequency of 1–100 Hz is used to drive the PMP platform to obtain a swept output signal. The frequency-response curves are shown in figure 4 and the transfer function of the dynamic linear part can be obtained as follows

$$\begin{eqnarray} &&{G_0}(s) = \frac{{15800~{\textrm s} + 10}}{{0.001{s^3} + 3.2{s^2} + 15800~{\textrm s} + 10}} \end{eqnarray} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

According to equation (2), the mathematical expression of the inverse model of the static non-linear part is as following

$$\begin{eqnarray} &&\begin{cases} {u(t) = \frac{1}{{d^{^{\prime}}}}({y_r}(t) + H^{^{\prime}}(t))}\\ &&{\dot H^{^{\prime}}(t) = \alpha^{^{\prime}} \dot y_r(t) - \beta^{^{\prime}} \left| {\dot y_r(t)} \right|H^{^{\prime}}(t) - \gamma^{^{\prime}} \dot y_r(t) \left| H^{^{\prime}}(t) \right|} \end{cases} \end{eqnarray} \tag{ 12 }$$

where x(t) is the output displacement, u(t) is the input voltage, $d^{^{\prime}}$ is the actuator factor, and y_r (t) is the expected displacement. The unknown parameters of the inverse model of the dynamic linear part are still identified by the CMTPDE algorithm, and the results are $d^{^{\prime}} = 0.9050$, $\alpha^{^{\prime}} = -0.4448$, $\beta^{^{\prime}} = 1.0806$, $\gamma^{^{\prime}} = -0.1527$. During the experiment, the input signal passes through the inverse model feed-forward controller first, and then through the PMP platform, the inverse model feed-forward compensation can be achieved.

The BW model-based ICC method compensates for the static hysteresis of the PMP plant. But the disturbance rejection ability of the feed-forward control method is weak. In order to avoid the effect of the system uncertainties and external disturbances [32, 33], to further improve the tracking accuracy, a model reference adaptive controller is designed based on the inverse compensation in this paper.

The PMP platform based on the inverse compensation-based controller can be regard as a linear system which is expressed as

$$\begin{eqnarray} &&{G_p}(s) = \frac{{y_p(s)}}{{u(s)}} = {k_p}\frac{{N(s)}}{{D(s)}} \end{eqnarray} \tag{ 13 }$$

where k_p is the actual gain, and $k_p\gt 0$. $D(s) = s^m + a_{m-1}s^{m_1}+\cdots+a_1s+a_0$, $N(s) = b_ns^n + b_{n-1}s^{n_1}+\cdots+b_1s+b_0$ and $m\gt n$. The purpose of control is to adjust the control rate automatically depended on the error between reference model and actual system, so that the actual output can track the expected output. Adopted (11) as the reference model

$$\begin{eqnarray} &&{G_0}(s) = \frac{{y_m(s)}}{{{y_r}(s)}} = {k_m}\frac{{N(s)}}{{D(s)}} \end{eqnarray} \tag{ 14 }$$

where k_m is the gain of the reference model, $k_m\gt0$.

Define the control rate as $u(s) = {k_c}{y_r}(s)$, where y_r (s) is the expected trajectory of the PMP platform. The error between the output of the reference model and the actual output of the system is expressed as

$$\begin{eqnarray} &&e(s) = {y_m}(s) - {y_p}(s) = ({k_m}-{k_c}{k_p})\frac{{N(s)}}{{D(s)}}{y_r}(s) \end{eqnarray} \tag{ 15 }$$

where k_c is the controller gain, which is adjusted by e(s) online. Then the method of gradient descent is adopted to adjust the control parameter k_c online. In order to make $\mathop {\lim }\limits_{t \to \infty } e(t) = 0$, define objective function as

$$\begin{eqnarray} &&J(t) = {\frac{1}{2}}e^2(t) \end{eqnarray} \tag{ 16 }$$

In this paper, the gradient descent method is adopted as the adaptive law to adjust k_c (t) online. The k_c (t) can be get from

$$\begin{align} k_c(t)& = k_c(t-\tau)-\eta \frac{\partial J}{\partial k_c}\nonumber\\ & = k_c(t-\tau)-\eta \frac{{\partial J}}{{\partial e}}\frac{\partial e}{\partial k_c}\nonumber\\ & = k_c(t-\tau)-\eta e(t)\frac{\partial e}{\partial k_c} \end{align} \tag{ 17 }$$

where the step size η is a positive constant. τ is the sampling time.

From equation (15)

$$\begin{eqnarray} &&\frac{\partial e(s)}{\partial k_c} = -k_p\frac{{N(s)}}{{D(s)}}{y_r}(s). \end{eqnarray} \tag{ 18 }$$

Substituted (14) into (18)

$$\begin{eqnarray} &&\frac{\partial e(s)}{\partial k_c} = -\frac{k_p}{k_m}{y_m}(s). \end{eqnarray} \tag{ 19 }$$

By using the inverse Laplace transform and the definition of derivative

$$\begin{eqnarray} &&\dot{k_c}(t) = \eta e(t)\frac{k_p}{k_m}{y_m}(t) = \mu e(t){y_m}(t) \end{eqnarray} \tag{ 20 }$$

where $\mu = \eta \frac{k_p}{k_m}$ is a positive constant. The control rate of the proposed method is expressed as

$$\begin{eqnarray} &&u(t) = k_c(t)y_r(t) \end{eqnarray} \tag{ 21 }$$

Translate (13) and (15) as the observable state-space equations

$$\begin{eqnarray} &&\begin{cases} {\dot{x} = Ax+{k_p}b{u}}\\ &&{{y_p} = cx} \end{cases} \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&\begin{cases} {\dot{x} = Ax+{k_o}b{y_r}}\\ &&{{e} = cx} \end{cases} \end{eqnarray} \tag{ 23 }$$

where ${k_o}(t) = {k_m}-{k_c}(t){k_p}$.

According to the Kalman–Yacubovich lemma [36, 37]: Define a linear stationary system

$$\begin{eqnarray} &&\begin{cases} \dot{x^{^{\prime}}} = A^{^{\prime}}x^{^{\prime}} + b^{^{\prime}}u^{^{\prime}}\\ &&y^{^{\prime}} = c^{^{\prime}}x^{^{\prime}} \end{cases} \end{eqnarray} \tag{ 24 }$$

which is controllable and considerable. Then the sufficient and necessary condition for the transfer function $G(s) = c^{^{\prime}}{(sI - A^{^{\prime}})^{- 1}}b^{^{\prime}}$ to be a strictly positive real function is that there are positive definite matrices P and Q such that

$$\begin{equation} \begin{cases} A^{^{\prime} T}P + PA^{^{\prime}} = - Q \\ b^{^{\prime} T}P = c^{^{\prime}} \end{cases} \end{equation} \tag{ 25 }$$

Define the Lyapunov function

$$\begin{eqnarray} &&V = \mu^{^{\prime}}{x^T}Px + {{k_o}^2} \end{eqnarray} \tag{ 26 }$$

where $\mu^{^{\prime}}\gt0$. The derivation of equation (26) can be expressed as

$$\begin{eqnarray} &&\dot {V} = \frac{{dV}}{{dt}} = \mu^{^{\prime}}\left(\frac{{d{x^T}}}{{dt}}Px + {x^T}P\frac{{dx}}{{xt}}\right) + 2k_o\frac{{dk_o}}{{dt}} \end{eqnarray} \tag{ 27 }$$

Substitute (24) and (25) into (27)

$$\begin{eqnarray} &&\dot{V} = - \mu^{^{\prime}}{x^T}Qx + 2k_o(\dot k_o + \mu^{^{\prime}} {y_r}{b^T}Px) \end{eqnarray} \tag{ 28 }$$

From (20) and (24), the derivative of k_o can be get from

$$\begin{align} {\dot {k_o}} & = - {{\dot k}_c}{k_p}\nonumber\\ & = - \mu e{y_m}{k_p}\nonumber\\ & = - \mu cx{y_m}{k_p} \end{align} \tag{ 29 }$$

Substitute (25) into (29)

$$\begin{eqnarray} &&\dot k_o = - \mu {b^T}Px{y_m}{k_p} \end{eqnarray} \tag{ 30 }$$

Assumption: For the bounded desired trajectory $y_r(t)\gt0$, the y_m (t) is uniformly bounded. For $t\gt0$, it has $y_m(t)\geq \lambda y_r(t)$, where $\lambda\gt0$.

Remark: From a practical viewpoint, the actual input voltage and output displacement of the PMP platform are bounded and non-negative. So the desired trajectory y_r (t) is bounded and non-negative. According to (14) and inverse Laplace transform, it has y_r (t) = 0 and $y_r(t)\gt0$ when y_m (t) = 0 and $y_m(t)\gt0$, respectively. Therefore, there is $y_m(t)\geq \lambda y_r(t)$ and λ is a positive constant.

On the basis of the Assumption and (30), the (28) can be rewritten as

$$\begin{eqnarray} {\dot{V}} & = {- \mu^{^{\prime}}{x^T}Qx + 2k_o(- \mu {b^T}Px{y_m}{k_p} + \mu^{^{\prime}} {y_r}{b^T}Px)}\nonumber\\ & \leq {- \mu^{^{\prime}}{x^T}Qx + 2k_o{b^T}Px(- \mu\lambda y_r{k_p} + \mu^{^{\prime}} {y_r})} \end{eqnarray} \tag{ 31 }$$

where k_p , µ, $\mu^{^{\prime}}$ and λ are positive constants. According to the Assumption, y_r ≥ 0. So that we can get $\dot{V}\lt0$, while $\mu^{^{\prime}} = \mu\lambda{k_p}$. The closed-loop stability of hysteresis inverse compensation-based MRAC is proved.

To demonstrate the effectiveness of the proposed model and controller, the experiment setup is conducted on a PMP platform (MPT-2MRL102A, Boshi, China). The experiment setup and its architecture are shown as figures 6 and 7, respectively. A computer host is used to provide the MATLAB/Simulink environment to compile the control program which drives the PMP platform. The positioning controller (PPC-2CR0150, Boshi Robotics, China) is consisted of a power amplifier and a position sensor, where the power amplifier is used to magnify the control voltage by 15 times to drive the PMP platform and the position sensor is used to measure the actual output displacement of the PMP platform. The multi-function data acquisition card (PCI-1716, Advantech, China) is employed to accomplish the digital-analog and analog-digital converter between the computer and the PMP platform. The sample frequency of the experiments deployed in this paper is set to 10 KHz.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the effectiveness of the RBW model, a series of comparative experiments between the BW model and the RBW model are carried out on the experiment platform. The sinusoid references with different frequency are employed. Figures 8 and 9 show the experiment results of the proposed model. Compared with the BW model, the output of the RBW model is more close to the real output displacement of the PMP platform. The performance evaluation indices of the proposed model are shown in table 1 and figure 10. Apparently, with the increase of frequency, both the maximum (MAX) error rate and the root-mean-square (RMS) error of RBW model are less than the modeling error of the BW model. It signifies that the RBW model has a stronger ability to approximate the character of the PMP platform than the BW model, especially to describe the rate-dependent behavior of the PMP platform.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to demonstrate the effectiveness of the proposed control method, the sinusoidal signals with frequencies of 1 Hz, 10 Hz, 50 Hz and 100 Hz are adopted to excite the PMP platform. The comparison results of the tracking performance for the IBW model-based ICC method, IBW model based PID control method, IPI model based PID control method, where the PI model is proposed in [38] and [39], and the hysteresis inverse compensation-based MRAC method are shown in the figure 11. In this figure, the tracking error of the proposed control method is less than the other control methods under different frequency sinusoidal expect displacement. The performance evaluation indices of the both controller are shown in the table 2 and figure 12. Compared with the the other control methods proposed in the figure 11, the proposed method has less MAX tracking error rate and the RMS tracking error. Especially, with the increase of frequency, the RMS tracking error of the proposed method is significantly decrease. And at the frequency of 100 Hz, the RMS error of the proposed method reduces by 53.6%, 44.5% and 29.6% compared with the ICC method, IBW model based PID control method and IPI model based PID control method, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, the compared experiments are conducted under the triangular signals with frequencies of 1, 10, 50 and 100 Hz to certify the performance of the proposed control method. The tracking results are shown in the figure 13. The MAX tracking error rates and the RMS errors are described in the table 3 and figure 14. We can observe that the tracking error of the proposed method is less than the other control methods of the triangular signals with different frequencies. Moreover, in contrast with the ICC method, the hysteresis inverse compensation-based MRAC reduces the MAX error by 42.4% and reduces the RMS error by 34.0% at 100 Hz, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further demonstrate the performance of the proposed method, the sinusoidal and triangular signals with damping amplitude and different frequencies are adopted as the desired trajectories for the PMP platform, respectively. The comparison experiment results are shown in the figure 15. Figure 15 (a) is the comparison tracking results with the desired complex sinusoidal signal, where the MAX error rate of the proposed method is 0.84%, which reduces by 45.1%, 14.0% and 35.7% compared with the ICC method, IBW model based PID control method and IPI model based PID control method, respectively. Figure 15(b) shows the comparison tracking results with the desired complex triangular signal, where the tracking error of inverse compensation-based MRAC method is less than the other control method proposed in the figure 15(b) distinctly. Hence, the effectiveness and feasibility of the hysteresis inverse compensation-based MRAC method proposed in this paper are verified.

Zoom In Zoom Out Reset image size