LSTM-Inversion-Based Feedforward–Feedback Nanopositioning Control

Abstract

This work proposes a two-degree of freedom (2DOF) controller for motion tracking of nanopositioning devices, such as piezoelectric actuators (PEAs), with a broad bandwidth and high precision. The proposed 2DOF controller consists of an inversion feedforward controller and a real-time feedback controller. The feedforward controller, a sequence-to-sequence LSTM-based inversion model (invLSTMs2s), is used to compensate for the nonlinearity of the PEA, especially at high frequencies, and is collaboratively integrated with a linear MPC feedback controller, which ensures the PEA position tracking performance at low frequencies. Therefore, the proposed 2DOF controller, namely, invLSTMs2s+MPC, is able to achieve high precision over a broad bandwidth. To validate the proposed controller, the uncertainty of invLSTMs2s is checked such that the integration of an inversion model-based feedforward controller has a positive impact on the trajectory tracking performance compared to feedback control only. Experimental validation on a commercial PEA and comparison with existing approaches demonstrate that high tracking accuracies can be achieved by invLSTMs2s+MPC for various reference trajectories. Moreover, invLSTMs2s+MPC is further demonstrated on a multi-dimensional PEA platform for simultaneous multi-direction positioning control.

1. Introduction

Nanopositioning devices, such as piezoelectric actuators (PEAs), are broadly used in many high-precision industries, products, and systems due to their advantages of fast response and high motion resolution. As one of the most popular nanopositioning devices, PEAs have been implemented in various applications, such as atomic force microscopes (AFMs) [1,2], micro forming [3], and adaptive optics [4]. However, the operation of such a system at high frequency and/or large motion range still remains challenging because of the PEA’s nonlinear dynamics (e.g., hysteresis and creep), which is more pronounced when operations are performed under these conditions. As a result, the system operation bandwidth and motion range are rather limited for commercial applications, especially for real-time motion control applications.

To address these issues, efforts in PEA position tracking control have been reported. Feedforward controllers have been proposed to increase the control bandwidth by compensating for the PEA nonlinear dynamics. For example, model-based feedforward control designs have been implemented to improve control performance by compensating for PEA hysteresis, through inverse dynamics approximation [5,6]. To capture the frequency-dependent PEA hysteresis dynamics, rate-independent hysteresis models, such as the improved Preisach model [7], the modified Prandtl–Ishlinskii model [8], and the Bouc–Wen hysteresis model [9], have been implemented in the feedforward controller designs. The feedforward controller can then be combined with a real-time feedback controller to form a two-degree-of-freedom (2DOF) controller to increase the overall tracking accuracy and robustness [8,10]. In these 2DOF controllers, the feedforward part accounts for the nonlinear and high-frequency PEA dynamics, and the feedback controller takes care of the remaining errors to improve the overall control precision. Theoretically, the 2DOF controller should be able to deliver satisfactory trajectory tracking and positioning performance if both feedforward and feedback parts are well designed. However, due to the complicated PEA nonlinear dynamics, model parameter identification is quite challenging, and the trajectory tracking improvements obtained are rather limited [11], especially when system uncertainties exist.

With the recent development of deep learning algorithms, neural networks (NNs) have been widely proposed for system identification and controller design. With adequate training data, NNs can capture the entire PEA dynamics presented in the data quite well, not just limited to hysteresis. Previous studies [12] show a vanilla recurrent neural network (RNN) performs well at handling the time series input and is increasingly popular in controller design. However, a vanilla RNN still has its limits. For example, an RNN can potentially suffer from gradient vanishing or exploding issues [13,14,15] and is prone to overfitting. Furthermore, compared to the long short-term memory (LSTM), it cannot consider and understand the system’s long-term dependencies subject to sequential input, so the accuracy of a vanilla RNN in system dynamic identification can be limited [16,17]. The advantages of LSTM come from the memory cells in network neurons. The memory cells inherit the time-relevant information from long-term and short-term memory together and help LSTM make future predictions based on important historical data [16,18]. This feature enables LSTM to perform better in handling the complexity of nonlinear systems. But LSTM still can possibly encounter overfitting issues [19,20,21]. The reason can be originated from the conventional training mechanism. To train an accurate LSTM, it conventionally uses a very long single time series that concatenates many training data samples to cover a broad frequency range [16,18]. As a result, the conventional LSTM training progress could be extremely time-consuming. To address this issue, sequence-to-sequence LSTM learning (namely, LSTMs2s) has been recently implemented for PEA system identification [22] with high accuracy. This provides potential for controller design in high accuracy PEA position tracking.

