Piezoelectric ceramic materials offer several advantages, including compact size, high precision, and rapid response times. These materials are ideally suited for use in microprocessing technology, precision flow control, and accurate positioning systems. They have been extensively employed across various sectors such as aerospace, medical, energy, and defense^[1]. Nonetheless, piezoelectric ceramics demonstrate hysteresis in their input-output relationship, which can compromise the control accuracy, degrade system stability, and potentially induce system oscillations, thereby hindering further advancements in these materials. Consequently, the research into modeling hysteresis in nonlinear systems and the development of high-performance controllers are of paramount importance.

In terms of hysteresis nonlinear modeling, there are currently three main types of hysteresis nonlinear models: (1) Mechanism-based physical models such as the Duhem model^[2], the Maxwell model^[3], and the Jiles-Atherton model^[4]; (2) Phenomenon-based mathematical models such as the Prandtl-Ishlinskii model^[5], the Krasnoselskii-Pokrovskii model^[6], and the Preisach model^[7]; (3) Intelligent computing models like support vector machines^[8], neural networks^[9], and fuzzy methods^[10]. Reference [11] solved the nonlinear parameters of the electromechanical model for piezoelectric actuators based on the Maxwell physical model, providing a nonlinear model for the actuator. However, this requires in-depth research into complex physical mechanisms and parameters with clear physical significance, making this approach less versatile. Reference [12] employed the Duhem model to establish the hysteresis model for piezoelectric actuators, achieving a model that could accurately describe the relationship between input voltage and output displacement. The Duhem model has a clear functional expression, which makes it more convenient than that in reference [11]. Reference [13] improved the traditional Duhem model by dividing it into two half loops for separate modeling and used spline interpolation and neural network methods for model parameter identification, yielding a more precise piezoelectric actuator hysteresis model. However, it involves many parameters, making the identification and computation process more complex. Reference [14] derived a generalized nonlinear Preisach model applicable to piezoelectric actuators based on the nonlinear Preisach model, which enhances generalizability through phenomenological mathematical modeling. Reference [15] identified the functional relationship between hysteresis elements and frequency for the model, establishing a frequency-dependent Prandtl-Ishlinskii (P-I) model. This model overcomes the disadvantage of the Preisach model lacking an analytical inverse, simplifying computation, but the rate-dependent model requires input frequency determination in advance.

In the control of systems with hysteresis nonlinearity, there are primarily two methods of hysteresis compensation: inverse compensation and closed-loop control. The inverse compensation approach involves mathematically constructing a hysteresis model and its inverse model. By incorporating the inverse model in series before the system, it decouples the hysteresis system to eliminate the effects of hysteresis^[16]. Inverse compensation falls under open-loop control. To suppress various disturbances present in the system, closed-loop feedback (such as PID feedback composite control^[17], internal model control^[18], etc.) can be added on top of inverse compensation. Closed-loop control does not require inverse compensation; instead, it involves considering the hysteresis nonlinearity directly during the design process of the controller^[19]. This direct control method increases the burden on the closed-loop control system to suppress disturbances and the design of nonlinear control methods is complex and difficult to implement. Currently, the controller design is only possible for some nonlinear systems with hysteresis characteristics, presenting certain limitations.

Sliding mode control is a variable structure control method that can design sliding modes independent of external disturbances and the controlled object. As a result, sliding mode controllers can effectively suppress various disturbances and the impact of model uncertainties. This control method features a simple and clear design approach and is easy to implement, making it widely used to solve control problems in complex nonlinear systems. Therefore, sliding mode control strategies can be used to implement tracking control in piezo-positioning systems. The control performance of sliding mode control is closely related to the accuracy of the identification model, and precise tracking of the control system requires an accurate system model.

This study examines a piezo-positioning actuator and introduces a Hammerstein model adept at describing dynamic hysteresis. The static nonlinear part of this model is characterized by a P-I model, and the dynamic linear part is derived through system identification via the Hankel matrix method. This model is phenomenological, offering a straightforward and practical modeling technique that does not necessitate the prior determination of input frequency.

Based on the establishment of the model, this paper first implements the P-I inverse model for serial compensation of the hysteresis characteristics of the piezo-positioning actuator. Then, using approximation methods and notch filter, the dynamic linear part of the Hammerstein model is tuned to facilitate controller design. Considering that the hysteresis nonlinearity is difficult to fully counteract, a sliding mode controller is designed to suppress the residual hysteresis and external disturbances. To enhance the tracking performance of the system, an integral augmentation method is used in the sliding mode control. Finally, real-time control experiments are conducted, and the results show that the proposed control strategy improves the control precision and stability of the piezo-positioning system, and exhibits strong adaptability to different input signals.

Hysteresis characteristics represent intrinsic nonlinear properties of piezoelectric ceramic materials, significantly affecting the positioning accuracy of piezo-positioning actuators in precision applications. These characteristics manifest as different displacement values when the same excitation voltage is applied to a piezo-positioning actuator during voltage ascent and descent.

The hysteresis characteristics of piezoelectric ceramics are categorized into static and dynamic hysteresis nonlinearities. Static hysteresis nonlinearity predominantly exhibits memory effects and multi-valued mapping properties, implying that the output of piezoelectric ceramics at any given time is influenced not only by the current input but also by previous inputs. Additionally, a single input can correspond to multiple outputs, as shown in Fig.1(a) (color online). Dynamic hysteresis nonlinearity, primarily characterized by its frequency-dependent properties, demonstrates minimal variations in nonlinearity at low input frequencies, as shown in Fig.1(b) (color online).

