Iterative Tuning of Feedforward Controller with Precise Time-Delay Compensation for Precision Motion System

Abstract

The accuracy of feedforward control model including system time-delay significantly affects the position tracking performance in a precision motion system. In this paper, an iterative tuning method for feedforward control with precise time-delay compensation is proposed. First, considering system time-delay from actuator, sensor, calculation, and communication in real platform, a feedforward control model with time-delay compensation is established, and a nonlinear objective function with time-delay is designed based on the measured data of a finite time task, to minimize the position tracking error. Second, in order to deal with both the nonlinear objective function and also unknown disturbances and noise in the real system, an optimization strategy combining the Gauss–Newton iterative (GNI) scheme and instrumental variable (IV) is proposed to realize the unbiased estimation of the feedforward parameters and precise delay time. Finally, with the identified feedforward control parameters, the precise system time-delay which is a nonintegral multiple of the sampling period is compensated accurately for the feedforward control with accurate path planning time-shift in the implementation. The effectiveness of the proposed feedforward parameter tuning and precise time-delay compensation scheme is verified by the simulation and also experimental result on a precision motion platform with obvious position tracking performance improvement.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In real control systems, the sensors, actuators, data communication, and signal processing can all generate some time-delay [1, 2]. Without compensation, the time-delay will lead to reduced system control bandwidth, will lead to slow response and even affects the system stability. On the one hand, the characteristics of time-delay can also be utilized in the active vibration absorber to obtain better vibration suppression performance [3, 4]. In precision motion systems, the feedforward control is introduced to compensate the time-delay by injecting the control signal in advance [5], which can significantly improve the tracking accuracy. This feedforward method is widely used in the high-speed and high-precision motion control systems such as photolithography equipment [6, 7], machine tool [8], and atomic force microscope system [9].

The parameter accuracy of the feedforward control significantly affects the control performance of the precision motion system [10]. In the feedback-feedforward control structure, when the feedforward model is equal to the inverse of the controlled plant, the position error can be effectively compensated. The model-based feedforward control is commonly used in the motion system, but it needs to determine the precise model of the plant in advance. When the system performs a finite time task, such as point-to-point motion, the measured data contains the knowledge of the plant. Based on this, the iterative feedforward tuning (IFT) [11, 12] method utilizes the measured data to optimize the feedforward parameters without the need for detailed knowledge of the plant. In the IFT method, the least square (LS) method based on instrumental variable (IV) is adopted to eliminate the effect of noise unrelated to the input signal, and the unbiased estimation of feedforward parameters is obtained. It can be seen that the IFT establishes a connection between the feedforward tuning and closed-loop system identification and clarifies the direction of feedforward parameter regulation. In terms of parameter estimation variance, an iterative refined instrumental variable is constructed to achieve optimal accuracy [13, 14], and the Kalman filtering (KF) approach is introduced into the IV-IFT framework, which enables unbiased parameter estimation with zero asymptotic variance [15]. Then, the IFT is extended to flexible motion systems. For the non-minimum phase system, the stability problem of model inversion is solved by the input shaping method [16, 17]. The unbiased parameter estimation with optimal accuracy in terms of variance is obtained for feedforward controllers with a rational basis [18], and the feedforward control with rational basis functions can enable higher performance and more enhanced extrapolation capabilities than polynomial basis functions [18, 19]. A high-order IFT algorithm is proposed by introducing the iterative domain into the IFT, and IV is also employed to tolerate the noise data [6].

The inherent time-delay in precision motion systems causes nonlinearity issue. In existing IFT methods, the linear feedforward model is adopted, which does not solve the coupling problem between the delay time and model parameters. Besides, the LS cannot be directly used in the nonlinear system with time-delay [20], so the tuning accuracy of feedforward parameters cannot be guaranteed. In addition to LS, the gradient-descent method (GD) is used to tune a feedforward controller with force ripple compensation [21]. The Newton iterative (NI) method is used to tune the feedforward controller iteratively for non-minimum phase systems [16] and design a feedforward controller for multi-input multioutput systems [7]. These two methods are suitable for the feedforward tuning of nonlinear systems, but they cannot deal with disturbances in real systems.

Many methods are proposed to identify the delay time. The Bode diagram of the plant can be used to fit the delay time [22, 23]. The time-delay term is linearly parameterized by Taylor expansion, and a new adaptive law is constructed to identify the delay time [20]. Then, an approximate nonlinear LS is proposed to simultaneously estimate the delay time and dynamic parameters in the continuous system [24]. Moreover, the Pade approximant is applied to replace the time-delay to realize parameter estimation [25]. In addition, the Newton iterative and separable methods are applied to identify the delay time of nonintegral multiple period for the discrete system [26]. Therefore, the approximate linear transformation is a key step to realize the identification of the time-delay.

As we can see, the existing IFT based on IV is not suitable for nonlinear systems with time-delay, and the feedforward adjustment method based on GD and NI cannot deal with disturbances in real systems. To solve these problems, this paper presents a precise parameter tuning method of feedforward control with time-delay compensation. The main contributions of this paper are as follows: (1) fully considering the comprehensive system time-delay from actuator, sensor, calculation, and communication in real platform, a nonlinear objective function is proposed based on the measured data of a finite time task for iterative tuning of the feedforward parameters, to minimize the position tracking error; (2) in order to handle the proposed nonlinear objective function and also tolerate unknown disturbances and noise in real system, a desired optimization strategy combining the Gauss–Newton iterative (GNI) scheme and instrumental variable (IV) is proposed in this paper to realize the unbiased estimation of the feedforward parameters and precise delay time, which is the key innovation of this paper; (3) the identified precise system time-delay which is a nonintegral multiple of the sampling period, is exactly compensated in the feedforward control with accurate path planning time-shift. The simulation illustration and experimental validation demonstrate the advantages of the proposed control strategy.

The paper is organized as follows: in Section 2, the mathematical model of the feedforward controller with time-delay compensation is established, and a feedforward parameter tuning method considering the time-delay is elaborated; in Section 3, a discrete realization method of feedforward control with time-delay compensation is proposed; in Sections 4 and 5, the effectiveness of the proposed method is verified in the simulation example and experiment on an air floating precision motion platform.

2. Iterative Tuning Method of Feedforward Parameters with Time-Delay

In order to achieve the position tracking control, a feedback-feedforward control system is established in this paper, as shown in Figure 1. The error signal between the reference position and the feedback position is processed by the feedback controller to generate the control signal, which is loaded on the motor to generate thrust and the controlled plant moves. The measured position is fed back to the loop to build the closed-loop control, which realizes the position tracking control on the premise of ensuring the stability of the system, and has the ability of restraining the disturbance. On this basis, the feedforward control improves the position tracking accuracy by adding an input signal in the forward channel.

[figure omitted; refer to PDF]

In Figure 1, the unknown controlled plant $Ps$ is a single-input single-output (SISO) and continuous time-invariant system with time-delay, and $Ps$ is considered with unit numerator for motion systems with dominant rigid-body dynamics [6, 11, 13]. $C_{fb}s$ is the feedback controller, and $C_{ff}s$ is the feedforward controller. $\mathbf{r}$ denotes the reference position and is a known third-order multisegment polynomial trajectory, ${\mathbf{u}}_{fb}$ is the feedback control signal, and ${\mathbf{u}}_{ff}$ is the feedforward control signal. The unknown disturbance $\mathbf{w}$ is assumed to be given by $\mathbf{w}=Hs\mathbf{\epsilon }$, where $\mathbf{\epsilon }$ is the normally distributed white noise with zero mean and variance ${\lambda }_{\epsilon }^2$, $Hs$ is the shaping filter that makes $\mathbf{w}$ a colored noise; and hence, $\mathbf{w}$ and $\mathbf{r}$ are uncorrelated. In a finite time task [11], the planned path $\mathbf{r}$ is input to the system and the system implements point-to-point motion, ${\mathbf{e}}_m$ is the measured position error signal, and ${\mathbf{y}}_m$ is the position signal.

From Figure 1, the position error is ${\mathbf{e}}_m=G1-{PC}_{ff}\mathbf{r}–GP\mathbf{w}$, where $G=1/1+{PC}_{fb}$, where the Laplace transform symbol $s$ is omitted. It shows that when the feedforward controller $C_{ff}$ is equal to the inversion of the plant $P^{-1}$, $G1-{PC}_{ff}\mathbf{r}$ in the error ${\mathbf{e}}_m$ is zero and this is the goal of the feedforward parameters tuning.

2.1. Feedforward Control Model

In order to realize the feedforward control and tune the feedforward parameters, the mathematical model of the feedforward controller $C_{ff}$ is established. Considering the time-delay in the system, the time-advance term $e^{\tau s}$ is introduced in the linear parametric model [11–17] and $C_{ff}$ is parameterized as follows:$$\begin{array}{}\text{(1)}& C_{ff}={\theta }_n{s}^n+{\theta }_{n-1}{s}^{n-1}S+\cdots +{\theta }_0{e}^{\tau s}={\displaystyle \sum _{i=0}^n{\theta }_i{s}^i}{e}^{\tau s},\end{array}$$where $n$ can be determined by means of a model order selection procedure [27].${\theta }_i$ are the feedforward parameters, and $\tau$ is the delay time. The feedforward parameters vector is $\mathbf{\theta }={{\theta }_0,{\theta }_1,\dots ,{\theta }_n,\tau }^T$.