Leveraging LSTMs2s in system identification, this work proposes a 2DOF controller (invLSTMs2s+MPC) consisting of a feedforward controller based on an LSTMs2s model of the system inverse dynamics (namely, invLSTMs2s) and a model predictive controller (MPC) as a feedback controller. The MPC was chosen primarily owing to its ability to handle constraints and capability for controlling multivariable plants, among other advantages. Considering the frequency and range dependent PEA nonlinear dynamics, the invLSTMs2s training set was generated using sinusoidal and triangle references with various frequency and amplitude components. To reduce the training data size and training time cost, invLSTMs2s was primarily designed to capture the PEA dynamics at mid- and high-frequency ranges (e.g., a range of tens to hundreds of hertz). In this work, the invLSTMs2s model was trained using a parallel training method, in which groups of training datasets were divided into mini-batches. After training, the modeling accuracy of invLSTMs2s was validated by comparing the model hysteresis prediction with the actual PEA response. With the accuracy of invLSTMs2s in nonlinear dynamics identification confirmed, an MPC based on a linear model was integrated as a feedback controller to take care of the low-frequency control and improve the overall tracking accuracy. Weighting factors were incorporated in the feedforward and feedback control paths to generalize the proposed 2DOF controller for tracking various references in different frequency regions. For validation and demonstration, the proposed 2DOF controller was implemented to control the displacement of a commercial PEA stage and compared with other approaches. The experiment results clearly demonstrate the superiority of invLSTMs2s+MPC for PEA trajectory tracking.

2. Feedforward–Feedback invLSTMs2s+MPC Controller Design

2.1. Feedforward Controller

2.1.1. System Inversion Model

In general, a dynamic system can be modeled in state-space as Equation (1), where $U=[u_k,u_{k+1},\dots \dots ]$ and $Y=[y_k,y_{k+1},\dots \dots ]$ are the input and output, respectively. The corresponding inversion model can be modeled as Equation (2). The inversion model’s (Equation (2)) output becomes U when the Y is fed as the inversion model’s input. If an accurate inversion model can be obtained for a time-invariant system and disturbances/uncertainties are not considered, precise control of the system can be achieved by connecting the inversion model and the system in series where the output of the former is the input of the latter and directly feeding the desired system reference as the input of the combination.

$$\begin{gathered} x_{k+1}=f(x_k,u_k) \\ {y}_k=h(x_k,u_k) \end{gathered}$$

(1)

$$\begin{gathered} z_{k+1}=f_{inv}(z_k,y_k) \\ {u}_k=h_{inv}(z_k,y_k). \end{gathered}$$

(2)

For a linear system, its inversion model can be directly determined through linear system theories. However, for a nonlinear system, such as PEA positioning systems in this work, obtaining an accurate inversion model is not straightforward. One way to solve this problem is to train and regress an inversion machine learning model using the data pair $\{Y,U\}$ from the original system.

2.1.2. LSTMs2s Inversion Model (invLSTMs2s)

The goal of designing an inversion model is to map the output time series of the original system back to the input time series of the original system. Previous studies have shown that LSTM networks are more potent for this task compared to other forms of RNNs [17,22].

In this work, an advanced development of LSTM, sequence-to-sequence LSTM (LSTMs2s), was selected to generate the nonlinear inversion model of PEA for its efficiency and accuracy [22,23,24]. As shown in Figure 1, the LSTMs2s contains two LSTM layers expanded timewise: an LSTM encoder and an LSTM decoder. Both LSTM layers have similar structures.

In LSTMs2s (see Figure 1), the encoder LSTM layer reads the entire input time series $u_{\left(r\right)}$ and its final hidden states $x_{\left(r\right)}$ are saved as context vector $C_v$. After the encoder processes the entire input time series, the decoder is coupled with a fully connected network to predict the output time series $y_{\left(r\right)}$ based on $C_v$. The entire LSTMs2s model can be formulated as follows:

Encoder:

$$x_{e,\left(k\right)}=f(x_{e,(k-1)},u_{\left(r\right),k}).$$

(3)

Decoder:

$$\begin{gathered} y_{\left(r\right),k-1}=g\left(x_{d,(k-1)}\right) \\ {x}_{d,\left(k\right)}=f(x_{d,(k-1)},y_{\left(r\right),k-1},C_v) \\ {y}_{\left(r\right),k}=g\left(x_{d,\left(k\right)}\right), \end{gathered}$$

(4)

where $x_{d,\left(k\right)}$ and $x_{e,\left(k\right)}$ are the hidden states in the encoder and decoder at the sampling instant k, respectively. $f(·)$ denotes the LSTM unit function. $g(·)$ represents the output layer that is implemented by a fully connected network. In this encoder–decoder model, it is clear that the decoder generates output data points $y_{\left(r\right),k}$ at each sampling instant k and uses it as an input parameter to the decoder for the next output prediction. Note that the context vector $C_v$ serves as an input parameter to the decoder for each step.

2.1.3. invLSTMs2s Training

invLSTMs2s aims to map the relation from the PEA output to the PEA input, hence the inverse PEA dynamics. Suppose the measured PEA system output is $Y_{\left(ts\right)}$ subject to the drive input $U_{\left(ts\right)}$, and $Y_{\left(rts\right)}$ is the inversion model output subject to $Y_{\left(ts\right)}$; a perfect inversion model should generate $Y_{\left(rts\right)}$ exactly the same as $U_{\left(ts\right)}$. In other words, such a perfect inversion model satisfies that the error $\left|\right|U_{\left(ts\right)}-Y_{\left(rts\right)}\left|\right|<ϵ$ for any $ϵ>0$. Therefore, the invLSTMs2s can be trained using the training set constructed with $(Y_{ts},U_{ts})$ by solving the following optimization problem:

$$\begin{array}{cc}\underset{W_f,b_f,W_i,b_i,W_o,b_o,W_C,b_C,W_{\left(r\right)},b_{\left(r\right)}}{min}& J_{\left(r\right)}=\left|\right|U_{\left(ts\right)}-Y_{\left(rts\right)}\left|\right|\hfill \\ \hfill \mathrm{subject}\mathrm{to}:& x_{e,\left(k\right)}=f(x_{e,(k-1)},y_{\left(ts\right),k})\hfill \\ \hfill & y_{\left(rts\right),k-1}=g\left(x_{d,(k-1)}\right)\hfill \\ \hfill & x_{d,\left(k\right)}=f(x_{d,(k-1)},y_{\left(r\right),k-1},s)\hfill \\ \hfill & y_{\left(rts\right),k}=g\left(x_{d,\left(k\right)}\right)\hfill \\ \hfill & x_{d,0}=x_{e,L}=s,k=1,2,3,\dots \dots ,L\hfill \end{array},$$

(5)

where $f(·)$ and $g(·)$ are the LSTM unit function and the fully connected output layer, respectively. L is the length of the time series.

We utilized parallel training to train this structure with multiple input–output time series samples for highly efficient system identification. Considering the frequency- and amplitude-dependent behavior of PEAs, each time series $U_{\left(ts\right),i}$ is generated using the multi-period sinusoidal function $S(f_i,A_i)$ and triangular function $T(f_i,A_i)$, with $f_i$ as the signal frequency and $A_i$ as the amplitude. Mathematically, the collection of k samples of $U_{\left(ts\right),i}$, $U_{\left(ts\right)}$ is formulated as