The P-I model belongs to a phenomenological hysteresis model that primarily employs a weighted superposition of finite linear Play operators or linear Stop operators to model hysteresis nonlinearity. This paper establishes a P-I model using the Play operator, thereby simulating the system’s hysteresis characteristics effectively.

The Play operator is defined as:

$ \begin{split} & y(t) = L\left( {x(t),y\left( {{t_i}} \right),r} \right) = \\ & \max \left\{ {x(t) - r,\min \left[ {x(t) + r,y\left( {{t_i}} \right)} \right]} \right\}\quad, \\ \end{split} $ (1)

where: t₀≤···≤t_i≤t≤t_i+1···≤t_m; x(t) is the input signal; r is the threshold of the Play operator.

The initial conditions for the Play operator are:

$ \begin{split} & {y\left( {{t_0}} \right) = L\left( {x\left( {{t_0}} \right),0,r} \right) = } \\ & {\max \left\{ {x\left( {{t_0}} \right) - r,\min \left[ {x\left( {{t_0}} \right) + r,0} \right]} \right\}\quad.} \end{split} $ (2)

Fig.2 illustrates the relationship between the input signal x and the output signal y of the Play operator, manifested as a parallelogram structure.

Upon the weighted superposition of various play operators, the P-I model with hysteresis characteristics is obtained, and its output formula is as follows:

$ \begin{split} & {{\textit{z}}(t) = \mathop \sum \limits_{j = 1}^n {\omega _j}{L_j}\left( {x(t),{y_j}\left( {{t_i}} \right),{r_j}} \right) = } \\ & { \sum \limits_{j = 1}^n {\omega _j}\max \left\{ {x(t) - {r_j},\min \left[ {x(t) + {r_j},{y_j}\left( {{t_i}} \right)} \right]} \right\}\quad,} \end{split} $ (3)

where: r_j denotes the threshold of the j-th Play operator, satisfying r_j≥0; ω_j represents the weight of the j-th Play operator; and n signifies the number of weighted superpositions.

For the P-I model, the weighting coefficients and threshold coefficients can be identified from the experimental data of the system's input and output. By selecting appropriate weighting coefficients and threshold coefficients, the P-I model can describe the actual hysteresis characteristic curve of piezoelectric ceramics.

It is important to note that the P-I model belongs to the category of static nonlinear models, while the actual hysteresis nonlinearity of the piezo-positioning actuator is frequency-dependent. Therefore, for the practical piezo-positioning system, it is essential to fully consider the variations in input frequency. It is necessary to establish a dynamic hysteresis model that is related to the input signal frequency on the basis of the static nonlinear model.

To establish a dynamic hysteresis model for the piezo-positioning system, this paper adopts the structure of the Hammerstein model to model the piezo-positioning system, connecting the static nonlinear model with the dynamic linear model in series. The structure of the Hammerstein model is shown in Fig.3. Here, u(t) represents the excitation voltage, and y(t) represents the generated displacement signal.

For the Hammerstein model of the piezo-positioning system, the static nonlinear part of the model is described by a P-I model, while the dynamic linear part is described by the system model obtained through the Hankel matrix system identification method.

The fundamental principle of system identification method is to utilize real input-output data of the system to obtain an equivalent model of the identified system^[20]. In this paper, the dynamic linear part of the piezo-positioning system's Hammerstein model is solved using the Hankel matrix constructed from the system's impulse response sequence, resulting in a more accurate model. The identification theory is as follows.

Given the input-output data of the system, the autocorrelation sequence and cross-correlation sequence are obtained as R_uu and R_yu, respectively, defined as:

$ \left\{ \begin{gathered} {R_{uu}}(kt) = \frac{1}{N}\sum\limits_{i = k}^{N - 1} {u(it)u(it - kt)} \\ {R_{yu}}(kt) = \frac{1}{N}\sum\limits_{i = k}^{N - 1} {y(it)u(it - kt)} \\ \end{gathered} \right.\quad, $ (4)

where, k = 0, 1, 2, ··· , N-1, N is the number of input sequences in one period, u represents the input data, and y represents the output data.

The state equation of a discrete system is given by:

$ \left\{ \begin{gathered} {\boldsymbol{x}}((k + 1)t) = {{\boldsymbol{A}}_h}{\boldsymbol{x}}(kt) + {{\boldsymbol{B}}_h}{\boldsymbol{u}}(kt) \\ {\boldsymbol{y}}(kt)= {{\boldsymbol{C}}_h}{\boldsymbol{x}}(kt) + {{\boldsymbol{D}}_h}{\boldsymbol{u}}(kt) \\ \end{gathered} \right.. $ (5)

If the impulse response of the discrete linear system is denoted as α(kt), where k = 0, 1, 2, ···, then the relation between the correlation function and the impulse response is:

$ {R_{yu}}(Nt) = \sum\limits_{j = 1}^{ + \infty } {\alpha (jt} ){R_{uu}}(Nt - jt)\quad. $ (6)

The Hankel matrix H is constructed from the impulse response as:

$ {\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {\alpha (t)}&{\alpha (2t)}& \cdots &{\alpha (nt)} \\ {\alpha (2t)}&{\alpha (3t)}& \cdots &{\alpha ((n + 1)t)} \\ \vdots & \vdots & \ddots & \vdots \\ {\alpha (nt)}&{\alpha ((n + 1)t)}& \cdots &{\alpha ((2n - 1)t)} \end{array}} \right]. $ (7)

The relationship between the Hankel matrix obtained from the system's impulse response and the state-space equation is:

$ {\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{C}}_h}} \\ {{{\boldsymbol{C}}_h}{{\boldsymbol{A}}_h}} \\ {{{\boldsymbol{C}}_h}{{\boldsymbol{A}}_h}^2} \\ \vdots \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_h}}&{{{\boldsymbol{A}}_h}{{\boldsymbol{B}}_h}}&{{{\boldsymbol{A}}_h}^2{{\boldsymbol{B}}_h}}& \cdots \end{array}} \right]. $ (8)