2.2. Feedforward Parameter Tuning Method

In this method, the feedforward parameters are tuned iteratively by using the measurement data of the finite time task of the system to achieve the feedforward control goal. In a task, the system starts from the static state and executes point-to-point motion to obtain the complete motion data. The measured signal vector $\mathbf{u}$ in each task is defined as $\mathbf{u}={u1,u2,\dots ,uN}^T\in {\mathbf{R}}^N$, where $ut$ is a measurement at time instant $t$ for $t=1,\dots ,N{T}_s$, where $T_s$ is the sampling time and N is the number of sampling.

In the $k$th task, the feedforward controller is $C_{ff}^k$; from equation (1), we obtain$$\begin{array}{}\text{(2)}& C_{ff}^k={\displaystyle \sum _{i=0}^n{\theta }_i^k{s}^i}{e}^{{\tau }^{k}s}.\end{array}$$

And ${\mathbf{\theta }}^k={{\theta }_0^k,{\theta }_1^k,\dots ,{\theta }_n^k,{\tau }^k}^T$ is known but may not be equal to the plant parameters. The task is executed, and the measured signals ${\mathbf{e}}_m^k$ and ${\mathbf{y}}_m^k$ are given by$$\begin{array}{}\text{(3)}& {\mathbf{e}}_m^k={\mathbf{e}}_r^k+{\mathbf{e}}_w=G1-{PC}_{ff}^k\mathbf{r}–GP\mathbf{w},\text{(4)}& {\mathbf{y}}_m^k={\mathbf{y}}_r^k+{\mathbf{y}}_w=GP{C}_{fb}+C_{ff}^k\mathbf{r}+GP\mathbf{w}.\end{array}$$

In ${\mathbf{e}}_m^k$, the error caused by the reference $\mathbf{r}$ is ${\mathbf{e}}_r^k=G1-{PC}_{\text{ff}}^k\mathbf{r}$ and the error caused by the disturbance is ${\mathbf{e}}_w=-GP\mathbf{w}$. In ${\mathbf{y}}_m^k$, the position caused by the reference $\mathbf{r}$ is ${\mathbf{y}}_r^k=GP{C}_{fb}+C_{ff}^k\mathbf{r}$ and the position caused by the disturbance is ${\mathbf{y}}_w=GP\text{w}$.

Then, as shown in Figure 2, the feedforward parameters are tuned, and the iterative format is$$\begin{array}{c}\text{(5)}& {\mathbf{\theta }}^{k+1}={\mathbf{\theta }}^k+{\mathbf{\theta }}^{\Delta },C_{ff}^{k+1}={\displaystyle \sum _{i=0}^n{\theta }_i^{k+1}{s}^i}{e}^{{\tau }^{k+1}s}={\displaystyle \sum _{i=0}^n{\theta }_i^k+{\theta }_i^{\Delta }{s}^i}{e}^{{\tau }^k+{\tau }^{\Delta }s},\end{array}$$where ${\mathbf{\theta }}^{\mathrm{\Delta }}={{\theta }_0^{\mathrm{\Delta }},\dots ,{\theta }_n^{\mathrm{\Delta }},{\tau }^{\mathrm{\Delta }}}^T$ is the feedforward correction vector, ${\theta }_i^{\mathrm{\Delta }}$ are the feedforward parameters corrections, and ${\tau }^{\mathrm{\Delta }}$ is the delay time correction, which constitute the feedforward correction model $C_{ff}^{\mathrm{\Delta }}={\displaystyle {\sum }_{i=0}^n{\theta }_i^{\mathrm{\Delta }}{s}^i}{e}^{{\tau }^{\mathrm{\Delta }}s}$. In a linear system, $C_{ff}^{k+1}$ is a linear superposition of $C_{ff}^k$ and $C_{ff}^{\mathrm{\Delta }}$ [11, 12]. Here, $C_{ff}$ is nonlinearly parameterized and the feedforward parameters and delay time need to be superimposed separately.

[figure omitted; refer to PDF]

Next, an objective function is established by the measurement data ${\mathbf{e}}_m^k$ and ${\mathbf{y}}_m^k$ to calculate the feedforward parameter correction vector ${\theta }^{\mathrm{\Delta }}$.

2.2.1. Objective Function for Feedforward Parameter Identification

The objective function $J{\mathbf{\theta }}^{\mathrm{\Delta }}$ is established based on the measured data, and the correction vector ${\mathbf{\theta }}^{\mathrm{\Delta }}$ can be obtained by minimizing $J{\mathbf{\theta }}^{\mathrm{\Delta }}$. Then, $C_{ff}^k$ is updated to $C_{ff}^{k+1}$ to minimize the position error.

First, the residual vector is established based on the measured data. In order to associate the measured data with ${\mathbf{\theta }}^{\mathrm{\Delta }}$, $C_{ff}^{k+1}-C_{ff}^k$ is introduced into equation (4), which contains ${\mathbf{\theta }}^{\mathrm{\Delta }}$. Assuming the disturbance $\mathbf{w}=0$, equations (3) and (4) can be transformed into$$\begin{array}{}\text{(6)}& {\mathbf{e}}_m^k=\mathbf{r}G1-{\text{PC}}_{ff}^k,\text{(7)}& {\mathbf{y}}_m^k{C_{fb}+C_{ff}^k}^{-1}{C}_{ff}^{k+1}-C_{ff}^k=\mathbf{r}GP{C}_{ff}^{k+1}-C_{ff}^k.\end{array}$$

The following formula can be obtained by subtracting equation (7) from equation (6):$$\begin{array}{}\text{(8)}& {\mathbf{e}}_m^k-{\mathbf{y}}_m^k{C_{fb}+C_{ff}^k}^{-1}{C}_{ff}^{k+1}-C_{ff}^k=\mathbf{r}GP{P}^{-1}-C_{ff}^{k+1}.\end{array}$$

Then, $C_{fb}+C_{ff}^{k}1/s^{n+1}{e}^{-{\tau }^k+{\tau }^{\mathrm{\Delta }}s}$ is introduced to filter the left polynomial of equation (8), which can avoid the calculation of the transfer function inversion, the differential term, and the advance term, and the residual vector is defined as$$\begin{array}{c}\text{(9)}& \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }={\mathbf{e}}_m^k-{\mathbf{y}}_m^k{C_{fb}+C_{ff}^k}^{-1}{C}_{ff}^{k+1}-C_{ff}^k{C}_{fb}+C_{ff}^k\frac{1}{s^{n+1}}{e}^{-{\tau }^k+{\tau }^{\Delta }s}={\mathbf{e}}_m^{k}A\frac{1}{s}{e}^{-{\tau }^{\Delta }s}-{\mathbf{y}}_m^k{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k+{\theta }_i^{\Delta }{s}^i}{s^{n+1}}-{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}\frac{1}{s}{e}^{-{\tau }^{\Delta }s},\end{array}$$where $A=C_{fb}{e}^{-{\tau }^{k}s}/s^n+{\displaystyle {\sum }_{i=0}^n{\theta }_i^k{s}^i/s^n}$ and the time-delay $e^{-{\tau }^{k}s}$ and $e^{-{\tau }^{\mathrm{\Delta }}s}$ can be approximated by first-order Pade polynomial [25].

Based on the residual vector, $J{\theta }^{\mathrm{\Delta }}$ is established as follows:$$\begin{array}{}\text{(10)}& J{\mathbf{\theta }}^{\Delta }={\mathbf{\epsilon }}^T{\mathbf{\theta }}^{\Delta }\mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }={\displaystyle \sum _{j=1}^N{\mathbf{\epsilon }}^{2}j}.\end{array}$$

The optimization problem is ${\mathbf{\theta }}^{\mathrm{\Delta },\text{opt}}={\min }_{{\mathbf{\theta }}^{\mathrm{\Delta }}}J{\mathbf{\theta }}^{\mathrm{\Delta }}$. When $J{\mathbf{\theta }}^{\mathrm{\Delta }}$ is minimized, the residual vector $\mathbf{\epsilon }{\mathbf{\theta }}^{\mathbf{\Delta }}$ will be zero without considering the disturbance and modeling error. Since $\mathbf{\epsilon }{\mathbf{\theta }}^{\mathrm{\Delta }}$ is derived from equation (8), $C_{ff}^{k+1}=P^{-1}$ can be obtained from the right side of equation (8) equal to 0; that is, the feedforward parameters after tuning can be equal to the true plant parameters. Then, in the next task, ${\mathbf{e}}_r^{k+1}=G1-{PC}_{ff}^{k+1}\mathbf{r}=0$, the position error is minimized, and the feedforward control goal is achieved.