$$\begin{array}{cc}\hfill & U_{\left(ts\right)}=\{U_{\left(ts\right),1},U_{\left(ts\right),2},\dots ,U_{\left(ts\right),k}\}\hfill \\ \hfill & \mathrm{with}{U}_{\left(ts\right),i}=\left\{\begin{array}{cc}{S}_i,i={\displaystyle \frac{k+1}{2}}\hfill & \mathrm{if}k\mathrm{is}\mathrm{odd}\mathrm{number}\hfill \\ T_i,i={\displaystyle \frac{k}{2}}\hfill & \mathrm{if}k\mathrm{is}\mathrm{even}\mathrm{number}\hfill \end{array}\right.\hfill \\ \hfill & S_i=S(f_i,A_i),T_i=T(f_i,A_i),\mathrm{where}(f_i,A_i)\in \mathrm{\Omega }\hfill \\ \hfill & S(f_i,A_i)=A_i[sin(2\pi f_{i}t+{\displaystyle \frac{3\pi }{2}})+1],t\in [0,{\displaystyle \frac{P}{f}}]\hfill \\ \hfill & T(f_i,A_i)=4·A_i|f_{i}t-⌊f_{i}t+{\displaystyle \frac{1}{2}}⌋|,t\in [0,{\displaystyle \frac{P}{f}}],\hfill \end{array}$$

(6)

where $\mathrm{\Omega }$ is the set of selected $(f_i,A_i)$ pairs, which can be determined using k-means clustering [16,25]. P is the number of periods the signal sample contains. The PEA system output time series set $Y_{\left(ts\right)}=\{Y_{\left(ts\right),1},Y_{\left(ts\right),2},\dots ,Y_{\left(ts\right),k}\}$ can be then measured subject to the input set $U_{\left(ts\right)}=\{U_{\left(ts\right),1},U_{\left(ts\right),2},\dots ,U_{\left(ts\right),k}\}$. Then, the entire training set can be constructed as $\{Y_{\left(ts\right)},U_{\left(ts\right)}\}$. The samples of the entire training set are then sorted based on their length and divided into mini-batches using a pre-defined batch size n (n samples per batch). Then, the inversion model is trained batch by batch, and the total training iteration needed is $\lambda K_s/n$, where $\lambda$ is the epoch number, and $K_s$ is the total sample size.

Note that the sample length can be very big for small $f_{i}s$. This can dramatically increase the training time. Given the PEA’s nonlinear dynamics mostly occur for high driving frequencies [22,25], the invLSTMs2s was primarily designed to capture the high-frequency and nonlinear dynamics of the PEA system [22].

2.2. Feedback Controller Design

Given the high-frequency dynamics of the PEA can be accurately compensated by the invLSTMs2s, what is left for the PEA positioning control is to design a controller for the PEA’s low-frequency dynamics, which is mostly linear. Therefore, a linear MPC is proposed for its efficiency and flexibility in handling control constraints and system disturbances. The MPC works in parallel with the invLSTMs2s feedforward controller. The block diagram of the entire controller design is shown in Figure 2.

The MPC utilizes a state observer and an estimated linear model for the PEA system state estimation and output prediction. Such a state observer, like the Kalman filter, can be expressed as,

$$\begin{array}{cc}\hfill x_{\left(c\right),k+1}& =A{x}_{\left(c\right),k}+B{u}_k\hfill \\ \hfill \widehat{y_k}& =C{x}_{\left(c\right),k},\hfill \end{array}$$

(7)

where $x_{\left(c\right),k}$ is the estimated state of the state observer, and $u_k$ is the system input. $\widehat{y_k}$ is the predicted system output. $A,B,C$ are the estimated state-space matrices that can be obtained using a linear system identification algorithm based on system frequency response fitting, and their dimensions depend on the order of the estimated system chosen by the algorithm and user preference.

State estimation: A steady-state Kalman filter [26] was implemented as the state estimator for the MPC. To design the Kalman filter, the process noise covariance and measurement noise covariance were estimated based on the system’s linear model identification accuracy, the PEA stage’s measurement noise level, and the manufacturer’s provided system’s noise-to-signal ratio.

For each control interval k, the system state was estimated as follows:

First, compute the innovation using the current system output measurement $y_k$

$$e_k=y_k-\left[C{x}_{\left(c\right),k|k-1}\right].$$

(8)

Second, update the current state

$$x_{\left(c\right),k|k}=x_{\left(c\right),k|k-1}+M{e}_k,$$

(9)

where $x_{\left(c\right),k|k-1}$ is the k instant controller state estimation at the previous instant, $k-1$. M is the Kalman innovation gain matrix. Then, the MPC uses the updated state $x_{\left(c\right),k|k}$ to predict the system output and to solve for the optimal drive input to the system, $u_k^{opt}$, with quadratic programming (QP).

Finally, the state estimator predicts the state for the next control interval $k+1$

$$x_{\left(c\right),k+1|k}=A{x}_{\left(c\right),k|k-1}+B{u}_k^{opt}+L{e}_k,$$

(10)

where L is the Kalman filter gain matrix.

Output prediction: The MPC computes the future system state within the pre-selected prediction horizon p as follows:

$$x_{\left(c\right),k+i|k}=A{x}_{\left(c\right),k+i-1|k}+B{u}_{k+i-1|k},i=1,2,3,\dots ,p.$$

(11)

For each step i, the PEA’s output is:

$${\widehat{y}}_{k+i}=C{x}_{\left(c\right),k+i|k}.$$

(12)

QP matrices: Let $U=[u_{k+1},u_{k+1},\dots ,u_{k+m}]$, $\mathrm{\Delta }U=[u_{k+1}-u_k,u_{k+2}-u_{k+1},\dots ,u_{k+m}-u_{k+m-1}]$, and $1_n={[1,1,\dots ,1]}^T$, where m is the control horizon. The aforementioned system output prediction can be rewritten in an alternative form:

$$\begin{array}{cc}\hfill \widehat{Y}& =\alpha x_{\left(c\right),k|k}+\beta U+F{u}_k\hfill \\ \hfill & =\alpha x_{\left(c\right),k|k}+\beta (S_1\mathrm{\Delta }U+1_n{u}_k)+F{u}_k\hfill \\ \hfill & =\alpha x_{\left(c\right),k|k}+\beta S_1\mathrm{\Delta }U+(\beta 1_n+F)u_k,\hfill \end{array}$$

(13)

where

$$\begin{array}{c}\hfill \widehat{Y}={\left[\begin{array}{c}{\widehat{y}}_{k+1}\\ {\widehat{y}}_{k+2}\\ ⋮\\ {\widehat{y}}_{k+p}\end{array}\right]}_{p\times 1},U={\left[\begin{array}{c}{u}_{k+1}\\ u_{k+2}\\ ⋮\\ u_{k+m}\end{array}\right]}_{m\times 1}.\end{array}$$

The terms $\alpha$ and $\beta$ are defined as

$$\begin{array}{cc}\hfill \alpha ={\left[\begin{array}{c}CA\\ C{A}^2\\ ⋮\\ C{A}^p\end{array}\right]}_{p\times 1},& \beta ={\left[\begin{array}{ccccc}0& 0& 0& \cdots & 0\\ CB& 0& 0& \cdots & 0\\ CAB& CB& 0& ·& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C{A}^{p-2}B& C{A}^{p-3}B& C{A}^{p-4}B& ·& 0\end{array}\right]}_{p\times m},\hfill \\ \hfill F=\left[\begin{array}{c}CB\\ CAB\\ C{A}^{2}B\\ ⋮\\ C{A}^{p-1}B\end{array}\right],& S_1=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 1& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 1\end{array}\right].\hfill \end{array}$$

Cost function: With the given QP matrices, the cost function can be mathematically expressed as

$$J={(\widehat{Y}-R)}^T(\widehat{Y}-R)+\rho \mathrm{\Delta }{U}^T\mathrm{\Delta }U,$$

(14)

where $R=[r_{k+1},r_{k+2},\dots ,r_{k+p}]$ is the reference of desired system output, $\rho$ is the weighting coefficient for the reference tracking error and MV movement suppression, respectively. Let $E=\alpha x_{\left(c\right),k|k}+(\beta 1_n+F)u_k-R$; the cost function can be rewritten as