The singular values of the Hankel matrix can indicate the importance of each mode of the system, thereby allowing the selection of the system's order based on these singular values. By arranging the singular values in descending order, the system order can be determined based on the locations of significant changes in the singular values.

Therefore, singular value decomposition of the Hankel matrix is performed, and further decomposition of the results yields:

$ \begin{split} {\boldsymbol{H}} = & {\boldsymbol{U}}diag\left\{ {\begin{array}{*{20}{c}} {{\sigma _1}} &\cdots &{{\sigma _n}} \end{array}} \right\}{{\boldsymbol{V}}^{\rm{T}}} =\\ & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_1}}&{{{\boldsymbol{U}}_2}} \end{array}} \right]diag\left\{ {{{\sum} _1},{{\sum} _2}} \right\}{\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{V}}_1}}&{{{\boldsymbol{V}}_2}} \end{array}} \right]^{\rm{T}}}= \\ & {{\boldsymbol{U}}_1}{{\sum} _1}{{\boldsymbol{V}}_1}^{\rm{T}} + {{\boldsymbol{U}}_2}{{\sum} _2}{\boldsymbol{V}}_2^{\rm{T}} \approx {{\boldsymbol{U}}_1}{{\sum} _1}{\boldsymbol{V}}_1^{\rm{T}}\quad, \end{split} $ (9)

where, σ₁ ··· σ_n are the singular values of the Hankel matrix, with σ_i > 0, U₁ = [u₁···u_r], $\displaystyle\sum\nolimits_1 $ = diag{σ₁, ···, σ_r}, and V₁ = [v₁···v_r].

In accordance with formula (8), the matrices C_h and B_h can be selected as follows:

$ \left\{ \begin{gathered} {{\boldsymbol{C}}_h} \approx {{\boldsymbol{U}}_1}diag\left\{ {\begin{array}{*{20}{l}} {\sqrt {{\sigma _1}} }& \cdots &{\sqrt {{\sigma _r}} } \end{array}} \right\}{\text{(First row)}} \\ {{\boldsymbol{B}}_h} \approx diag\left\{ {\begin{array}{*{20}{c}} {\sqrt {{\sigma _1}} }& \cdots &{\sqrt {{\sigma _r}} } \end{array}} \right\}{\boldsymbol{V}}_1^{\mathrm{T}}{\text{(First column)}} \\ \end{gathered} \right.. $ (10)

Furthermore, a new Hankel matrix H₁ can be constructed based on the impulse response：

$ \begin{split} {{\boldsymbol{H}}_1} =& \left[ {\begin{array}{*{20}{c}} {\alpha (2t)}&{\alpha (3t)}& \cdots &{\alpha ((n + 1)t)} \\ {\alpha (3t)}&{\alpha (4t)}& \cdots &{\alpha ((n + 2)t)} \\ \vdots & \vdots & \ddots & \vdots \\ {\alpha ((n + 1)t)}&{\alpha ((n + 2)t)}& \cdots &{\alpha (2nt)} \end{array}} \right] =\\ & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{C}}_h}} \\ {{{\boldsymbol{C}}_h}{{\boldsymbol{A}}_h}} \\ \vdots \\ {{{\boldsymbol{C}}_h}{{\boldsymbol{A}}_h}^{n - 1}} \end{array}} \right]{{\boldsymbol{A}}_h}\left[ {\begin{array}{*{20}{l}} {{{\boldsymbol{B}}_h}}&{{{\boldsymbol{A}}_h}{{\boldsymbol{B}}_h}}& \cdots &{{{\boldsymbol{A}}_h}^{n - 1}{{\boldsymbol{B}}_h}} \end{array}} \right] \approx \\ & {{\boldsymbol{U}}_1}\sqrt {{{\sum} _1}} {{\boldsymbol{A}}_h}\sqrt {{{\sum} _1}} {\boldsymbol{V}}_1^{\mathrm{T}}\quad. \\[-8pt] \end{split} $ (11)

By decomposing and representing the matrix A_h as:

$ {{\boldsymbol{A}}_h} \approx {\left(\sqrt {{{\sum} _1}} \right)^{ - 1}}{\boldsymbol{U}}_1^{\rm{T}}{{\boldsymbol{H}}_1}{{\boldsymbol{V}}_1}{\left(\sqrt {{{\sum} _1}} \right)^{ - 1}}\quad. $ (12)

Usually, in engineering systems, it is assumed that D_h = 0. By following the aforementioned steps, a mathematical model based on Hankel matrix identification method can be obtained. With the matrices A_h, B_h, C_h, and D_h calculated, the transfer function can be determined, establishing the mathematical model of the system and obtaining the dynamic linear part of the piezo-positioning system's Hammerstein model.