2.2.2. Iterative Identification of Parameters

After the definition of the objective function, the Gauss–Newton iteration method based on instrumental variable (IV-GNI) is proposed to realize the unbiased estimation of the feedforward correction vector ${\mathbf{\theta }}^{\mathrm{\Delta }}$.

According to GNI, $J{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ is constructed by introducing the parameter variation $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ into $J{\mathbf{\theta }}^{\mathrm{\Delta }}$. Through the first-order Taylor expansion, $J{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ becomes a linear function of $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$:$$\begin{array}{c}\text{(11)}& J{\mathbf{\theta }}^{\Delta }+\Delta {\mathbf{\theta }}^{\Delta }={\mathbf{\epsilon }}^T{\mathbf{\theta }}^{\Delta }+\Delta {\mathbf{\theta }}^{\Delta }\mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }+\Delta {\mathbf{\theta }}^{\Delta }={\mathbf{\epsilon }}^T{\mathbf{\theta }}^{\Delta }\mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }+2{\mathbf{\epsilon }}^T{\mathbf{\theta }}^{\Delta }\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}\Delta {\mathbf{\theta }}^{\Delta }+{\Delta {\mathbf{\theta }}^{\Delta }}^T{\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}}^T\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}\Delta {\mathbf{\theta }}^{\Delta },\end{array}$$where $\partial \mathbf{\epsilon }{{\mathbf{\theta }}^{\mathrm{\Delta }}/\partial {\mathbf{\theta }}^{\mathrm{\Delta }}}^T$ is the derivative of $\mathbf{\epsilon }{\mathbf{\theta }}^{\mathrm{\Delta }}$ with respect to ${\mathbf{\theta }}^{\mathrm{\Delta }}$:$$\begin{array}{}\text{(12)}& \frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}=-{\mathbf{y}}_m^k\frac{1}{s^{n+1}},\dots ,{\mathbf{y}}_m^k\frac{s^n}{s^n},{\mathbf{e}}_m^{k}A+{\mathbf{y}}_m^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}{e}^{-{\tau }^{\Delta }s}.\end{array}$$

To minimize $J{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathbf{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$, $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ can be calculated by making the derivative of $J{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ zero:$$\begin{array}{}\text{(13)}& \Delta {\mathbf{\theta }}^{\Delta }=-{{\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}}^T\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}}^{-1}{\frac{\partial \mathbf{\epsilon }{\mathbf{\theta }}^{\Delta }}{\partial {{\mathbf{\theta }}^{\Delta }}^T}}^T\mathbf{\epsilon }{\mathbf{\theta }}^{\Delta },\end{array}$$which is in the form of LS [27, 28]. Then, ${\mathbf{\theta }}^{\mathrm{\Delta }}$ can be identified iteratively based on $\mathbf{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$. Let ${\mathbf{\theta }}^{\mathrm{\Delta },l}$ be the $l$th estimate of ${\mathbf{\theta }}^{\mathbf{\Delta }}$, and the iterative form is given by$$\begin{array}{}\text{(14)}& {\mathbf{\theta }}^{\Delta ,l+1}={\mathbf{\theta }}^{\Delta ,l}+\Delta {\mathbf{\theta }}^{\Delta ,l},\end{array}$$where $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta },l}$ can be obtained by equation (13) based on the known ${\mathbf{\theta }}^{\mathrm{\Delta },l}$:$$\begin{array}{}\text{(15)}& \Delta {\mathbf{\theta }}^{\Delta ,l}={{{\mathbf{X}}_m^l}^T{\mathbf{X}}_m^l}^{-1}{{\mathbf{X}}_m^l}^T{\mathbf{Y}}_m^l.\end{array}$$In the above formula,$$\begin{array}{c}\text{(16)}& {\mathbf{X}}_m^l={\mathbf{y}}_m^k\frac{1}{s^{n+1}},\dots ,{\mathbf{y}}_m^k\frac{s^n}{s^{n+1}},{\mathbf{e}}_m^{k}A+{\mathbf{y}}_m^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}{e}^{-{\tau }^{\Delta ,l}s}\in {\mathbf{R}}^{N\times n+2},{\mathbf{Y}}_m^l={\mathbf{e}}_m^{k}A+{\mathbf{y}}_m^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}\frac{1}{s}{e}^{-{\tau }^{\Delta ,l}s}-{\mathbf{y}}_m^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k+{\theta }_i^{\Delta ,l}{s}^i}{s^{n+1}}\in {\mathbf{R}}^{N\times 1}.\end{array}$$

In this method, the optimization of $J{\mathbf{\theta }}^{\mathrm{\Delta }}$ is transformed into the optimization of $J{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$, which is an ingenious transformation.

From equations (3) and (4), ${\mathbf{e}}_m^k$ and ${\mathbf{y}}_m^k$ are constructed based on $\mathbf{r}$, which provides a persistence of excitation condition and ensures that ${{\mathbf{X}}_m^l}^T{\mathbf{X}}_m^l$ is nonsingular, so $\Delta {\mathbf{\theta }}^{\mathrm{\Delta },l}$ exists. As the motion starts from the static state and the initial state is 0, it will not affect the identification of parameters. In the iteration, the initial value is selected as ${\mathbf{\theta }}^{\mathrm{\Delta }}=0$, which is located in the convex interval of the objective function containing the required extreme points, so that the algorithm can converge to the correct value. See Appendix for the proof process.

The influence of the disturbance in the measured data on $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta },l}$ is explained below. In Figure 1, $\mathbf{w}$ contains the random disturbance in the system and the unmodeled factors when establishing the feedforward model. From equation (16), we can get ${\mathbf{X}}_m^l={\mathbf{X}}_r^l+{\mathbf{X}}_w^l$ and ${\mathbf{Y}}_m^l={\mathbf{Y}}_r^l+{\mathbf{Y}}_w^l$, where$$\begin{array}{c}\text{(17)}& {\mathbf{X}}_r^l={\mathbf{y}}_r^k\frac{1}{s^{n+1}},\dots ,{\mathbf{y}}_r^k\frac{s^n}{s^{n+1}},{\mathbf{e}}_r^{k}A+{\mathbf{y}}_r^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}{e}^{-{\tau }^{\Delta ,l}s},{\mathbf{Y}}_r^l={\mathbf{e}}_r^{k}A+{\mathbf{y}}_r^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}\frac{1}{s}{e}^{-{\tau }^{\Delta ,l}s}-{\mathbf{y}}_r^k\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k+{\theta }_i^{\Delta ,l}{s}^i}{s^{n+1}},\\ {\mathbf{X}}_w^l=\mathbf{w}GP\frac{1}{s^{n+1}},\dots ,\mathbf{w}GP\frac{s^n}{s^{n+1}},\mathbf{w}{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}-AGP{e}^{-{\tau }^{\Delta ,l}s},\\ {\mathbf{Y}}_w^l=\mathbf{w}{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k{s}^i}{s^n}-AGP\frac{1}{s}{e}^{-{\tau }^{\Delta ,l}s}-\mathbf{w}\cdot {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k+{\theta }_i^{\Delta ,l}{s}^i}{s^{n+1}}GP.\end{array}$$

According to equation (15), based on the measured data with disturbance, we can get $\mathrm{\Delta }{\mathbf{\theta }}_m^{\mathrm{\Delta },l}=\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}+\mathbf{A}$, and the matrix inversion formula is used in the derivation [27, 28], where$$\begin{array}{c}\text{(18)}& \Delta {\mathbf{\theta }}_r^{\Delta ,l}={{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}{{\mathbf{X}}_r^l}^T{\mathbf{Y}}_r^l,\mathbf{A}=-{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}\mathbf{B}{\mathbf{I}+{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}\mathbf{B}}^{-1}{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}{{\mathbf{X}}_r^l}^T{\mathbf{Y}}_r^l\\ +\mathbf{I}-{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}\mathbf{B}{\mathbf{I}+{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}\mathbf{B}}^{-1}{{{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l}^{-1}\mathbf{C},\\ \mathbf{B}={{\mathbf{X}}_r^l}^T{\mathbf{X}}_w^l+{{\mathbf{X}}_w^l}^T{\mathbf{X}}_r+{{\mathbf{X}}_w^l}^T{\mathbf{X}}_w^l,\\ \mathbf{C}={{\mathbf{X}}_r^l}^T{\mathbf{Y}}_w^l+{{\mathbf{X}}_w^l}^T{\mathbf{Y}}_r^l+{{\mathbf{X}}_w^l}^T{\mathbf{Y}}_w^l,\end{array}$$in which $\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}$ is determined by the system parameters and is not affected by the disturbance, which is accurate. However, due to the presence of $\mathbf{w}$, $\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}$ cannot be obtained directly. The reasons are as follows: the mathematical expectation of $\mathrm{\Delta }{\mathbf{\theta }}_m^{\mathrm{\Delta },l}$ is $E\mathrm{\Delta }{\mathbf{\theta }}_m^{\mathrm{\Delta },l}=\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}+E\mathbf{A}$. It can be seen from equations (17) and (18) that $\mathbf{A}$ contains ${{\mathbf{X}}_r^l}^T{\mathbf{X}}_r^l$ and ${{\mathbf{X}}_r^l}^T{\mathbf{Y}}_r^l$, where the reference $\mathbf{r}$ is auto-correlated, and $\mathbf{A}$ contains ${{\mathbf{X}}_w^l}^T{\mathbf{X}}_w^l$ and ${{\mathbf{X}}_w^l}^T{\mathbf{Y}}_w^l$, where $\mathbf{w}$ is autocorrelated, so that there exists $E\mathbf{A}\ne 0$. Therefore, when the disturbance exists, $\mathrm{\Delta }{\mathbf{\theta }}_m^{\mathrm{\Delta },l}$ obtained by the Gauss–Newton iterative method cannot be guaranteed to be the unbiased estimate of $\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}$.