$$\begin{array}{cc}\hfill J& ={(\beta S_1\mathrm{\Delta }U+E)}^T(\beta S_1\mathrm{\Delta }U+E)+\mathrm{\Delta }{U}^T\left(\rho I\right)\mathrm{\Delta }U\hfill \\ \hfill & =\mathrm{\Delta }{U}^T{S}_1^T{\beta }^T\beta S_1\mathrm{\Delta }U+E^T\beta S_1\mathrm{\Delta }U+{\left(E^T\beta S_1\mathrm{\Delta }U\right)}^T+E^{T}E+\mathrm{\Delta }{U}^T\left(\rho I\right)\mathrm{\Delta }U\hfill \\ \hfill & =\mathrm{\Delta }{U}^T(\rho I+S_1^T{\beta }^T\beta S_1)\mathrm{\Delta }U+2E^T\beta S_1\mathrm{\Delta }U+E^{T}E.\hfill \end{array}$$

(15)

Now, it becomes a convex optimization problem for the QP progress to compute the optimal $\mathrm{\Delta }U$ based on the known value $u_k$ and $x_{\left(c\right),k|k}$, and add the first element of $\mathrm{\Delta }U$ with $u_k$ to obtain $u_k^{opt}$. The convex optimization problem can be formulated as follows:

$$\begin{array}{cc}\hfill \underset{\mathrm{\Delta }U}{argmin}& J=\mathrm{\Delta }{U}^T(\rho I+S_1^T{\beta }^T\beta S_1)\mathrm{\Delta }U+2E^T\beta S_1\mathrm{\Delta }U+E^{T}E.\hfill \end{array}$$

(16)

Solution and stability analysis: The solution for $\mathrm{\Delta }U$ that can minimize the cost function J should be

$$\mathrm{\Delta }U=arg\underset{\mathrm{\Delta }U}{min}J=-{(\rho I+S_i^T{\beta }^T\beta S_i)}^{-1}{E}^T\beta S_i.$$

(17)

According to the MPC algorithm, only the first step of the manipulated variable would be taken for system input. Setting $G=[1,0,0,\dots ,0]$, we can compute

$$\begin{array}{cc}\hfill u_{k+1}& =u_k-G{(\rho I+S_i^T{\beta }^T\beta S_i)}^{-1}{E}^T\beta S_i\hfill \\ \hfill & =M_1{x}_c\left(k\right)+M_2{u}_k+M_{3}R,\hfill \end{array}$$

(18)

where $M_1=-G{(\rho I+S_i^T{\beta }^T\beta S_i)}^{-1}{\alpha }^T\beta S_i$, $M_2=1-G{(\rho I+S_i^T{\beta }^T\beta S_i)}^{-1}{\left(\beta 1_n\right)}^T\beta S_i$, and $M_3=G{(\rho I+S_i^T{\beta }^T\beta S_i)}^{-1}\beta S_i$. In this case, $M_1$, $M_2$, and $M_3$ are all constant matrices and can be regarded as the coefficient to $x_{\left(c\right),k|k}$, $u_{k-1}$, and R, so that the feedback system can be rewritten as

$$\begin{array}{c}\hfill \left[\begin{array}{c}{x}_{\left(c\right),k+1|k}\\ u_{k+1}\end{array}\right]=\left[\begin{array}{cc}A& B\\ M_1& M_2\end{array}\right]\left[\begin{array}{c}{x}_{\left(c\right),k|k}\\ u_k\end{array}\right]+\left[\begin{array}{c}0\\ M_3\end{array}\right]R.\end{array}$$

(19)

Based on the alternative form above, it is obvious that the feedback system stability condition needs to satisfy the following:

$${\left|\lambda \left(\left[\begin{array}{cc}A& B\\ M_1& M_2\end{array}\right]\right)\right|}_{max}<1,$$

(20)

where $\lambda (·)$ denotes the eigenvalues of the matrix.

2.3. Impact of the invLSTMs2s Uncertainty

Previous studies have shown that the inversion’s modeling uncertainty must be small such that the inversion-based feedforward input improves the output tracking performance when compared to the case with feedback control only (regardless of the type of feedback controller used) [27]. Thus, it is necessary to evaluate the uncertainty of invLSTMs2s in this work. According to [27], for the SISO PEA system, the criteria that need to be satisfied such that integrating the inversion feedforward with the MPC improves the tracking performance can be formulated as

$${\left|\mathrm{\Delta }\left(j\omega \right)\right|}_2\le {\left|G_o\left(j\omega \right)\right|}_2,$$

(21)

where $\mathrm{\Delta }\left(j\omega \right)$ is the modeling uncertainty at different frequencies $\omega$, and it is defined as $\mathrm{\Delta }:=G_o-G$, with G and $G_o$ denote the PEA’s true dynamics and the captured dynamics, respectively. Therefore, the inverse model $G_o^{-1}\left(j\omega \right)$ is exactly the trained inversion model, invLSTMs2s. Then, (21) can be rewritten as

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\left|\mathrm{\Delta }\left(j\omega \right)\right|}_2}{|G_0{\left(j\omega \right)|}_2}}\le & 1\hfill \\ \hfill {\left|\mathrm{\Delta }\left(j\omega \right)\right|}_2{\left|G_o^{-1}\left(j\omega \right)\right|}_2\le & 1\hfill \\ \hfill |G_0{\left(j\omega \right)-G\left(j\omega \right)|}_2{\left|G_o^{-1}\left(j\omega \right)\right|}_2\le & 1\hfill \\ \hfill |1-G\left(j\omega \right)G_o^{-1}{\left(j\omega \right)|}_2\le & 1.\hfill \end{array}$$