In accordance with the characteristics of the piezo-positioning system, an experimental platform as depicted in Fig.4 was set up. The entire experimental platform comprises a piezo-positioning actuator, an AD conversion module, a DA conversion module, a dSPACE simulation platform, a system controller, and an upper computer.

As deduced from the previous description, when a low-frequency voltage signal is applied to the piezo-positioning actuator, the output hysteresis curve changes very minimally, making it approximate to a static hysteresis system. Therefore, in this study, the static non-linear part of the Hammerstein model is identified at a frequency of 1 Hz, namely the P-I model. The accuracy of the P-I model is closely related to the number of Play operators. To balance the model accuracy and computational speed, 20 Play operators were chosen in the experiment to fit the static hysteresis curve. The fitting results of the threshold coefficient r and the weight coefficient ω are as follows:

$ \left\{ \begin{gathered} r = {\left[ {\begin{array}{*{20}{c}} 0&{0.15} &{0.30}&{\cdots}&{2.85} \end{array}} \right]^{\rm{T}}} \\ \omega = {\left[ {\begin{array}{*{20}{c}} {0.559\;7} & {0.137\;4} & {0.018\;0} & {\cdots} & { - 0.005\;8} \end{array}} \right]^{\rm{T}}} \\ \end{gathered} \right.. $ (13)

The model accuracy is represented by the relative error (RE) and the root mean square error (RMSE), defined as:

$ {\mathrm{RE}} = \sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^N {{{(\hat y(i) - y(i))}^2}} }}{{\displaystyle\sum\limits_{i = 1}^N {{{(y(i))}^2}} }}}\quad , $ (14)

$ {\mathrm{RMSE}} = \sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^N {{{(\hat y(i) - y(i))}^2}} }}{N}} \quad, $ (15)

where $ \hat y(i) $ is the model output and $ y(i) $ is the actual system output.

When the input sinusoidal signal has a frequency of 1 Hz, the hysteresis curve output by the P-I model matches the hysteresis curve measured from the piezo-positioning system, as shown in Fig.5 (color online). The modeling errors, RMSE and RE, are 0.0051 μm and 0.0034 μm, respectively, indicating excellent modeling results for the static nonlinear part.

Then, proceed with the identification of the dynamic linear part of the Hammerstein model. Once the P-I model is obtained, its inverse model can be calculated (the method for solving the inverse model is detailed in Section 4.1), then apply the derived P-I inverse model to the system's feedforward channel, and further identify the entire system using the Hankel matrix system identification method, thereby obtaining the mathematical model of the dynamic linear part.

When collecting input-output data of the system, the input voltage signal is a pseudo-random signal with an amplitude of 5 V and a length of 14 bytes. Based on the input-output data, the system model is identified. The input-output data obtained from the experiment is shown in Fig.6 (color online).

The system's impulse response can be obtained based on the input-output data of the system, then a Hankel matrix is constructed based on the impulse response, and its singular value decomposition is performed. The resulting singular value distribution is shown in Fig.7.

Upon observing Fig.7, it is evident that there is a significant jump between the first 6 data points, and from the 7th data point onwards, the singular values almost approach zero. After repeated experimental comparisons, the optimal choice for the system order is determined to be 6th order. Based on the determined system order and the previous outlined steps for Hankel matrix identification, for the system of 6th order, the identified transfer function of the dynamic linear part is as follows:

$ G(s) = \frac{{3.860\;8(s + 3.03 \times {{10}^4})({s^2} + 213s + 3.24 \times {{10}^6})}}{{(s + 2\;001)(s + 628.4)({s^2} + 1\;884s + 3.489 \times {{10}^6})}} \cdot \frac{{({s^2} - 1.383 \times {{10}^4}s + 1.141 \times {{10}^8})}}{{({s^2} + 138.8s + 9.957 \times {{10}^6})}}\quad. $ (16)

Firstly, the accuracy of the identified model is verified in the frequency domain. The comparison graph of the numerical solution obtained by the impulse response method with the frequency response of the identified system, as shown in Fig.8 (color online), reveals a good fitting effect before the frequency of 1000 Hz. This validates the accuracy of the identification results obtained by the Hankel matrix system identification method.

Further validating the accuracy of the model by inputting specific frequency sinusoidal signals ranging from 1 to 200 Hz into both the actual system and the identified system, and examining the model in the time domain. Tab.1 shows the errors between the output of the Hammerstein model and the actual output of the piezo-positioning system. By observing the data in the table, it can be noticed that the error between the model output and the system output is very small, thus further confirming the effectiveness of the constructed model.

To mitigate the impact of the hysteresis characteristic of the piezo-positioning actuator, a feedforward control approach can be employed. In this paper, the inverse model of the established P-I nonlinear model is derived, and this inverse model is used as a feedforward controller in series with the piezo-positioning actuator to compensate for the nonlinear error caused by the hysteresis characteristic. The structure of the feedforward compensation is illustrated in Fig.9. Here, u represents the input signal, x denotes the output of the feedforward controller, and y signifies the output signal.

The P-I model itself possesses the characteristic of analytic inverse. The inverse model of P-I, when connected in series with the piezo-positioning system model, can be constructed into a pseudo-linear functional relationship, thereby presenting a linearized relationship between the input and output of the entire system. The P-I inverse model also falls under the category of nonlinear models, with its parameters derivable through mathematical calculations based on the parameters of the P-I model. The organizational structure of the P-I inverse model remains a combination of a finite number of hysteresis operators and corresponding weight coefficients. Its expression is as follows:

$ \begin{split} & {u(t) = \sum\limits_{i = 1}^n {\omega _i' } \left\{ { \max \left\{ {y(t) - r_i' ,} \right.} \right.} \\ & {\left. {\left. {\min \left[ {y(t) + r_i' ,y\left( {{t_i}} \right)} \right]} \right\}} \right\}\quad,} \end{split} $ (17)

in the expression:

$ \omega _1' = \frac{1}{{{\omega _1}}}\quad, $ (18)

$ \omega _i' = \frac{{ - {\omega _i}}}{{\left(\displaystyle\sum\limits_{j = 1}^i {{\omega _j}} \right)\left(\displaystyle\sum\limits_{j = 1}^{i - 1} {{\omega _j}} \right)}}\quad, $ (19)

$ r_i' = \sum\limits_{j = 1}^i {{\omega _j}({r_i} - {r_j})} \quad, $ (20)

where ω_i' and r_i' represent the weight coefficients and threshold coefficients of the P-I inverse model.

The feedforward control structure obtained from the P-I inverse model belongs to open-loop control, with poor disturbance rejection capability and stability. In order to enhance the reliability of practical engineering applications, and considering the presence of external disturbances, a control method combining feedforward control with sliding mode closed-loop control will be further employed.

The transfer function of the dynamic linear part of the piezo-positioning system model is depicted as in equation (16), revealing the presence of unstable zeros in the right half-plane, indicating that the identified function is non-minimum phase. Designing a controller for non-minimum phase systems is notably challenging, particularly as the unstable zeros of this transfer function are significantly far from the imaginary axis. Therefore, this study employs an approximation approach by directly eliminating the unstable zeros while ensuring minimal alteration in system characteristics. By approximating the unstable zero in equation (16) as s approaching 0, the resulting transfer function becomes minimum phase:

$ {G_n}(s) = \frac{{4.405 \times {{10}^8}(s + 3.03 \times {{10}^4})}}{{(s + 2\;001)(s + 628.4)({s^2} + 1\;884s + 3.489 \times {{10}^6})}} \cdot \frac{{({s^2} + 213s + 3.24 \times {{10}^6})}}{{({s^2} + 138.8s + 9.957 \times {{10}^6})}}\quad. \\ $ (21)

We proceed with the subsequent design based on this.

From Fig.8, it can be observed that the system's frequency response curve exhibits prominent peak and valley. These can be mitigated by employing a notch filter for bandwidth adjustment, effectively eliminating the peak and valley in the high-frequency range, resulting in a smoother overall curve and thus enhancing the system's bandwidth. The principle behind the notch filter devised in this study is based on the positions of the zeros and poles in the frequency response plot (at 288 Hz and 502 Hz), allowing for the direct cancellation of the system's zeros and poles, thereby optimizing the model. The designed notch filter, denoted as G_f, is as follows:

$ {G_f}(s) = \frac{{3.24({s^2} + 138.8s + 9.957 \times {{10}^6})}}{{9.957({s^2} + 213s + 3.24 \times {{10}^6})}}. $ (22)

The transfer function of the system after performance adjustment by cascading the notch filter with the system is:

$ \begin{split} &{G_{nf}}(s) =\\ &\frac{{1.433 \times {{10}^8}(s + 3.03 \times {{10}^4})}}{{(s + 2\;001)(s + 628.4)({s^2} + 1\;884s + 3.489 \times {{10}^6})}}. \end{split}$ (23)

The performance adjustment effect of the model after incorporating the notch filter is illustrated in Fig.10 (color online). The notch filter not only optimizes the model but also achieves a reduction in the order of the piezo-positioning system's linear model, from the original 6th order to the 4th order, thereby facilitating the subsequent design of sliding mode controllers. Equation (23) is transformed into state space representation as:

$ \begin{gathered} {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} { - 4513}&{ - 9.7 \times {{10}^6}}&{ - 1.154 \times {{10}^{10}}}&{ - 4.387 \times {{10}^{12}}} \\ 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1&0 \end{array}} \right], \\ {\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \end{array}} \right]^{\rm{T}}}, \\ {\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} 0&0&{1.433 \times {{10}^8}}&{4.343 \times {{10}^{12}}} \end{array}} \right], \\ {\boldsymbol{D}} = \left[ 0 \right]. \\ \end{gathered} $ (24)

After the implementation of feedforward control, the hysteresis nonlinearity of the piezo-positioning actuator is effectively counteracted. However, two issues persist in this hysteresis compensation strategy: firstly, it is challenging to achieve complete compensation of the hysteresis nonlinearity in practice, and secondly, the piezo-positioning actuator itself exhibits parameter uncertainties and is susceptible to disturbances. Therefore, this paper proposes a control strategy based on P-I inverse compensation. By utilizing feedforward control of the inverse model to compensate for the hysteresis characteristics, and following model performance adjustment, a sliding mode controller is designed to suppress the remaining hysteresis characteristics, existing disturbances, and model uncertainties. As shown in Fig.11, the sliding mode inverse compensation control block diagram is depicted.

The state equation of the piezo-positioning actuator after feedforward compensation and performance adjustment is as follows:

$ \left\{ {\begin{split} & {\dot {\boldsymbol{x}} = {\boldsymbol{Ax}} + {\boldsymbol{Bu}}} \\ & {{\boldsymbol{y}} = {\boldsymbol{Cx}}} \end{split}} \right.\quad, $ (25)

where A is the n×n matrix, B is the n×m matrix, and C is the k×n matrix.

For the linear dynamic model of the system obtained through system identification and model performance adjustment, the state variables do not have actual physical significance. In other words, if the input voltage signal has actual physical meaning, it is not possible to directly control the state variables to achieve tracking performance. Therefore, according to the internal model principle, an integral component can be introduced to construct an integral augmented system for achieving tracking performance.

The system's output y(t) is required to track the reference input r(t), with the error signal e(t) given by:

$ e(t) = y(t) - r(t)\quad. $ (26)

Introducing the integral of the error signal, denoted as p(t):

$ p(t) = \int_0^{\rm{t}} e (\tau ){\mathrm{d}}\tau = \int_0^{\rm{t}} {\left[ {y(\tau ) - r(\tau )} \right]} {\mathrm{d}}\tau \quad. $ (27)

Taking the derivative of the above formula gives:

$ \dot p(t) = e(t) = Cx(t) - r(t)\quad. $ (28)

Connecting the integral component will increase the system order, hence considering the integral p as an augmented state variable, the state equation of the integral augmented system is formulated as follows:

$ \left\{ \begin{gathered} \left[ {\begin{array}{*{20}{c}} {\dot {\boldsymbol{x}}} \\ {\dot {\boldsymbol{p}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\boldsymbol{A}}&0 \\ {\boldsymbol{C}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\boldsymbol{x}} \\ {\boldsymbol{p}} \end{array}} \right] + \left[ {\begin{array}{*{20}{l}} {\boldsymbol{B}} \\ 0 \end{array}} \right]{\boldsymbol{u}} + \left[ {\begin{array}{*{20}{c}} 0 \\ { - {\boldsymbol{r}}} \end{array}} \right] \\ {\boldsymbol{y}} = \left[ {\begin{array}{*{20}{l}} {\boldsymbol{C}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\boldsymbol{x}} \\ {\boldsymbol{p}} \end{array}} \right] \\ \end{gathered} \right.\quad. $ (29)

The above state equation can also be expressed as:

$ \left\{ {\begin{split} & {{{\dot {\boldsymbol{x}}}_d} = {{\boldsymbol{A}}_d}{{\boldsymbol{x}}_d} + {{\boldsymbol{B}}_d}{\boldsymbol{u}} + {\boldsymbol{Er}}} \\ & {{\boldsymbol{y}} = {{\boldsymbol{C}}_d}{{\boldsymbol{x}}_d}} \end{split}} \quad.\right. $ (30)

After the augmentation of integrals, the system transitioned from its original n dimensions to n+k dimensions. Moreover, rank [B_dA_dB_dA_d²B_dA_d³B_dA_d⁴B_d] = 5 = n+k. According to the necessary and sufficient conditions for controllability of the system, it is evident that the augmented system is controllable. Therefore, it is possible to design a sliding mode controller for the augmented controlled system.

In designing a sliding mode controller, the initial step involves designing the sliding mode switching function, ensuring that it guarantees the asymptotic stability of the sliding motion and possesses excellent dynamic characteristics^[21]. Let the switching function be:

$ {\boldsymbol{S}} = {\boldsymbol{M}}{{\boldsymbol{x}}_d}\quad, $ (31)

where M is the k×(n+k) coefficient matrix.

Therefore, the sliding mode equation for the system during sliding motion on the switching surface is:

$ {\dot {\boldsymbol{x}}_d} = \left[ {{\boldsymbol{I}} - {{\boldsymbol{B}}_d}{{({\boldsymbol{M}}{{\boldsymbol{B}}_d})}^{ - 1}}{\boldsymbol{M}}} \right]{{\boldsymbol{A}}_d}{{\boldsymbol{x}}_d}\quad. $ (32)

Due to the rank of B_d being m, there exists a non-singular linear transformation x_d = T$ \tilde{{\boldsymbol{x}}} $, which transforms the augmented system state equation into the following form:

$ \left[ {\begin{array}{*{20}{c}} {{{\dot {\tilde {\boldsymbol{x}}}}_1}} \\ {{{\dot {\tilde {\boldsymbol{x}}}}_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\tilde {\boldsymbol{A}}}_{11}}}&{{{\tilde {\boldsymbol{A}}}_{12}}} \\ {{{\tilde {\boldsymbol{A}}}_{21}}}&{{{\tilde {\boldsymbol{A}}}_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde {\boldsymbol{x}}}_1}} \\ {{{\tilde {\boldsymbol{x}}}_2}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} 0 \\ {{{\tilde {\boldsymbol{B}}}_2}} \end{array}} \right]{\boldsymbol{u}}\quad, $ (33)

where, $ {\tilde {\boldsymbol{x}}_1} \in {R^{n + k - m}} $, $ {\tilde {\boldsymbol{x}}_2} \in {R^m} $, $ {\tilde {\boldsymbol{B}}_2} $ are m×m invertible matrices.

The corresponding switching surface becomes:

$ {\boldsymbol{S}} = {\boldsymbol{MT}}\tilde {\boldsymbol{x}} = {{\boldsymbol{M}}_1}{\tilde {\boldsymbol{x}}_1} + {{\boldsymbol{M}}_2}{\tilde {\boldsymbol{x}}_2} = 0\quad. $ (34)

Furthermore, based on equation (33), it can be obtained that:

$ {\dot {\tilde {\boldsymbol{x}}}_1} = {\tilde {\boldsymbol{A}}_{11}}{\tilde {\boldsymbol{x}}_1} - {\tilde {\boldsymbol{A}}_{12}}{\boldsymbol{M}}_2^{ - 1}{{\boldsymbol{M}}_1}{\tilde {\boldsymbol{x}}_1} = ({\tilde {\boldsymbol{A}}_{11}} - {\tilde {\boldsymbol{A}}_{12}}{\boldsymbol{F}}){\tilde {\boldsymbol{x}}_1}. $ (35)

The above equation is fully equivalent to the sliding mode equation, indicating that the sliding mode control system can be represented as an n+k-m dimensional subsystem described by the above equation. Moreover, since $ ({{\boldsymbol{A}}_d},{{\boldsymbol{B}}_d}) $ is controllable, it follows that $ ({\tilde {\boldsymbol{A}}_{11}},{\tilde {\boldsymbol{A}}_{12}}) $ is also controllable, thus allowing the poles of the aforementioned subsystem to be arbitrarily placed by F.

According to equation (34), the coefficient matrix M of the switching function is derived as:

$ {\boldsymbol{M }}= \left[ {\begin{array}{*{20}{l}} {\boldsymbol{F}}&{{{\boldsymbol{M}}_2}} \end{array}} \right]{{\boldsymbol{T}}^{ - 1}}\quad. $ (36)

Thus, the sliding mode switching function can be determined, ensuring that the sliding mode motion is asymptotically stable.

Next comes the design of the sliding mode control law. In order to perform sliding mode motion for the system, the control law must satisfy the reachability condition:

$ {\boldsymbol{S}}\dot {\boldsymbol{S}} \leqslant 0\quad. $ (37)

This study employs the hyperbolic tangent function tanh(αS) to replace the sign function sgn(S) typically used in exponential reaching laws, which effectively reduces chattering. The improved exponential reaching law is:

$ \dot {\boldsymbol{S}} = - \varepsilon \tanh (\alpha {\boldsymbol{S}}) - k{\boldsymbol{S}}\quad, $ (38)

where, ε > 0, α > 0, k > 0.

Subsequently, from the improved exponential reaching law, the sliding mode control law u satisfying the reachability condition can be derived as:

$ {\boldsymbol{u}} = - {({\boldsymbol{M}}{{\boldsymbol{B}}_d})^{ - 1}}({\boldsymbol{M}}{{\boldsymbol{A}}_d}{{\boldsymbol{x}}_d} + \varepsilon \tanh ({\boldsymbol{\alpha}} {\boldsymbol{S}}) + {\boldsymbol{kS}}). $ (39)

The motion of the sliding mode control system consists of two parts. The first part is the motion stage where the system enters the switching surface from the initial point. By designing a Lyapunov function V(x_d) = S², since the sliding mode control law u always satisfies the reachability condition, it follows that $ \dot V({x_d}) \leqslant 0 $. The second part is the motion stage of the system on the switching surface, where the system's stability can be ensured as long as the sliding mode switching function is appropriately designed. In conclusion, it can be seen that the sliding mode control system designed in this study is asymptotically stable.

To achieve sliding mode inverse compensation control of the piezo-positioning system, it is necessary to obtain the full-dimensional state variable x_d of the integral augmented system. The state variable p in x_d represents the integral of the error signal, which can be directly measured, while the state variable x does not have a physical meaning and thus cannot be directly measured. In order to obtain the state variable x, it is essential to design a state observer for the original system.

For the original system (25), rank [C^TA^TC^T (A^T)²C^T (A^T)³C^T] = 4, it can be concluded from the necessary and sufficient conditions for observability that the original system is observable. The equation for the designed state observer is as follows:

$ \dot {\hat {\boldsymbol{x}}} = \left( {{\boldsymbol{A}} - {\boldsymbol{GC}}} \right)\hat {\boldsymbol{x}} + {\boldsymbol{Bu}} + {\boldsymbol{Gy}}\quad, $ (40)

where G represents the feedback gain matrix.

Constructing an asymptotic switching function based on the estimated state $ \hat{{\boldsymbol{x}}} $:

$ \hat {\boldsymbol{S}} = {\boldsymbol{M}}{\hat {\boldsymbol{x}}_d}\quad. $ (41)

Simultaneously, the sliding mode control law of the system is also transformed to:

$ {\boldsymbol{u}} = - {({\boldsymbol{M}}{{\boldsymbol{B}}_d})^{ - 1}}({\boldsymbol{M}}{{\boldsymbol{A}}_d}{\hat {\boldsymbol{x}}_d} + \varepsilon \tanh ({\boldsymbol{\alpha \hat S}}) + {\boldsymbol{k}}\hat {\boldsymbol{S}}). $ (42)

At this point, the design of the piezo-positioning actuator integral augmented sliding mode inverse compensation control system is complete, with its system structural diagram shown in Fig.12.

To verify the hysteresis compensation effect of the P-I inverse model, an open-loop inverse compensation control experiment is first conducted. The input is a sinusoidal signal with a frequency of 1 Hz and an amplitude of 15 V. The obtained fitting effect of the compensated input and output is shown in Fig.13 (color online), where it can be observed that after the introduction of the feedforward controller, the input-output curves of the system basically overlap, indicating a good compensation for the hysteresis nonlinearity. It is also evident that hysteresis nonlinearity is highly intricate, making it impossible to achieve extremely precise modeling in practice and completely compensate for all hysteresis characteristics, thus necessitating the use of sliding mode control to further suppress the incompletely compensated hysteresis nonlinearity.

Building upon the open-loop inverse compensation control, the sliding mode closed-loop control is then implemented based on the previously designed sliding mode controller. In order to validate the effectiveness and superiority of sliding mode inverse compensation control, this study designs a comparative experiment involving PID inverse compensation control and sliding mode control without inverse compensation.

During the experiment, through repeated parameter simulations and physical tuning, taking into account the system's stability and robustness, parameter settings and adjustments were carried out. The finalized pole configuration of the sliding mode control equivalent subsystem (36) was determined to be [−3000 ± 200i, −2500 ± 150i], with its corresponding state feedback matrix F = [5.67×10¹³, 8.28×10¹⁰, 4.53×10⁷, 1.10×10⁴]. Consequently, the switching function coefficient matrix M = [1.46, 1.35×10⁴, 5.15×10⁷, 8.91×10¹⁰, 14.07]. The pole configuration of the state observer was set as [−5000 ± 500i, −3000 ± 200i], with its corresponding feedback gain matrix G = [11.76, −0.01, 1.90×10⁻⁶, 1.83×10⁻⁹]^T. The parameters for the improved exponential reaching law were set as ε = 110, α = 2, k = 1200, under which better control performance can be achieved.

For the PID inverse compensation control as part of the comparative experiment, considering the system control performance in both the time domain and frequency domain, the parameters of the PID controller were ultimately set as K_p = 0.42, K_i = 344, K_d = 0.

Selecting a step signal as the system input to test the control effect of sliding mode inverse compensation control in the time domain, with a step signal amplitude of 10 V and an expected output displacement of 10 μm. The step responses of the three control schemes are shown in Fig.14 (color online).

The system under PID inverse control exhibits an overshoot of 10.5% in its step response, with a settling time of 13.9 ms, using a ±2% error band as the criterion. On the other hand, the sliding mode control without inverse compensation shows no overshoot and a settling time of 9.0 ms. Similarly, the sliding mode inverse control also demonstrates no overshoot and a settling time of 6.2 ms. It is evident that compared to PID inverse control, both sliding mode inverse control and sliding mode control exhibit a significant reduction in overshoot. The settling time of sliding mode inverse control is approximately 55.4% shorter than PID inverse control and 31.1% shorter than sliding mode control. Furthermore, the sliding mode inverse control shows almost no sign of oscillation, while the sliding mode control exhibits slight oscillations due to the system's non-linear characteristics.

To assess the tracking control performance in the frequency domain, a sinusoidal sweep signal with frequencies ranging from 0.1 to 500 Hz, an amplitude of 15 V, and a bias of 15 V is used as the system input. The obtained closed-loop frequency characteristic curves for the three control schemes are depicted in Fig.15 (color online).

Based on the experimental results, the closed-loop tracking bandwidth (defined at −3 dB) for the system under PID inverse control is 91.5 Hz, for sliding mode control is 110.1 Hz, and for sliding mode inverse control is 119.9 Hz. Comparatively, sliding mode inverse control shows an increase in closed-loop tracking bandwidth of approximately 8.9% and 31.0% when compared to sliding mode control and PID inverse control, respectively. It is evident that sliding mode inverse compensation control demonstrates significant advantages, showcasing a more stable tracking control performance under varying frequency sweep input signal.

To evaluate the disturbance rejection control performance in the frequency domain, a sinusoidal sweep signal with frequencies ranging from 0.1 to 500 Hz, an amplitude and bias of 1.5 V, is introduced as the disturbance signal. The system input is set to 0. The disturbance rejection magnitude-frequency characteristic curves for the three control schemes are illustrated in Fig.16 (color online).

After testing, it is determined that the disturbance rejection bandwidth (defined at 0 dB) for the system under PID inverse control is 59.9 Hz, for sliding mode control is 63.9 Hz, and for sliding mode inverse control is 86.2 Hz. Comparatively, sliding mode inverse control shows an increase in disturbance rejection bandwidth of approximately 34.9% and 43.9% when compared to sliding mode control and PID inverse control, respectively. It is evident that sliding mode inverse compensation control exhibits a more effective disturbance rejection performance.

This study examines the system complexity characteristics of piezo-positioning systems and their demanding control bandwidth. Our research endeavors to develop a precise Hammerstein model comprising a P-I model and a dynamic linear model in series. Building upon this, we propose a sliding mode inverse compensation control method that integrates a P-I inverse model with integral augmentation.

The experimental results demonstrate that the dynamic hysteresis model of the piezo-positioning system developed in this study offers excellent generalization across typical input frequencies below 200 Hz. The implemented sliding mode inverse compensation control effectively mitigates hysteresis nonlinearity, significantly enhancing the system's control precision and stability. Compared to PID inverse control and sliding mode control without inverse compensation, sliding mode inverse compensation control achieves a step response free of overshoot and oscillations, with a reduced settling time of 6.2 ms. This method exhibits robust resistance to hysteresis nonlinearity and external disturbances, showing enhanced adaptability to various input signals. The control system attains a closed-loop tracking bandwidth of 119.9 Hz and a disturbance rejection bandwidth of 86.2 Hz. This study validates the effectiveness and superiority of the sliding mode inverse compensation control approach in improving tracking accuracy and disturbance rejection capabilities of piezo-positioning systems in practical engineering applications.