Since equation (15) has the least square calculation format, the instrumental variable method in linear system identification can be employed to eliminate the deviation caused by the disturbance $\mathbf{w}$. The instrumental variable $\mathbf{Z}$ is introduced, which is the same dimension as ${\mathbf{X}}_m^l$ and satisfies the following limit characteristics: ${\lim }_{N⟶\infty }1/N{\mathbf{Z}}^T\mathbf{w}=0$ and ${\lim }_{N⟶\infty }1/N{\mathbf{Z}}^T{\mathbf{X}}_r^l=\mathbf{Q}$, in which $\mathbf{Q}$ is nonsingular. That is, $\mathbf{Z}$ needs to be independent of $\mathbf{w}$ and is strongly correlated to ${\mathbf{e}}_r^k$ and ${\mathbf{y}}_r^k$ in ${\mathbf{X}}_r^l$. Obviously, the ideal choice is $\mathbf{Z}={\mathbf{X}}_r^l$, but ${\mathbf{X}}_r^l$ cannot be obtained when $\mathbf{w}\ne 0$. Since both ${\mathbf{e}}_r^k$ and ${\mathbf{y}}_r^k$ are correlated to $\mathbf{r}$, one choice of $\mathbf{Z}$ is $\mathbf{Z}=\mathbf{r},\dots ,\mathbf{r}{s}^n,\mathbf{r}{s}^{n+1}$, which can meet the limit characteristics.

From equation (15), we can get ${\mathbf{Y}}_m^l={\mathbf{X}}_m^l\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta },l}$, and ${\mathbf{Z}}^T$ is introduced to obtain ${\mathbf{Z}}^T{\mathbf{Y}}_m^l={\mathbf{Z}}^T{\mathbf{X}}_m^l\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}$. Then, the IV estimate of $\mathbf{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta },l}$ is $\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}={{\mathbf{Z}}^T{\mathbf{X}}_m^l}^{-1}{\mathbf{Z}}^T{\mathbf{Y}}_m^l$. According to the limit characteristics of $\mathbf{Z}$, we can get$$\begin{array}{c}\text{(19)}& \underset{N⟶\infty }{\lim }\Delta {\mathbf{\theta }}_{\text{IV}}^{\Delta ,l}={\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}}^{-1}\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}\cdot \Delta {\mathbf{\theta }}_r^{\Delta ,l}-\underset{N⟶\infty }{\lim }{\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}}^{-1}\frac{{\mathbf{Z}}^T{\mathbf{X}}_w^l}{N}{\mathbf{I}+{\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}}^{-1}\frac{{\mathbf{Z}}^T{\mathbf{X}}_w^l}{N}}^{-1}{\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}}^{-1}\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}\Delta {\mathbf{\theta }}_r^{\Delta ,l}\\ +\underset{N⟶\infty }{\lim }{\frac{{\mathbf{Z}}^T{\mathbf{X}}_r^l}{N}+\frac{{\mathbf{Z}}^T{\mathbf{X}}_w^l}{N}}^{-1}\frac{{\mathbf{Z}}^T{\mathbf{Y}}_w^l}{N}\\ =\Delta {\mathbf{\theta }}_r^{\Delta ,l}.\end{array}$$

That is, the instrumental variable $\mathbf{Z}$ eliminates the deviation caused by the disturbance $\mathbf{w}$ and $\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}$ is the consistent estimate of $\mathrm{\Delta }{\mathbf{\theta }}_r^{\mathrm{\Delta },l}$. When the number of measured data $N$ is large enough, $\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}$ has an unbiased estimation characteristic. Therefore, ${\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta }}$ obtained by the iteration in equation (14) based on $\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}$ is also the unbiased estimation of the feedforward parameter correction vector ${\mathbf{\theta }}^{\mathrm{\Delta }}$, which can minimize $J{\mathbf{\theta }}^{\mathrm{\Delta }}$. Since ${\mathbf{\theta }}^{\mathrm{\Delta }}$ is calculated by a complete finite time task data, the whole tuning procedure is conducted offline.

Based on the above optimization for $J{\mathbf{\theta }}^{\mathrm{\Delta }}$, the iterative tuning of feedforward controller can be executed. The flow chart of iterative tuning of feedforward parameters is shown in Figure 3, and the procedure is as follows:

(1) $k=0$, the current feedforward controller is $C_{ff}^k=C_{ff}^0$, measures ${\mathbf{e}}_m^k$ and ${\mathbf{y}}_m^k$ with $C_{ff}^k$ applied to the system.

(2) Iterative calculation of ${\mathbf{\theta }}^{\mathrm{\Delta }}$. As shown in Figure 3, when $l=0$, set the initial correction vector ${\mathbf{\theta }}^{\mathrm{\Delta },l}={\mathbf{\theta }}^{\mathrm{\Delta },0}=\mathbf{0}$. Construct ${\mathbf{X}}_m^l$ and ${\mathbf{Y}}_m^l$ from equation (16) and introduce the instrumental variable $\mathbf{Z}=\mathbf{r},\dots ,\mathbf{r}{s}^n,\mathbf{r}{s}^{n+1}$; solve $\mathrm{\Delta }{\mathbf{\theta }}_{\text{IV}}^{\mathrm{\Delta },l}={{\mathbf{Z}}^T{\mathbf{X}}_m^l}^{-1}{\mathbf{Z}}^T{\mathbf{Y}}_m^l$, and then, in the next iteration, ${\mathbf{\theta }}^{\mathrm{\Delta },l+1}={\mathbf{\theta }}^{\mathrm{\Delta },l}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta },l}$. If each element in ${\mathbf{\theta }}^{\mathrm{\Delta }}$ converges to its stable value, the iterative calculation of ${\mathbf{\theta }}^{\mathbf{\Delta }}$ ends and proceeds to the next step. Otherwise, $l=l+1$, and repeat this step.

(3) Construct the new feedforward controller $C_{ff}^{k+1}={\displaystyle {\sum }_{i=0}^n{\theta }_i^k+{\theta }_i^{\mathrm{\Delta }}{s}^i}{e}^{{\tau }^k+{\tau }^{\mathrm{\Delta }}s}$, and this adjustment is completed. Then, set $k=k+1$ and proceed from step 1.

[figure omitted; refer to PDF]

In this section, the iterative tuning method of the feedforward parameters in the motion system with time-delay is elaborated. A nonlinear objective function $J{\mathbf{\theta }}^{\mathrm{\Delta }}$ is established based on the measured data ${\mathbf{e}}_m^k$ and ${\mathbf{y}}_m^k$ of a single finite time task, and the IV-GNI method is proposed to identify unbiased feedforward parameters and delay time with the presence of disturbance. Then, the feedforward parameters are tuned iteratively to match the tuned feedforward controller with the inverse model of the plant, so as to minimize the position error and achieve the optimal feedforward control performance.

3. Realization Method of Feedforward Control with Precise Time-Delay Compensation

Discrete signals are utilized in the digital control system, and the control period is $T_s$. The realization method of feedforward control is explained by taking $n=2$ as an example in equation (1). The feedforward control signal is ${\mathbf{u}}_{ff}=\mathbf{r}{e}^{\tau s}{\theta }_2{s}^2+{\theta }_{1}s+{\theta }_0={\mathbf{a}}_{ff}{\theta }_2+{\mathbf{v}}_{ff}{\theta }_1+{\mathbf{r}}_{ff}{\theta }_0$, where ${\mathbf{r}}_{ff}$, ${\mathbf{v}}_{ff}$, and ${\mathbf{a}}_{\text{ff}}$ are the position, velocity, and acceleration inputs in the feedforward control, respectively, and $\tau$ is the delay time. The control block diagram is shown in Figure 4. At time $t$, the system reference position is the planned signal $rt$, and ${\mathbf{r}}_{ff}$, ${\mathbf{v}}_{ff}$, and ${\mathbf{a}}_{ff}$ are $rt+\tau$, $vt+\tau$, and $at+\tau$, respectively, which are ahead of time period $\tau$. Then, they are introduced into the feedforward model to generate the feedforward force, which can make the plant reach the reference position, and thus realize the feedforward control with time-delay compensation, so as to improve the control performance.