(22)

Note that $G\left(j\omega \right)G_0^{-1}\left(j\omega \right)$ represents the model of the cascaded invLSTMs2s+PEA system, i.e., when the invLSTMs2s feedforward controller is applied to the PEA alone. Thus, the tracking error subject to reference R in this case is defined as $E_{ff-only}:=R-Y=[1-G\left(j\omega \right)G_o^{-1}\left(j\omega \right)]R$, and (22) becomes

$$\underset{{\left|R\left(j\omega \right)\right|}_2\ne 0}{max}{\displaystyle \frac{|E_{ff-only}{\left(j\omega \right)|}_2}{{\left|R\left(j\omega \right)\right|}_2}}\le 1,$$

(23)

which matches the $E_{rms}$ in (24) as discussed in next subsection.

Note that ${\displaystyle \frac{|E_{ff-only}{\left(j\omega \right)|}_2}{{\left|R\left(j\omega \right)\right|}_2}}$ is the RMS tracking error (see (24)). Therefore, as long as the RMS tracking error satisfies the above (23) when only the invLSTMs2s feedforward controller is used, it is guaranteed that the proposed 2DOF controller performs better than the MPC alone.

Due to the training data limit in practice, such as data length, frequency, and amplitude limits, invLSTMs2s may not capture the PEA dynamics outside of the training frequency and/or amplitude ranges. In this case, the modeling uncertainty $\mathrm{\Delta }$ may not satisfy condition (21). To solve this issue, we introduced frequency-dependent weighting factors ${\kappa }_{ff}\left(\omega \right)$ and ${\kappa }_{fb}\left(\omega \right)$, such that the contribution of the inversion feedforward part could be flexibly adjusted based on the frequency range of the reference (see Figure 2). For example, for low-frequency (lower than the invLSTMs2s lower training frequency bounds) tracking, ${\kappa }_{ff}$ can be set small so the MPC feedback control is dominant to minimize the uncertainty effect. Also, a low-pass filter (LPF) was used in the MPC part to avoid high-frequency data that could be fed into the feedback loop.

3. Experiment Results and Discussion

The proposed invLSTMs2s+MPC 2DOF controller was implemented on a PEA stage (Nano-OP30, Mad City Labs, Madison, WI, USA) for the experimental validation. The performance was compared with that of feedforward invLSTMs2s only, an MPC only, and a traditional PI controller. All the signals were acquired using a data acquisition (DAQ) system (NI PCIe-6353, National Instruments, Austin, TX, USA) installed on a desktop workstation (Intel Xeon W-2125, RAM 32 GB), as shown in Figure 3. All the models and controllers were trained/designed using MATLAB (R2023a) and Simulink (10.7) (MathWorks, Inc., Natick, MA, USA). For all the experiments, the sampling frequency was set at 20 kHz.

3.1. invLSTMs2s Training Set Construction and Training Process

The prerequisite step for generating the invLSTMs2s training set $\{Y_{\left(ts\right)},U_{\left(ts\right)}\}$ is finding sufficient f–A pairs. To do so, 10,000 points (blue dots in Figure 4) in the f range of 25–325 Hz and amplitude range of 0–4.5 V were first randomly generated in the f–A plane, and the k-means algorithm [22] picked 1000 points (red points in Figure 4) as $\mathrm{\Omega }$’s for the training set generation. Then, (6) was used to generate $U_{\left(ts\right)}$ with $P=3$. Using $U_{\left(ts\right)}$ as the drive input of the PEA, the corresponding PEA output $Y_{\left(ts\right)}$ was collected by the DAQ device. Thus, the entire training set $(Y_{\left(ts\right)}$, $U_{\left(ts\right)})$ contained 2000 samples.

The invLSTMs2s model with 40 hidden states was used to model the PEA inverse dynamics and trained with $\lambda =32$ epochs using parallel training. A batch size of $n=20$ was selected. Therefore, the entire training process took $\lambda K_s/n=3200$ iterations. Figure 5 shows the invLSTMs2s root-mean-square error (RMSE) learning curve as the epoch number increases. It clearly shows the convergence of the invLSTMs2s training process with very low training error.

3.2. Accuracy of invLSTMs2s