[figure omitted; refer to PDF]

The ahead planning of the position ${\mathbf{r}}_{ff}$ by accurate path planning time-shift is explained below. The delay time $\tau$ may be a nonintegral multiple of period, which satisfies $n_0{T}_s\le \tau <n_0+1T_s$ and $n_0$ is a nonnegative integer. As shown in Figure 5(a), the continuous time position $\mathbf{r}$ (the thick black line) is obtained by the position planning function [29] and it is shifted to the left by $\tau$ along the time axis to obtain ${\mathbf{r}}_{ff}$ (the thick blue line), which realizes ${\mathbf{r}}_{ff}=\mathbf{r}{e}^{\tau s}$. At the sampling time $t=n{T}_s,n=\mathrm{0,1,2},\dots ,$, ${\mathbf{r}}_d$ (the fine black line) is $rt=rn{T}_s$, and ${\mathbf{r}}_{ffd}$ (the fine blue line) is $rt+\tau =rn{T}_s+\tau$, which are the discrete position signals received by the system. In fact, the desired signal processed by the actual lead-time term is ${\mathbf{r}}_{\text{lead}}$ (the dotted red line), i.e., shifting ${\mathbf{r}}_d$ to the left by $\tau$, but each point of ${\mathbf{r}}_{\text{lead}}$ is not at the sampling time, which cannot be obtained. Here, there are errors in approximating ${\mathbf{r}}_{\text{lead}}$ with ${\mathbf{r}}_{ffd}$, but since the system runs periodically, the errors cannot be eliminated. In [30], the position signal with lead time of integer multiple of period is obtained first and then delayed by means of calculation to realize the lead-time of noninteger multiple of period. In comparison, the method by path planning time-shift does not need delay calculation, which is more convenient to implement in one step. The lead or lag is a relative relationship in time and can be converted into each other. Thus, ${\mathbf{r}}_{ff}=\mathbf{r}{e}^{\tau s}$ can be realized by the delay planning of $\mathbf{r}$ relative to ${\mathbf{r}}_{ff}$; that is, $\mathbf{r}={\mathbf{r}}_{ff}{e}^{-\tau s}$. As shown in Figure 5(b), at the sampling time $t$, ${\mathbf{r}}_{\text{ff}}$ is planned as $rn{T}_s$ and $\mathbf{r}$ is planned as $rn{T}_s-\tau$. Thus, the lag of the reference position signal is realized, and the same effect can be achieved. The lead time of the feedforward signal ${\mathbf{a}}_{ff}$ and ${\mathbf{v}}_{ff}$ relative to the reference signal $\mathbf{r}$ can also be implemented by this method to complete the feedforward control.

[figures omitted; refer to PDF]

4. Simulation Analysis

In this section, the simulation is implemented to verify whether the feedforward model obtained by the proposed algorithm matches the inverse of the controlled plant. Three controlled plants with time-delay, a mass model, a mass damping model, and a mass stiffness damping model, are considered separately in the simulation, and the feedback-feedforward control system is shown in Figure 1. The feedforward model is established, and the initial parameters are given. Then, the Gauss–Newton iterative method based on instrumental variable is used to identify and tune the feedforward parameters with the measurement data of the finite time task. Finally, the feedforward control model with delay compensation is introduced into the system for control performance improvement. Meanwhile, the parameter accuracy of the proposed algorithm is compared with several existing tuning algorithms.

4.1. Simulation 1

In simulation 1, a mass model with time-delay is considered and the plant is given by$$\begin{array}{}\text{(20)}& P=\frac{1}{s^2}{e}^{-0.0017s}.\end{array}$$

It is introduced into the control loop, and the control period is set as 1 ms. The feedback controller is a PID controller, and $C_{fb}=k_{p}1+2\pi f_i/s+s/2\pi f_d$, where $k_p=22000.0$, $f_i=5.0\text{Hz}$, and $f_d=10.0\text{\hspace{0.17em}Hz}$, which enables the system with bandwidth of 129.6 Hz, phase margin of 35.6°, and amplitude margin of 6.3 dB.

4.1.1. Simulation Illustration of Feedforward Parameter Tuning Method

Firstly, according to the parameterization method in Section 2 and the model of plant, the feedforward controller is parameterized as $C_{ff}={\theta }_2{s}^2{e}^{\tau s}$. The initial feedforward controller is $C_{ff}^0=0.5s^2{e}^{0.0006s}$, and the initial parameters are ${\theta }_2^0=0.5$ and ${\tau }^0=0.0006$, which are given based on the rule that it should be smaller than the plant parameters according to Appendix.

In a finite time task, the system reference position is a third-order point-to-point motion path, as shown in Figure 6. Since $C_{ff}^0$ contains $e^{0.0006s}$, $\mathbf{r}$ (the solid black line) is delayed by 0.6 ms relative to ${\mathbf{r}}_{ff}$ (the dashed blue line) by the method of path planning in Section 3, so the acceleration ${\mathbf{a}}_{ff}$ in feedforward control is advanced relative to $\mathbf{r}$. Then, the control system with acceleration and time-delay feedforward is implemented.

[figures omitted; refer to PDF]

The task is executed, and the feedforward correction parameters are identified based on the measured data. Since the key step is the calculation of the feedforward correction parameters ${\mathbf{\theta }}^{\mathrm{\Delta }}$ in the proposed tuning method and the existing IFT methods, ${\mathbf{\theta }}^{\mathrm{\Delta }}$ is calculated by IV-GNI in the proposed method and by NI and IV-LS in the existing IFT, respectively. Then, the results are compared to elaborate the applicability of these methods.

In order to verify the ability of the proposed identification algorithm to tolerate the disturbance, the simulation analysis is carried out in two cases:

(1) In the case of disturbance $\mathbf{w}=0$

  No disturbance is introduced into the system, and ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ (the black line in Figure 7(a)) is obtained by IV-GNI. For comparison, the Gauss–Newton iteration method is adopted to calculate $\mathrm{\Delta }{\mathbf{\theta }}_m^{\mathrm{\Delta }}$ directly by the measurement data and ${\mathbf{\theta }}_{\text{GNI}}^{\mathrm{\Delta }}$ (the blue line) is obtained by equation (14). Then, ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$ (the magenta line) and ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ (the green line) are obtained by the Newton iterative method and the least square method based on instrumental variable in existing IFT. Since the linear model is established in IV-LS, the delay time cannot be calculated. The results are also shown in Table 1.

(2) In the case of $\mathbf{w}\ne 0$

  The disturbance $\mathbf{w}=Hs\mathbf{\epsilon }$ is introduced into the system. Similar to the disturbance as present in [11, 12], $\mathbf{\epsilon }$ is a normally distributed white noise with zero mean and variance ${\lambda }_{\epsilon }^2=0.02$, and $Hs=0.05s+1/10s+1$. ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$, ${\mathbf{\theta }}_{\text{GNI}}^{\mathrm{\Delta }}$, ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$, and ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ are calculated, and the results are shown in Figure 7(b) and Table 2.

[figures omitted; refer to PDF]

Table 1

Iterative calculation results of feedforward parameters correction without disturbance in simulation 1.

	${\theta }_2^{\mathrm{\Delta }}$ (1 × 10−1)				${\tau }^{\mathrm{\Delta }}$ (1 × 10−3)
Iterations	IV-GNI	GNI	NI	IV-LS	IV-GNI	GNI	NI
1	5.00	5.00	5.02	4.83	1.08	1.10	0.216
2	5.00	5.00	5.00	4.83	1.09	1.10	1.09
3	5.00	5.00	5.00	4.83	1.09	1.10	1.10
4	5.00	5.00	5.00	4.83	1.09	1.10	1.10
5	5.00	5.00	5.00	4.83	1.09	1.10	1.10

Table 2

Iterative calculation results of feedforward parameters correction under disturbance in simulation 1.

	${\theta }_2^{\mathrm{\Delta }}$ (1 × 10−1)				${\tau }^{\mathrm{\Delta }}$ (1 × 10−3)
Iterations	IV-GNI	GNI	NI	IV-LS	IV-GNI	GNI	NI
1	5.01	4.79	4.94	4.84	1.09	9.13	1.22
2	5.00	4.79	4.86	4.84	1.10	9.29	5.75
3	5.00	4.79	4.85	4.84	1.10	9.29	5.84
4	5.00	4.79	4.85	4.84	1.10	9.29	5.84
5	5.00	4.79	4.85	4.84	1.10	9.29	5.84

From the above results, we can get the following: (1) when $\mathbf{w}=\mathbf{0}$, each value in ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ converges to the corresponding correct value, ${\theta }_2^{\mathrm{\Delta }}=0.5$ and ${\tau }^{\mathrm{\Delta }}=0.0011$, which are the difference values between the plant parameters and the initial feedforward parameters. It is shown that the optimization scheme in Section 2 is effective and the IV-GNI method is also applicable in the case without disturbance. As the path planning time-shift method is employed to realize $e^{0.0006s}$ in $C_{ff}^0$, its validity is indirectly proved. ${\mathbf{\theta }}_{\text{GNI}}^{\mathrm{\Delta }}$ and ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$ are equal to the accurate values, indicating that GNI and NI can also be used in this case. The stable value of ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ is 0.483, the deviation from its accurate value is 0.017, and the delay time cannot be calculated, which shows that IV-LS is not suitable for the identification of the nonlinear system with time-delay. (2) When $\mathbf{w}\ne 0$, the stable values of ${\theta }_{2\_\text{GNI}}^{\mathrm{\Delta }}$, ${\theta }_{2\_\text{NI}}^{\mathrm{\Delta }}$, and ${\theta }_{2\_\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ are 0.479, 0.485, and 0.484, and the deviations from the accurate value are 0.021, 0.015, and 0.016. The stable values of ${\tau }_{\text{GNI}}^{\mathrm{\Delta }}$ and ${\tau }_{\text{NI}}^{\mathrm{\Delta }}$ are 0.00929 and 0.00584, and the deviations from the accurate value are 0.00819 and 0.00474. The stable values of ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ are equal to the accurate values. It shows that GNI and NI cannot deal with the disturbance, which will lead to the identification deviation. The stable values of ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ with and without disturbances are basically equal, which means that the instrumental variable eliminates the effect of disturbance, but IV-LS cannot calculate the delay time. The introduction of IV in GNI is feasible, which has the ability to tolerate the disturbance, and the unbiased estimation of ${\mathbf{\theta }}^{\mathrm{\Delta }}$ can be obtained.

Finally, the feedforward parameters are tuned. The feedforward controller is updated to $C_{ff}^1=s^2{e}^{0.0017s}$ according to equation (5) with ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$, which is equal to the inverse model of the plant in equation (20), so the feedforward tuning is completed and the effectiveness of the feedforward parameter tuning method is verified.

4.1.2. Feedforward Control Results in Simulation 1

Based on the realization method of the feedforward control in Section 3, the tuned feedforward controller $C_{ff}^1$ is introduced into the system to obtain the position error without disturbance $\mathbf{w}$ injection.

In Figure 8, the position error of the feedback control is shown by the black line, and the position error after introducing acceleration as well as velocity feedforward is shown by the dashed blue line; it can be seen that the position error decreases obviously after the feedforward control is introduced and the maximum error is reduced from $3.292\times {10}^{-4}\text{\hspace{0.17em}m}$ to $4.172\times {10}^{-5}\text{\hspace{0.17em}m}$.

[figure omitted; refer to PDF]

In Figure 9, the position error without the time-delay compensation in the feedforward control is shown by the black line, and the position error with the time-delay compensation is shown by the dashed blue line. It can be seen that the tracking accuracy is further improved by introducing the time-delay compensation. The error of nonzero jerk segment is greatly reduced, and the maximum error is reduced from $4.172\times {10}^{-5}$ m to $4.944\times {10}^{-6}$ m, which means the impact effect caused by jerk segment is mitigated and the time-delay feedforward compensation method in Section 3 is effective.

[figure omitted; refer to PDF]

In Figure 13, the tracking accuracy is further improved by introducing the time-delay compensation. The error of nonzero jerk segment is greatly reduced, and the maximum error is reduced from $5.231\times {10}^{-5}\text{\hspace{0.17em}m}$ to $7.707\times {10}^{-6}\text{\hspace{0.17em}m}$, which means the impact effect caused by jerk segment is mitigated and the time-delay feedforward compensation method in Section 3 is effective.

[figure omitted; refer to PDF]

In Figure 17, the tracking accuracy is further improved by introducing the time-delay compensation. The error of nonzero jerk segment is greatly reduced, and the maximum error is reduced from $2.36\times {10}^{-5}\text{\hspace{0.17em}m}$ to $1.338\times {10}^{-6}\text{\hspace{0.17em}m}$, which means the impact effect caused by jerk segment is mitigated and the time-delay feedforward compensation method in Section 3 is effective.

[figure omitted; refer to PDF]

The frequency response curve of the controlled plant is shown by the fine black line in Figure 19. In the frequency range below 100 Hz, the amplitude-frequency characteristic indicates that the controlled plant can be identified as a second-order rigid-body model approximately. The phase-frequency characteristic does not remain −180° but decreases linearly with the increase in the frequency and the slope is constant, which is considered to be caused by the time-delay [23].

[figures omitted; refer to PDF]

The feedback controller is a PID controller, and $C_{fb}=k_{p}1+2\pi f_i/s+s/2\pi f_d$, where $k_p=900.0$, $f_i=5.0\text{\hspace{0.17em}Hz}$, and $f_d=10.0\text{\hspace{0.17em}Hz}$, which enables the system with a bandwidth of 64.7 Hz, phase margin of 29.6°, and amplitude margin of 5.6 dB.

5.2. Experimental Verification of Feedforward Parameter Tuning Method

Firstly, the feedforward control model is established. The air floating guide in the platform eliminates the influence of friction, and the plant is a second-order rigid-body motion model with time-delay, so the feedforward controller is parameterized as $C_{ff}={\theta }_2{s}^2{e}^{\tau s}$.

In the finite time task, the system reference position is a third-order point-to-point motion path, as shown in Figure 20, and the maximum values for position, velocity, acceleration, and jerk are 0.1 m, 0.1 m/s, 1 m/s², and 10.0 m/s³.

[figures omitted; refer to PDF]

Then, the feedforward correction parameters are calculated. The initial feedforward controller is set to $C_{ff}^0=0$, so the initial task is executed without the feedforward control. ${\mathbf{e}}_m^0$ and ${\mathbf{y}}_m^0$ are measured, and the objective function is established. The same as the simulation, ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$, ${\mathbf{\theta }}_{\text{GNI}}^{\mathbf{\Delta }}$, ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$, and ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ are calculated by four algorithms, respectively. The results are shown in Figure 21 and Table 7. In ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$, ${\theta }_{2\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ converges to 0.080 and ${\tau }_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ converges to 0.0029.

[figures omitted; refer to PDF]

Table 7

Iterative calculation results of feedforward parameters correction in experiment.

	${\theta }_2^{\mathrm{\Delta }}$ (1 × 10−2)				${\tau }^{\mathrm{\Delta }}$ (1 × 10−3)
Iterations	IV-GNI	GNI	NI	IV-LS	IV-GNI	GNI	NI
1	8.0	7.4	8.1	8.4	2.70	6.84	0.015
2	8.0	7.4	7.6	8.4	2.90	7.17	5.59
3	8.0	7.4	7.5	8.4	2.90	7.17	6.25
4	8.0	7.4	7.5	8.4	2.90	7.17	6.37
5	8.0	7.4	7.5	8.4	2.90	7.17	6.37

Finally, the feedforward parameters are tuned. With the contribution of ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$, the new feedforward controller is obtained as $C_{\text{ff}}^1=0.08s^2{e}^{0.0029s}$. The inverse model of $C_{ff}^1$ is $1/0.08s^2{e}^{-0.0029s}$, and the frequency response curve is shown as the thick blue line in Figure 19, which is consistent with the plant. So the tuning of the feedforward parameters is completed, and the stable values of ${\mathbf{\theta }}_{\text{IV}-\text{GNI}}^{\mathrm{\Delta }}$ are correct. The stable values of ${\theta }_{2\_\text{GNI}}^{\mathrm{\Delta }}$, ${\theta }_{2\_\text{NI}}^{\mathrm{\Delta }}$, and ${\theta }_{2\_\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ are 0.074, 0.075, and 0.084, and the stable values of ${\tau }_{\text{GNI}}^{\mathrm{\Delta }}$ and ${\tau }_{\text{NI}}^{\mathrm{\Delta }}$ are 0.00717 and 0.00637, which means that the results of ${\mathbf{\theta }}_{\text{GNI}}^{\mathrm{\Delta }}$, ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$, and ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ are biased. The cable forces and other unknown disturbances during the experiment are not dealt within GNI and NI, which leads to the bias of ${\mathbf{\theta }}_{\text{GNI}}^{\mathrm{\Delta }}$ and ${\mathbf{\theta }}_{\text{NI}}^{\mathrm{\Delta }}$. The results of ${\mathbf{\theta }}_{\text{IV}-\text{LS}}^{\mathrm{\Delta }}$ indicate that IV-LS is not applicable to nonlinear systems, and the identification accuracy of dynamic parameters is affected. The GNI based on IV can achieve unbiased estimation under disturbance and can be applied to feedforward parameter tuning with time-delay.

5.3. Experimental Results

According to the realization method of feedforward control with precise time-delay compensation, the new feedforward controller $C_{ff}^1$ is introduced into the control system, and the following experiments are carried out to obtain the position error.

In Figure 22, the position error under the feedback control is the black line and the position error under the feedback and acceleration feedforward control is the dashed blue line, we can see that the tracking performance of nonuniform velocity section is improved obviously under the acceleration feedforward control, and the maximum error is reduced from $3.085\times {10}^{-5}\text{\hspace{0.17em}m}$ to $4.028\times {10}^{-6}\text{\hspace{0.17em}m}$. In Figure 23, the feedback control signal ${\mathbf{u}}_{fb}$ under the feedback control is shown as the black line. In Figure 24, the feedback control signal ${\mathbf{u}}_{fb}$ and the feedforward control signal ${\mathbf{u}}_{ff}$ under the feedback and acceleration feedforward control are shown as the black line and blue line, respectively. After the introduction of the feedforward control, the thrust required for system motion is provided by the sum of ${\mathbf{u}}_{fb}$ and ${\mathbf{u}}_{ff}$. As shown in Figure 24, ${\mathbf{u}}_{ff}$ provides most of the thrust required in the acceleration section, so ${\mathbf{u}}_{fb}$ is greatly reduced in this region correspondingly. Since the position error is processed by the controller to get ${\mathbf{u}}_{fb}$, the position error is also greatly reduced.

[figure omitted; refer to PDF]

Based on the acceleration feedforward, the position errors with and without the time-delay feedforward are shown as the black line and dashed blue line in Figure 25, respectively, and the feedback control signal ${\mathbf{u}}_{fb}$ is shown in Figure 26. Since the time-delay feedforward is realized by the delay planning of the reference position, the feedforward control signal ${\mathbf{u}}_{ff}$ under the two feedforward control are the same and shown by ${\mathbf{u}}_{ff}$ in Figure 24. It can be seen that the position error caused by the impact of jerk segment is further reduced under the time-delay feedforward and the maximum error is reduced from $4.028\times {10}^{-6}\text{\hspace{0.17em}m}$ to $2.057\times {10}^{-6}\text{\hspace{0.17em}m}$, and the fluctuation of ${\mathbf{u}}_{fb}$ decreases in this region, which shows that the time-delay compensation method in Section 3 is effective. We can see that ${\mathbf{u}}_{fb}$ still exists after the introduction of acceleration and time-delay feedforward, which should be the compensation signal generated by the control loop for the cable force and other disturbances in the system.

[figure omitted; refer to PDF]

The experimental results show that the proposed method can be applied in the precision motion control system with unknown time-delay. The feedforward parameters and delay time can be identified by IV-GNI based on the measurement data of a single finite time task without detailed knowledge of the plant, and the time-delay can be compensated with path planning time-shift, which can obviously improve the control performance of the system.

6. Conclusions

This paper proposes an iterative tuning method of feedforward parameters considering the time-delay. The key essentials of the proposed method lie in the following: (1) a nonlinear objective function with time-delay suitable for iterative feedforward tuning is established; (2) the Gauss–Newton iterative method and instrumental variables are combined to optimize the nonlinear objective function, and the unbiased estimation of the feedforward parameters and delay time are obtained in the presence of disturbance; and (3) the precise time-delay compensation can be realized in feedforward control with accurate path planning time-shift. The Newton iterative method and the least square method based on instrumental variable in existing IFT and the conventional GNI are compared with the proposed IV-GNI scheme. The identification results of the feedforward parameters in the simulation and experiment on an air floating precision motion platform show that the proposed method can be applied to the nonlinear system with time-delay and realize the unbiased parameter estimation with disturbance and noise. The results of feedforward control in simulation and experiment prove that the control performance of the precision system can be obviously improved with the proposed method.

This method can be extended to precision motion control MIMO systems, the decoupling method is needed to make the coupling MIMO become decoupled multiple SISO systems, and the corresponding controlled plant should have a mathematical model without zero.

Acknowledgments

The authors deeply acknowledge the financial support from the National Science and Technology Major Project of China (no. 2017ZX02101007-002) and the National Natural Science Foundation of China (nos. 51675195, 51705163, 51721092, and 51975234). The authors also thank Yan Zhu in Micro and Nano Fabrication and Measurement Laboratory for the support in experimental verification.

Appendix

Appendix

Proof of convexity of the objective function

The controlled plant is assumed to be $P=e^{-{\tau }^{\ast }s}/{\displaystyle {\sum }_{i=0}^n{\theta }_i^{\ast }{s}^i}$ and ${\theta }_i^{\ast }$, and ${\tau }^{\ast }$ are the true values. In the tuning procedure, the feedforward parameters should satisfy ${\theta }_i^k,{\tau }^k\in 0\le {\theta }_i^k\le {\theta }_i^{\ast },0\le {\tau }^k\le {\tau }^{\ast }$. A new objective function is established by the polynomial on the right side of equation (8):$$\begin{array}{}\text{(A.1)}& J_1{\mathbf{\theta }}^{\Delta }={\displaystyle \sum }_{j=1}^N{rj\text{GP}{C}_{\text{fb}}+C_{\text{ff}}^k{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^{n+1}}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}-{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^k+{\theta }_i^{\Delta }{s}^i}{s^{n+1}}}^2,\end{array}$$which is equivalent to the original objective function $J{\mathbf{\theta }}^{\mathrm{\Delta }}$. In the interval $D={\theta }_i^{\mathrm{\Delta }},{\tau }^{\mathrm{\Delta }}{-\theta }_i^k\le {\theta }_i^{\mathrm{\Delta }}\le {\theta }_i^{\ast }-{\theta }_i^k,-{\tau }^k\le {\tau }^{\mathrm{\Delta }}\le {\tau }^{\ast }-{\tau }^k$, $J_1{\mathbf{\theta }}^{\mathrm{\Delta }}$ is differentiable. When choosing and ${\theta }_i^{\mathrm{\Delta }}+\mathrm{\Delta }{\theta }_i^{\mathrm{\Delta }},{\tau }^{\mathrm{\Delta }}+\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}\in D$, we can obtain$$\begin{array}{}\text{(A.2)}& J_1{\mathbf{\theta }}^{\Delta }+\Delta {\mathbf{\theta }}^{\Delta }-J_1{\mathbf{\theta }}^{\Delta }=-{\displaystyle \sum _{j=1}^N\begin{array}{c}rjGP{C}_{fb}+C_{ff}^k\begin{array}{c}{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}\frac{1}{s^2}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}-e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }+\Delta {\tau }^{\Delta }s}\\ +{\displaystyle \sum }_{i=0}^n\frac{\Delta {\theta }_i^{\Delta }{s}^i}{s^{n+2}}\end{array}\cdot \\ rjGP{C}_{fb}+C_{ff}^k\begin{array}{c}{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}+e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }+\Delta {\tau }^{\Delta }s}\\ -{\displaystyle \sum }_{i=0}^n\frac{2{\theta }_i^k+2{\theta }_i^{\Delta }+\Delta {\theta }_i^{\Delta }{s}^i}{s^n}\end{array}\end{array}}.\end{array}$$