The model accuracy was tested by checking how well invLSTMs2s could compensate for the PEA hysteresis, i.e., whether the hysteresis could be mostly removed when invLSTMs2s was cascaded with the PEA (connected in series). As shown in Figure 6, the original PEA hysteresis curves were measured using different sinusoidal driving voltages with the frequencies 30, 120, and 240 Hz and amplitudes of 2 and 4 V, respectively. Please note that these frequencies were intentionally selected to avoid overlapping with those in the training set to minimize the effect of data discrepancy. It can be seen that the PEA hysteresis tended to be more pronounced with increasing frequency and/or amplitude. However, the corresponding hysteresis curves of the cascaded invLSTMs2s+PEA system were barely notable for all measured frequencies and amplitudes (see Figure 6). Therefore, the invLSTMs2s inversion model was indeed accurate in accounting for the PEA’s nonlinear dynamics. This accuracy ensured that the modeling uncertainty stayed small (i.e., (21) was satisfied) and guaranteed that the output tracking performance of the proposed 2DOF invLSTMs2s+MPC controller was better than that of the MPC alone. More detailed comparisons are presented next.

3.3. Tracking Performance Comparison

The tracking performance of the proposed invLSTMs2s+MPC was demonstrated as follows. For the MPC, the prediction horizon was set to 60, and the stability condition (20) was satisfied. To design the Kalman state estimator, the covariance for the process noise and measurement noise were selected as $Q={10}^{-4}$ and $R={10}^{-6}$, respectively, based on the system’s linear model identification accuracy, the PEA stage’s measurement noise level, and the system’s noise-to-signal ratio provided by the manufacturer. Also, to avoid high-frequency noise being fed into the feedback loop, a third-order Butterworth LPF with the cut-off frequency of 600 Hz was used to avoid high-frequency noise/disturbance being fed into the feedback loop.

Although the invLSTMs2s was trained up to 300 Hz, the control bandwidth was expected to reach 350 Hz due to the generalizability of LSTM. Based on the PEA behavior observed, whose nonlinear dynamics was highly nonlinear at frequencies close or above 200 Hz, frequencies ≥200 Hz were considered as high frequencies. To comprehensively test the controller performance over broad frequency/amplitude ranges, the reference trajectories were designed to have low (≤50 Hz), mid (∼100 Hz), and high frequency (≥200 Hz) profiles: sinusoidal signals (at 30, 120, and 360 Hz), triangle signals (at 25, 100, and 300 Hz), and a multi-frequency reference $\mathrm{\Gamma }$, where $\mathrm{\Gamma }\left(t\right)=0.8sin(2\pi 5t+1.5\pi )+0.43sin\left(2\pi 50t\right)+0.12sin(2\pi 120t+1.2\pi )+0.3sin(2\pi 180t+\pi )$.

First, the accuracy of invLSTMs2s as the feedforward controller only was evaluated. Then, the proposed invLSTMs2s+MPC 2DOF controller was compared with the cases of using invLSTMs2s, the MPC, and PI feedback. The PI controller parameters were tuned to $P=0.312$ and $I=1465$ to achieve the best performance possible. The tracking errors were computed as [28]:

$$E_{rms}={\displaystyle \frac{{|r(·)-y(·)|}_2}{{\left|r(·)\right|}_2}}\times 100\%,E_{max}={\displaystyle \frac{{|r(·)-y(·)|}_{\infty }}{{\left|r(·)\right|}_{\infty }}}\times 100\%,$$

(24)

where $r(·)$ and $y(·)$ are the complex vectors of reference and output obtained through discrete Fourier transforms, respectively.

Table 1. Performance comparison.

Controller	Refs.	Sinusoidal			Triangular			$\mathbf{\Gamma }$
		30 Hz	120 Hz	360 Hz	25 Hz	100 Hz	300 Hz
invLSTMs2s	$E_{rms}$	9.35%	7.25%	3.10%	9.04%	6.00%	4.29%	5.64%
	$E_{max}$	7.91%	6.14%	2.30%	7.13%	4.74%	2.51%	3.49%
invLSTMs2s+MPC	$E_{rms}$	3.19%	2.97%	7.53%	4.05%	1.94%	6.32%	3.63%
	$E_{max}$	2.09%	2.07%	6.25%	2.85%	1.05%	4.31%	2.39%
MPC	$E_{rms}$	1.66%	6.57%	27.18%	3.28%	6.80%	21.66%	10.06%
	$E_{max}$	1.12%	4.82%	23.26%	2.63%	3.77%	12.68%	7.57%
PI	$E_{rms}$	6.37%	24.34%	45.89%	5.08%	18.73%	38.01%	8.56%
	$E_{max}$	5.50%	21.04%	39.47%	3.77%	14.88%	30.60%	5.81%

shows the tracking errors for all the controllers. Figure 7 shows the detailed tracking performance in the time domain for the selected references of (120 Hz sinusoidal signal, 100 Hz triangular signal, and $\mathrm{\Gamma }$).