The derivative of $J_1{\mathbf{\theta }}^{\mathrm{\Delta }}$ to ${\mathbf{\theta }}^{\mathrm{\Delta }}$ is multiplied by $\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}$ to get$$\begin{array}{}\text{(A.3)}& \nabla J_1{\mathbf{\theta }}^{\Delta }\cdot \Delta {\mathbf{\theta }}^{\Delta }=-{\displaystyle \sum _{j=1}^N\begin{array}{c}rjGP{C}_{fb}+C_{ff}^k{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}\frac{1}{s^2}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}\cdot \Delta {\tau }^{\Delta }s+{\displaystyle \sum }_{i=0}^n\frac{\Delta {\theta }_i^{\Delta }{s}^i}{s^{n+2}}\cdot \\ rjGP{C}_{fb}+C_{ff}^k{\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}2e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}-{\displaystyle \sum }_{i=0}^n\frac{2{\theta }_i^k+{\theta }_i^{\Delta }{s}^i}{s^n}\end{array}}.\end{array}$$

In equations (A.2) and (A.3), we can obtain$$\begin{array}{}\text{(A.4)}& {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}2e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}-{\displaystyle \sum }_{i=0}^n\frac{2{\theta }_i^k+{\theta }_i^{\Delta }{s}^i}{s^n}\ge {\displaystyle \sum }_{i=0}^n\frac{{\theta }_i^{\ast }{s}^i}{s^n}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}+e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }+\Delta {\tau }^{\Delta }s}-{\displaystyle \sum }_{i=0}^n\frac{2{\theta }_i^k+2{\theta }_i^{\Delta }+\Delta {\theta }_i^{\Delta }{s}^i}{s^n}\ge 0.\end{array}$$

Using the first-order Taylor expansion with Lagrange remainder, $e^{-\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}s}=1-\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}s+s^2{e}^{-\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}\theta s}\cdot {\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}}^2/2$, $\theta \in \mathrm{0,1}$, so $1-e^{-\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}s}/s^2=\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}s/s^2-e^{-\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}\theta s}\cdot {\mathrm{\Delta }{\tau }^{\mathrm{\Delta }}}^2/2\le \mathrm{\Delta }{\tau }^{\mathrm{\Delta }}s/s^2$; we can get$$\begin{array}{}\text{(A.5)}& \frac{1}{s^2}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}\cdot \Delta {\tau }^{\Delta }s\ge \frac{1}{s^2}{e}^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }s}-e^{{\tau }^{\ast }-{\tau }^k+{\tau }^{\Delta }+\Delta {\tau }^{\Delta }s}\ge 0.\end{array}$$

Then, $J_1{\mathbf{\theta }}^{\mathrm{\Delta }}+\mathrm{\Delta }{\mathbf{\theta }}^{\mathrm{\Delta }}-J_1{\mathbf{\theta }}^{\mathrm{\Delta }}\ge \nabla J_1{\mathbf{\theta }}^{\mathrm{\Delta }}\cdot \mathrm{\Delta }{\theta }^{\mathrm{\Delta }}$ can be obtained, so $J_1{\mathbf{\theta }}^{\mathrm{\Delta }}$ and $J{\mathbf{\theta }}^{\mathrm{\Delta }}$ are the convex functions in the interval $D$. $D$ contains the extreme point ${\theta }_i^{\mathrm{\Delta }}={\theta }_i^{\ast }-{\theta }_i^k\text{\hspace{0.17em}and\hspace{0.17em}}{\tau }^{\mathrm{\Delta }}={\tau }^{\ast }-{\tau }^k$ that makes the derivative zero. In the iterative calculation, the initial value is ${\theta }_i^{\mathrm{\Delta }}=0,{\tau }^{\mathrm{\Delta }}=0\in D$, which can converge to the extreme point.

References

[1] J.-P. Richard, "Time-delay systems: an overview of some recent advances and open problems," Automatica, vol. 39 no. 10, pp. 1667-1694, DOI: 10.1016/s0005-1098(03)00167-5, 2003.

[2] A. San-Millan, D. Russell, V. Feliu, S. S. Aphale, "A modified positive velocity and position feedback scheme with delay compensation for improved nanopositioning performance," Smart Mater. Struct., vol. 24 no. 7,DOI: 10.1088/0964-1726/24/7/075021, 2015.

[3] H. Meng, X. Sun, J. Xu, F. Wang, "The generalization of equal-peak method for delay-coupled nonlinear system," Physica D: Nonlinear Phenomena, vol. 403,DOI: 10.1016/j.physd.2020.132340, 2020.

[4] F. Wang, X. T. Sun, H. Meng, J. Xu, "Time-delayed feedback control design and its application for vibration absorption," IEEE Transactions on Industrial Electronics, vol. 40, 2020.

[5] M. Tomizuka, "Zero Phase Error Tracking Algorithm for Digital Control," The Journal of Dynamic Systems, vol. 24 no. 1, pp. 65-68, 1987.

[6] F. Z. Song, Y. Liu, J. X. Xu, X. F. Yang, Q. Zhu, "Data-driven iterative feedforward tuning for a wafer stage: a high-order approach based on instrumental variables," IEEE Transactions on Industrial Electronics, vol. 66 no. 4, pp. 3106-3116, 2018.

[7] M. Heertjes, D. Hennekens, M. Steinbuch, "MIMO feed-forward design in wafer scanners using a gradient approximation-based algorithm," Control Engineering Practice, vol. 18 no. 5, pp. 495-506, DOI: 10.1016/j.conengprac.2010.01.006, 2010.

[8] E. D. Tung, M. Tomizuka, "Feedforward tracking controller design based on the identification of low frequency dynamics," Journal of Dynamic Systems, Measurement, and Control, vol. 115 no. 3, pp. 348-356, DOI: 10.1115/1.2899109, 1993.

[9] K. K. Leang, Q. Zou, S. Devasia, "Feedforward control of piezoactuators in atomic force microscope systems," IEEE Control Systems, vol. 29 no. 1, pp. 70-82, DOI: 10.1109/mcs.2008.930922, 2009.

[10] S. Devasia, "Should model-based inverse inputs be used as feedforward under plant uncertainty?," IEEE Transactions on Automatic Control, vol. 47 no. 11, pp. 1865-1871, DOI: 10.1109/tac.2002.804478, 2002.

[11] F. Boeren, T. Oomen, M. Steinbuch, "Iterative motion feedforward tuning: a data-driven approach based on instrumental variable identification," Control Engineering Practice, vol. 37, pp. 11-19, DOI: 10.1016/j.conengprac.2014.12.015, 2015.

[12] F. Boeren, T. Oomen, "Iterative Feedforward Control: A Closed-Loop Identification Problem and a Solution," Proceedings of the 52nd IEEE Conference on Decision and Control, pp. 6694-6699, .

[13] F. Boeren, D. Bruijnen, T. Oomen, "Enhancing feedforward controller tuning via instrumental variables: with application to nanopositioning," International Journal of Control, vol. 90 no. 4, pp. 746-764, DOI: 10.1080/00207179.2016.1219921, 2016.

[14] F. Boeren, T. Oomen, M. Steinbuch, "Accuracy Aspects in Motion Feedforward Tuning," Proceedings of the 2014 American Control Conference, pp. 2178-2183, .

[15] L. Li, Y. Liu, L. Li, J. Tan, "Kalman-filtering-based iterative feedforward tuning in presence of stochastic noise: with application to a wafer stage," IEEE Transactions on Industrial Informatics, vol. 15 no. 11, pp. 5816-5826, DOI: 10.1109/tii.2019.2906331, 2019.

[16] D. Bruijnen, N. V. Dijk, "Combined Input Shaping and Feedforward Control for Flexible Motion Systems," Proceedings of the 2012 American Control Conference, pp. 2473-2478, .

[17] F. Boeren, D. Bruijnen, N. van Dijk, T. Oomen, "Joint input shaping and feedforward for point-to-point motion: automated tuning for an industrial nanopositioning system," Mechatronics, vol. 24 no. 6, pp. 572-581, DOI: 10.1016/j.mechatronics.2014.03.005, 2014.

[18] F. Boeren, L. Blanken, D. Bruijnen, T. Oomen, "Optimal estimation of rational feedforward control via instrumental variables: with application to a wafer stage," Asian Journal of Control, vol. 20 no. 3, pp. 975-992, DOI: 10.1002/asjc.1595, 2018.

[19] W. Huang, K. Yang, Y. Zhu, X. Li, H. Mu, M. Li, "Data-driven rational feedforward tuning: with application to an ultraprecision wafer stage," Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 234 no. 6, pp. 748-758, DOI: 10.1177/0959651819874563, 2020.

[20] J. Na, X. Ren, Y. Xia, "Adaptive parameter identification of linear SISO systems with unknown time-delay," Systems & Control Letters, vol. 66, pp. 43-50, DOI: 10.1016/j.sysconle.2014.01.005, 2014.

[21] H. Stearns, S. Mishra, M. Tomizuka, "Iterative Tuning of Feedforward Controller with Force Ripple Compensation for Wafer Stage," Proceedings of the 2008 10th IEEE International Workshop on Advanced Motion Control, pp. 234-239, .

[22] S. V. Meulen, R. Tousain, O. Bosgra, "Fixed structure feedforward controller design exploiting iterative trials: application to a wafer stage and a desktop printer," The Journal of Dynamic Systems, vol. 130, 2008.

[23] F. Z. Song, Y. Liu, J. X. Xu, X. F. Yang, P. He, Z. L. Yang, "Iterative learning identification and compensation of space-periodic disturbance in PMLSM systems with time delay," The Journal of Dynamic Systems, vol. 65 no. 9, pp. 7579-7589, 2017.

[24] X. M. Ren, A. B. Rad, P. T. Chan, W. L. Lo, "Online identification of continuous-time systems with unknown time delay," IEEE Transactions on Automatic Control, vol. 50 no. 9, pp. 1418-1422, DOI: 10.1109/tac.2005.854640, 2005.

[25] R. W. Jones, Recursive Estimation of Time-Delay Using a High-Order Pade Approximant, 2004.

[26] F.-w. Chen, T. Liu, "Iterative identification of discrete-time output-error model with time delay," Journal of Central South University, vol. 24 no. 3, pp. 647-654, DOI: 10.1007/s11771-017-3465-1, 2017.

[27] L. Lennart, System Identification: Theory for the User, 1999.

[28] T. C. Hsia, S. Q. Liu, K. Y. Zhu, System Identification, Least Square Method, 1983.

[29] P. Lambrechts, M. Boerlage, M. Steinbuch, "Trajectory planning and feedforward design for electromechanical motion systems," Control Engineering Practice, vol. 13 no. 2, pp. 145-157, DOI: 10.1016/j.conengprac.2004.02.010, 2005.

[30] W. Teng, Y. F. Zhou, X. Q. Li, L. L. Yang, Design of Acceleration Feedforward and Study of Calibration Algorithm, 2008.