invLSTMs2s+MPC vs. invLSTMs2s: First of all, condition (23) was satisfied for the case when only invLSTMs2s was used when tracking all the designed references: the RMS tracking errors were all much smaller than one. Hence, invLSTMs2s could accurately capture the PEA nonlinear dynamics and produced very small errors for tracking high-frequency references; however, relatively higher errors at the low-frequency range were observed. For example, in Table 1, the tracking error of the invLSTMs2s for both sinusoidal and triangular references at 30 and 25 Hz were a few percent higher than those for 360 and 300 Hz, respectively. This was mainly because these frequencies were very close to the lower bound of the training frequency range. Thus, limited knowledge of the system dynamics at these frequencies was available when training invLSTMs2s. This shortage was compensated for by the MPC, which was designed using a linear model that primarily captured the PEA dynamics at low frequencies. The invLSTMs2s+MPC outperformed invLSTMs2s at the low-frequency range. The tracking errors of the invLSTMs2s+MPC for both sinusoidal and triangular signals at 30 and 25 Hz were reduced by more than 50% compared to those of invLSTms2s. As shown in Figure 7, the low-frequency tracking error was greatly reduced by invLSTMs2s+MPC compared to invLSTMs2s. As the PEA nonlinear dynamics became more pronounced at higher frequencies, the advantage of invLSTMs2s became more significant: invLSTMs2s achieved the highest accuracy for 360 Hz sinusoidal and 300 Hz triangular trajectories. Meanwhile, as the linear MPC was not designed for the PEA’s high-frequency control, it negatively impacted the performance of invLSTMs2s. Therefore, bigger errors were observed (see Table 1). However, compared to the errors of MPC, the performance was still greatly improved, and the achieved accuracy was satisfactory.

invLSTMs2s+MPC vs. MPC: The proposed invLSTMs2s+MPC was compared with the use of the MPC alone to clearly demonstrate the necessity of integrating the invLSTMs2s feedforward control. As the performance improvement condition (23) was satisfied for all the tracking cases ($E_{rms}<<1$ for the use of invLSTMs2s alone in Table 1), it was expected that the proposed 2DOF invLSTMs2s+MPC controller would outperform the MPC. Table 1 shows that the invLSTMs2s+MPC had less tracking error than the MPC, especially for tracking mid-, high-frequency, and the multi-frequency $\mathrm{\Gamma }$ trajectories. However, it made sense that the tracking errors of invLSTMs2s+MPC at 30 and 25 Hz were a little bit higher than those of the MPC, since the MPC is more accurate than invLSTMs2s at low frequencies, as aforementioned. Note that the slight error increase can be eliminated by setting the feedforward weighting coefficient ${\kappa }_{ff}\left(\omega \right)$ to a very small value to make sure the 2DOF controller takes full advantage of the MPC at low frequencies.

invLSTMs2s+MPC vs. PI: To demonstrate the significance of this work, invLSTMs2s+MPC was compared to the case of using a PI controller alone. In Table 1, it can be observed that the tracking error of the PI controller was more pronounced with increasing frequency. For the mid and high frequencies, the PI presented the highest tracking error among all controllers. This is mainly because the PI feedback controller was driven by the errors between references and PEA outputs, and the control bandwidth was low due to the sluggish response of this controller subject to fast error changes.

In summary, invLSTMs2s+MPC is capable of achieving high tracking accuracy for various reference trajectories over a broad bandwidth. The superiority of this 2DOF controller comes from taking the advantages of feedforward invLSTMs2s for the PEA’s high-frequency nonlinear dynamics compensation and MPC for low-frequency position tracking precision.

3.4. Application Demonstration of the invLSTMs2s+MPC

As an application demonstration, invLSTMs2s+MPC was further adopted to control a two-dimensional (X–Y lateral plane motion) PEA nanopositioning platform to track the pattern of character string “ISU” along the lateral X–Y plane. The tracked patterns achieved by different controllers are compared in Figure 8. It was designed to complete the motion within 0.11 s, with the stage average traveling speed of 320 $\mathrm{\mu }$m/sec. The designed X and Y direction trajectories and tracking results for different controllers are shown in Figure 9a,b, respectively. According to the results, it is clear that the proposed invLSTMs2s+MPC performed better and could achieve very high accuracy. This clearly demonstrated the capability of invLSTMs2s+MPC in two-dimensional positioning tracking control in nanofabrication applications, such as nano/micromanufacturing systems in which the PEA stages are required to track the desired fabrication pattern [29,30].

As invLSTMs2s does not assume any specific forms of nonlinearity, the outstanding performance of invLSTMs2s+MPC directly implies that it can be employed in various nonlinear systems for applications where position tracking is essential, such as precise control of friction force generated by deformable actuators in human–robot reaction [31], AFM applications for material mechanical characterization [32], as well as adaptive optics for high-resolution imaging [33], just to name a few.

4. Conclusions

In this work, we proposed a 2DOF controller, invLSTMs2s+MPC, for nanopositioning control of PEA devices with a broad bandwidth, high precision, and fast response. The key novelty of the proposed controller includes (1) the idea of implementing an LSTMs2s structure to model the PEA inverse dynamics for system nonlinearity compensation and (2) integrating invLSTMs2s and a linear MPC to improve the overall tracking performance. As the invLSTMs2s model can be trained offline and a linear MPC is used during real-time control, the proposed controller can achieve broad-bandwidth output tracking with its low computation requirement.

In the future, the work will focus on designing a dynamics separation mechanism to efficiently optimize the collaboration between the feedforward and feedback components. Also, we will explore other modeling network structures to improve the training speed further and increase the modeling bandwidth, robustness, and accuracy.