Dimensional Error Compensation Based on Data-Driven Sliding Mode Terminal Iterative Learning Control for CNC Batch Grinding

Abstract

Dimensional error of workpiece is an important index that indicts the grinding accuracy of machine tools. Conventional error compensation strategies focus on modeling one or several error factors for a particular CNC machining type, which are time consuming and costly. Based on data-driven sliding mode terminal iterative learning control, a new dimensional error compensation method is proposed in this paper. The compensation algorithm based on an iterative sliding surface and an objective function is established using the terminal experimental data. A modified NC program is fed to the machine tool to push the workpiece dimension towards desired. The theoretical analysis and experimental results show the proposed compensation method can effectively improve the batch grinding accuracy of the machine tool, which shows good perspective in manufacturing industry.

1. Introduction

As the manufacturing technology has continuously been developed, the accuracy requirements of computer numerical control (CNC) machining are becoming increasingly important. Error compensation has become a hot issue for actively improving the machining accuracy, which is influenced by error sources such as geometric error of machine tools, thermal error, cutting force, tool’s wear error, etc. [1,2,3]. Software compensating is a convenient solution to error compensation, which consists of three crucial steps, measurement strategy, identification approach, and compensation method. It improves machining accuracy and surface quality of workpiece through regenerating numerical control (NC) program for next compensation cycle, which is characterized by high efficiency and low cost.

There are many studies that have focused on one error source to compensate, which is the main cause of machine inaccuracy under certain machining. For example, the geometric error makes up a significant portion of the inaccuracy of a machine tool, accounting for more than 30% [4]. The authors of [5] presented a novel 13-line identification algorithm by employing a laser interferometer, and by correcting compensation parameters in NC instructions the linear axes are compensated. In [6], by using double ball bar (DBB), a comprehensive geometric error measurement and identification method is presented. However, for the experiments, on a tilt table of five-axis machine tool, the authors assumed that the geometric errors of the translational axes were abandoned. Wu et al. [7] presented an integrated error prediction and compensation method for the translational axes’ geometric errors. With this integrated approach, the positioning inaccuracy of tool ball for a double ball bar is greatly decreased in a multi-axis machine tool.

On the other hand, there are many studies reporting that the error caused by thermal error source constitutes another vital portion, up to 25–35% of total errors [8,9,10]. There are three categories of methods to eliminate thermal errors in current research, which are structural optimization of machine tool featured by thermal characteristics, temperature control, and thermal error compensation. The first two kinds are error prevention, which are high-cost. There is a large quantity of research focused on thermal error compensation for machining accuracy. On the basis of least square support vector machine (LS-SVM), a thermal error model of motorized spindle is proposed in [11]. The results showed that the predictive ability of this model was high. According to temperature variation, a thermal growth model of spindle is presented by Liu in [12], and a parameter identification method of the model is also given in this article. At present, almost all of the thermal error compensation method conducted on the machine tools is based on modeling error sources. On the contrary, the constructed various thermal error models are rarely applied to structural optimization or design of machine tools in no matter what manner [13].

Cutting force is another factor influencing machining accuracy of workpiece, which contributes about 10% of inaccuracy. Many researchers have focused on dealing with the reduction in tool deflection error induced by cutting force [14,15] and this kind compensation has widely been conducted in milling process.

Moreover, in practical application, after a long time machining, the tool’s wear increases, which is inevitable. Therefore, it is difficult to ensure that the precision of machine tools is still within acceptable limits.

From the above literature review and to our best knowledge of the research, it can be seen that most of the studies have concentrated on measuring and modeling of error sources, as well as parameter identification. As known to all, there are various error sources on a machine tool, not just one. In order to improve the machining precision more effectively, the error compensation system needs to possess the capability to deal with this combination of error sources rather than a single one [16]. It is obvious that the complexity of this work would be multiplied.

Apart from the above error compensation methods, there is another kind of method, which actively compensate errors from data analysis. For example, dimensional error of workpiece is the result of various error sources. The authors of [17] presented an intelligent pre-compensation method for dimensional error. A PD-type rule library is set up to compute the compensation of next workpiece according to the previous workpieces’ dimensions. Similarly, a P-type self-learning updating law is proposed to compute compensation according to the previous dimension error in [18]. Both methods have set up the iterative updating law with the previous three-dimensional errors or just the previous one. The form of the two updating laws is similar to proportion integration differentiation (PID) control law.

Our work focuses on the problem of dimensional error compensation in real CNC batch grinding, which is a meaningful and practical issue. For the data-driven methods, especially sliding mode terminal iterative learning methods, it is the first time they have been introduced in a CNC error compensation field. As the above mentioned, error sources are various and complex; furthermore, they are coupled to each other. Completely separating every error source is unrealistic. It is difficult to sense and model all error sources and model parameter identification is not easy. Its costs, including time and money, are very high. In CNC batch grinding, workpieces are grinded one by one, and terminal dimension data are easily obtained. With the data-driven compensation method, only one grating displacement sensor is used to measure the dimension, which does not need so many various sensors. Moreover, it does not need complex modeling and parameter identification. In practical application, the compensation is also easier to implement. Whether or not compensation is needed is determined before the workpiece machined. If it needs compensation, the compensation module is introduced and the compensated NC program is executed. Thus, it can be seen that compared with conventional compensation methods, data-driven compensation method has numerous advantages.

The terminal iterative learning control (TILC) as a data-driven algorithm can converge to the desired output. It does not need to track the entire output trajectory in the finite time, only the terminal output needs to be tracked [19,20]. Moreover, in view of its simple algorithm, strong robustness, and fast response, the sliding mode control (SMC) has received more and more attention [21,22]. Therefore, it is necessary to apply the two control strategies into batch grinding to improve dimension precision.

This paper presents a dimensional error compensation method based on sliding mode terminal ILC under the data-driven design framework. Different from the upper methods in [17,18], this method introduces a sliding surface, which is designed by the previous two-dimensional error. With the sliding surface, the dimensional error will quickly converge to zero after a number of iteration cycles. Consider the batch grinding system as a single input single output (SISO) nonlinear discrete-time system. Firstly, the nonlinear system is transferred into a linear system by using the dynamic linearization (DL) method [23]. Secondly, under the defined sliding surface, the terminal iterative learning law [24,25,26] is derived for the compensation system. Measurements of dimensions are used to compute compensation value according to the learning law. Finally, the NC program is modified to compensate the subsequent machining cycle. With this data-driven compensation method, the dimensional accuracy will be improved. The organization of the paper is shown as follows. Section 2 presents the grinding protocol and problem formulation. Section 3 presents the sliding mode terminal iterative learning compensation method. Simulation and experimental verification are drawn in Section 4. Finally, the paper is concluded in Section 5.

2. Grinding Protocol and Problem Formulation

2.1. Grinding Model

The structure of the CNC batch grinding machine tool is shown in Figure 1. As shown, the spindle drives the grinding wheel to rotate. The C axis is the rotation axis, which controls the workpiece to rotate. The B axis drives the workbench to rotate. The X axis is the linear axis, which is the feed shaft of the grinding wheel carriage. The Y axis is also the linear axis, which drives the grinding wheel carriage side to side. The peripheral surface of workpiece is machined through the combination of rotation axis and linear axis.

The motion model is shown in Equation (1):

(1)$\{\begin{array}{c}\mathsf{\omega }=\mathrm{l}(\mathsf{\theta })\\ \mathrm{X}=\mathrm{h}(\mathsf{\theta })\end{array}$

where $\mathsf{\theta }$ is the rotation angle of workpiece; $\mathsf{\omega }$ is the angular speed of workpiece rotation; and $\mathrm{X}$ is linear displacement of feed movement for grinding wheel carriage.

In the point of processing technology, grinding amount and grinding force are key factors influencing grinding quality. The grinding amount is indicated in terms of equivalent grinding thickness ${\mathrm{a}}_{\mathrm{eq}}$. It is shown in Equation (2):

(2)$\{\begin{array}{c}{\mathrm{a}}_{\mathrm{eq}}={\mathrm{a}}_{\mathrm{p}}\frac{{\mathrm{v}}_{\mathsf{\omega }}}{{\mathrm{v}}_{\mathrm{s}}}\\ {\mathrm{F}}_{\mathrm{t}}=\mathrm{k}{\mathrm{a}}_{\mathrm{p}}^{{\mathsf{\alpha }}_1}{\left(\frac{{\mathrm{v}}_{\mathsf{\omega }}}{{\mathrm{v}}_{\mathrm{s}}}\right)}^{{\mathsf{\alpha }}_2}\end{array}$

Here, ${\mathrm{a}}_{\mathrm{p}}$ is the grinding depth, which is in positive proportion to grinding wheel feed rate ${\mathrm{v}}_{\mathrm{r}}$; ${\mathrm{v}}_{\mathsf{\omega }}$ is speed of grinding point on the workpiece; ${\mathrm{v}}_{\mathrm{s}}$ is the speed of grinding point on the grinding wheel; and $\mathrm{k},{ \mathsf{\alpha }}_1,{ \mathsf{\alpha }}_2$ are parameters.

In order to improve the grinding quality, the linear speed of grinding point kept constant is introduced into processing technology. With the constant linear speed, the stock removal rate of workpiece and the grinding force would change stably. Therefore, the contour accuracy of workpiece is improved. On the other hand, with the constant linear speed, the linear speed of the grinding point on the machined peripheral surface is identical, so the consistency of the surface grinding quality can be ensured.

Take the square workpiece for example, the constant linear speed grinding model shown in Figure 2 is established. In Figure 2, point O1 is the rotation center of grinding wheel, and point O is the rotation center of workpiece. Set up the coordinate system X-Y with O as the origin. Point P is the grinding point, at which the workpiece and the grinding wheel are tangent. The angel between PO and X axis is $\mathsf{\tau }$. The radius of grinding wheel is ${\mathrm{R}}_{\mathrm{s}}$. Point N is an intersection point of the tangent line and X axis, and point M is an intersection point of the vertical line and X axis. When the angular speed of grinding wheel is 0, the linear speed at grinding point P is equal to the relative speed of point P on the workpiece to point P on the grinding wheel. To discuss the relation between the instantaneous linear speed of grinding point and rotation speed of workpiece, the following calculations are made.

(3)$\tan \mathsf{\gamma }=\raisebox{1ex}{$\mathrm{PM}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{MN}$}\right.$

Here, $\mathsf{\gamma }$ is the acute angel between the tangent line and X axis; PM is the vertical line, and MN is line from vertical point to intersection point N. From the triangular relation shown in Figure 2, the following Equation (4) can be easily obtained:

(4)$\{\begin{array}{c}\mathsf{\gamma }=\mathsf{\tau }+\mathsf{\delta }\\ \mathsf{\sigma }=\raisebox{1ex}{$\mathsf{\pi }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.-\mathsf{\delta }\end{array}$

Here, $\mathsf{\delta }$ is the angle between PO and PA. Note that V1 is the linear speed at grinding point P, V2 is the linear speed of point P on the workpiece, and V3 is the speed of grinding wheel carriage. $\mathsf{\sigma }$ is the angle between V1 and V2, and $\mathsf{\sigma }=\mathsf{\tau }+\raisebox{1ex}{$\mathsf{\pi }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.-\mathsf{\gamma }$.

The relation between V1, V2, and V3 satisfies the following Equation (5):

(5)$\{\begin{array}{c}\mathrm{V}2\left(\mathrm{t}\right)=\mathrm{V}\mathrm{1}\left(\mathrm{t}\right)·\cos \mathsf{\sigma }+\mathrm{V}3\left(\mathrm{t}\right)·\sin \mathsf{\tau }\\ \mathrm{V}1\left(\mathrm{t}\right)·\sin \mathsf{\sigma }-\mathrm{V}3\left(\mathrm{t}\right)·\cos \mathsf{\tau }=0\end{array}$

Then the instantaneous rotation speed of the workpiece is derived, which is shown in Equation (6),

(6)$\mathsf{\omega }\left(\mathrm{t}\right)=\frac{\mathrm{V}1\left(\mathrm{t}\right)·\left(\cos \mathsf{\sigma }+\sin \mathsf{\sigma }·\tan \mathsf{\tau }\right)}{\mathrm{OP}}$

Here, OP is the rotation diameter of workpiece.

With Equation (6), the function $\mathrm{h}(·)$ described the relationship between the rotation angle and the displacement of grinding wheel carriage can be deduced as following.

(7)$\{\begin{array}{c}\mathsf{\theta }\left(\mathrm{t}\right)={{\displaystyle \int }}_0^{\mathrm{t}}\mathsf{\omega }\left(\mathrm{t}\right)\mathrm{dt}\\ \mathrm{X}\left(\mathrm{t}\right)=\mathrm{h}\left(\mathsf{\theta }\right)\end{array}$

As mentioned above, the workpiece rotates with varying speed to keep the linear speed of grinding point constant.

2.2. Dimensional Error Analysis

For the indexable insert grinding process, it is complex and particular. Therefore, the error sources for dimensional precision decrease are much more complex, and many of them are coupled. Moreover, it is difficult to separate them from each other. For example, after long time grinding, the grinding wheel becomes blunt and the cutting reaction force increases, so the elastic deformation of feed system increases and the grinding wheel diameter decreases. Both of these increase the dimensional error. Dressing the grinding wheel further decreases the grinding wheel diameter. On the other hand, the feed system elastic deformation decreases. The influence of the grinding wheel diameter decrease on dimensional error is shown in Figure 3.

As it is shown in Figure 3, the solid circle is the real grinding wheel, and the dotted circle is the theoretical grinding wheel. $\mathrm{O}$ and ${\mathrm{O}}^{\prime }$ are the circle centers, respectively.$\mathrm{O}1$ is the rotation center of the workpiece. $\mathrm{Rs}$ and ${\mathrm{Rs}}^{\prime }$ represent the actual value and theoretical value of the grinding wheel radii, respectively. Point P is the theoretical tangent point. In fact, the workpiece is not tangent to the grinding wheel at point P. In the triangle ${\mathrm{O}}^{\prime }\mathrm{PO}$, with cosine law, Equation (8) can be derived.

(8)${\left(\mathrm{Rs}+\mathrm{e}\right)}^2={\left({\mathrm{Rs}}^{\prime }-\mathrm{Rs}\right)}^2+{{\mathrm{Rs}}^{\prime }}^2-2·\left({\mathrm{Rs}}^{\prime }-\mathrm{Rs}\right)·{\mathrm{Rs}}^{\prime }·\cos \mathsf{\vartheta }$

Here, $\mathrm{e}$ is the machining error. In this situation, the dimension of the workpiece will be large than the desired dimension.

Moreover, the spindle axis, ball screw, grinding wheel, and other parts generate much quantity of heat during grinding process, which result in thermal deformation of drive system. Finally, the dimensional precision decreases. In addition to, the static error of the machine tool, such as assembly error, also decreases the dimensional precision.

2.3. Problem Formulation

Consider the batch grinding system as a SISO nonlinear discrete-time system, which can be shown in Equation (9):

(9)${\mathrm{D}}_{\mathrm{k}}\left(\mathrm{t}+1\right)=\mathrm{f}({\mathrm{D}}_{\mathrm{k}}\left(\mathrm{t}\right),{\mathrm{u}}_{\mathrm{k}})$

where $\mathrm{t}=0,1,2,\dots ,\mathrm{N}$ is the sampling time index and N is the finite time interval of batch grinding system; $\mathrm{k}$ indicates the system repetition number and ${\mathrm{D}}_{\mathrm{k}}\left(\mathrm{t}\right)\in \mathrm{R}$ is the system output; ${\mathrm{u}}_{\mathrm{k}}\in \mathrm{R}$ denotes the system input, which is constant at all sampling time in the same iteration; and $\mathrm{f}(·)$ is an unknown function and continuously differentiable.

In the batch grinding system, only the terminal dimension ${\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)$ is measurable at the end of each run. Over iterations, the Equation (9) can be rewritten as a continuously differential function $\mathrm{g}(·)$ [25]. It is shown in Equation (10):

(10)$$\begin{gathered} {\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)=\mathrm{f}\left({\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}-1\right),{\mathrm{u}}_{\mathrm{k}}\right)=\mathrm{f}\left(\mathrm{f}\left({\mathrm{D}}_{\mathrm{k}}(\mathrm{N}-2\right),{\mathrm{u}}_{\mathrm{k}}\right),{\mathrm{u}}_{\mathrm{k}})=\mathrm{f}\left(\dots \mathrm{f}\left(\mathrm{f}\left({\mathrm{D}}_{\mathrm{k}}\left(0\right),{\mathrm{u}}_{\mathrm{k}}\right),{\mathrm{u}}_{\mathrm{k}}\right),\dots ,{\mathrm{u}}_{\mathrm{k}}\right) \\ =\mathrm{g}\left({\mathrm{D}}_{\mathrm{k}}\left(0\right),{\mathrm{u}}_{\mathrm{k}}\right) \end{gathered}$$

Here, ${\mathrm{D}}_{\mathrm{k}}\left(0\right)$ is the initial value of Equation (9). The system is controllable and the initial value ${\mathrm{D}}_{\mathrm{k}}\left(0\right)$ is identical for every iteration, ${\mathrm{D}}_{\mathrm{k}}\left(0\right)={\mathrm{D}}_{\mathrm{k}-1}\left(0\right)=\mathrm{c}$. $\mathrm{c}$ is a constant for the batch grinding system, which is the dimension of workpiece before grinding discussed in Section 4.

Assumption 1. 

The partial derivative of $f(\dots )$ with respect to control input $u_k$ is continuous.

Assumption 2. 

The Equation (9) satisfies the generalized Lipschitz condition. That is, $\left|\Delta D_k\left(N\right)\left|\le b\right|\Delta u_{k-1}\right|$ , here $\Delta D_k\left(N\right)=D_k\left(N\right)-D_{k-1}\left(N\right)$ , $\Delta u_{k-1}=u_{k-1}-u_{k-2}$ and b is positive.

Applying the mean value theorem into Equation (10), Equation (10) can be transformed to Equation (11),

(11)$\Delta {\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)=\frac{\partial \mathrm{g}}{\partial \mathrm{D}}\Delta {\mathrm{D}}_{\mathrm{k}}\left(0\right)+\frac{\partial \mathrm{g}}{\partial \mathrm{u}}\Delta {\mathrm{u}}_{\mathrm{k}}$

Since ${\mathrm{D}}_{\mathrm{k}}\left(0\right)={\mathrm{D}}_{\mathrm{k}-1}\left(0\right)$, then $\Delta {\mathrm{D}}_{\mathrm{k}}\left(0\right)$ is 0. Therefore, Equation (11) can be simplified as Equation (12),

(12)$\Delta {\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)=\frac{\partial \mathrm{g}}{\partial \mathrm{u}}\Delta {\mathrm{u}}_{\mathrm{k}}={\mathsf{\phi }}_{\mathrm{k}}\Delta {\mathrm{u}}_{\mathrm{k}}$

Here, ${\mathsf{\phi }}_{\mathrm{k}}=\frac{\partial \mathrm{g}}{\partial \mathrm{u}}$. According to Assumption 2, $\left|{\mathsf{\phi }}_{\mathrm{k}}\right|\le \mathrm{b}$. It is called pseudo-partial-derivative (PPD), which is a DL parameter. For more details, please see [23].

Then Equation (10) can be transformed to Equation (13),

(13)${\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)={\mathrm{D}}_{\mathrm{k}-1}\left(\mathrm{N}\right)+{\mathsf{\phi }}_{\mathrm{k}}\Delta {\mathrm{u}}_{\mathrm{k}}$

Remark 1. 

From a practical point, these assumptions imposed on the industrial system are reasonable and acceptable. Assumption 1 is a typical condition of control system design for general nonlinear systems. Assumption 2 imposes an upper bound limitation on the change rate of the system output driven by the changes in the control inputs. From the energy conservation point, the energy change rate in an industrial system is not infinite [24].

Remark 2. 

${\phi }_k$ is a differential signal and bounded to any k, which is unknown and iteratively changes.

An index function $\mathrm{J}={\left|\left|{\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)-{\mathrm{D}}_{\mathrm{k}-1}\left(\mathrm{N}\right)\right|\right|}^2+\mathsf{\mu }{\left|\left|{\mathrm{D}}_{\mathrm{k}}\left(\mathrm{N}\right)\right|\right|}^2$ [23] is introduced to estimate the value of ${\mathsf{\phi }}_{\mathrm{k}}$. By optimal condition, an estimation algorithm for ${\mathsf{\phi }}_{\mathrm{k}}$ is derived as Equation (14):

(14)${\hat{\mathsf{\phi }}}_{\mathrm{k}}={\hat{\mathsf{\phi }}}_{\mathrm{k}-1}+\frac{\mathsf{\rho }\left(\Delta {\mathrm{D}}_{\mathrm{k}-1}-{\hat{\mathsf{\phi }}}_{\mathrm{k}-1}\Delta {\mathrm{u}}_{\mathrm{k}-1}\right)\Delta {\mathrm{u}}_{\mathrm{k}-1}}{\mathsf{\mu }+{\left|\left|\Delta {\mathrm{u}}_{\mathrm{k}-1}\right|\right|}^2}$

Here it is bounded. The weighting factor $\mathsf{\mu }>0, 0<\rho <2$.

3. Sliding Mode Terminal Iterative Learning Compensation Method

• 1.. Define ${\mathrm{E}}_{\mathrm{k}}$ as the terminal dimensional error of the ${\mathrm{k}}_{\mathrm{th}}$ workpiece, which is shown as Equation (15):

  (15)${\mathrm{E}}_{\mathrm{k}}={\mathrm{D}}_{\mathrm{desired}}-{\mathrm{D}}_{\mathrm{k}}$

  Here,${ \mathrm{D}}_{\mathrm{k}}$ is the terminal dimension of the ${\mathrm{k}}_{\mathrm{th}}$ workpiece, which is measurable by laser sensor. ${\mathrm{D}}_{\mathrm{desired}}$ is the desired dimension of the workpiece. If ${\mathrm{E}}_{\mathrm{k}}$ is within the accepted error limit ${\mathrm{E}}_{\lim }$, the error compensation is ignored; otherwise, the compensation is operated.

  Substituting Equation (13) into Equation (15) yields Equation (16). Equation (16) is transformed as follows:

  (16)${\mathrm{E}}_{\mathrm{k}}={\mathrm{D}}_{\mathrm{desired}}-{\mathrm{D}}_{\mathrm{k}-1}-{\mathsf{\phi }}_{\mathrm{k}}\Delta {\mathrm{u}}_{\mathrm{k}}$

• 2.. Introduce an iteration-dependent slide surface function and an index function [26], which are two common functions in sliding mode control. The Equations are shown as follows:

  (17)${\mathrm{S}}_{\mathrm{k}}=\left({\mathrm{E}}_{\mathrm{k}}-{\mathrm{E}}_{\mathrm{k}-1}\right)+{\mathsf{\alpha }\mathrm{E}}_{\mathrm{k}}$

  (18)$\mathrm{J}={\left|\left|{\mathrm{S}}_{\mathrm{k}}-{\mathrm{S}}_{\mathrm{k}-1}\right|\right|}^2+\mathsf{\beta }{\left|\left|{\mathrm{S}}_{\mathrm{k}}\right|\right|}^2$

  where $\mathsf{\alpha }>0$, $\mathsf{\beta }>0$. They are weighting factors.

  Substituting Equation (16) and Equation (17) into Equation (18), the Equation (18) is transformed into Equation (19). The Equation (19) is shown as follows:

  (19)$$\begin{gathered} \mathrm{J}={\left|\left|-\left(1+\mathsf{\alpha }\right){\mathsf{\phi }}_{\mathrm{k}}\left({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{k}-1}\right)-{\mathrm{E}}_{\mathrm{k}-1}+{\mathrm{E}}_{\mathrm{k}-2}\right|\right|}^2 \\ +\mathsf{\beta }{\left|\left|{\mathsf{\alpha }\mathrm{E}}_{\mathrm{k}-1}-\left(1+\mathsf{\alpha }\right){\mathsf{\phi }}_{\mathrm{k}}\left({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{k}-1}\right)\right|\right|}^2 \end{gathered}$$

  In optimal conditions, the control input can be obtained in Equation (20), which is shown as follows:

  (20)${\mathrm{u}}_{\mathrm{k}}={\mathrm{u}}_{\mathrm{k}-1}+\frac{{\mathsf{\phi }}_{\mathrm{k}}\left({\mathrm{D}}_{\mathrm{desired}}-{\mathrm{D}}_{\mathrm{k}-1}\right)}{{\left|\left|{\mathsf{\phi }}_{\mathrm{k}}\right|\right|}^2}-\frac{\left(2+\mathsf{\alpha }+\mathsf{\beta }\right){\mathsf{\phi }}_{\mathrm{k}}{\mathrm{E}}_{\mathrm{k}-1}-{\mathsf{\phi }}_{\mathrm{k}}{\mathrm{E}}_{\mathrm{k}-2}}{\left(1+\mathsf{\alpha }\right)\left(1+\mathsf{\beta }\right){\left|\left|{\mathsf{\phi }}_{\mathrm{k}}\right|\right|}^2}$

• 3.. In the upper part, the compensation control input ${\mathrm{u}}_{\mathrm{k}}$ is deduced, which is the function of the previous two-dimensional errors, the previous dimension, and the estimation parameter. In this part, the convergence analysis of the grinding system with the control input is derived.

  In combination with Equation (20) and Equation (13), in virtue of Equation (16), then Equation (21) can be deduced. It is shown as follows [26]:

  (21)${\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}}=\frac{\left(2+\mathsf{\alpha }+\mathsf{\beta }\right){\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}-1}}{\left(1+\mathsf{\alpha }\right)\left(1+\mathsf{\beta }\right)}-\frac{{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}-2}}{\left(1+\mathsf{\alpha }\right)\left(1+\mathsf{\beta }\right)}$

  According to the above Equation (21), it is easy to derive the Equation (22),

  (22)$\left(1+\mathsf{\alpha }\right){\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}}-{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}-1}=\frac{\left(1+\mathsf{\alpha }\right){\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}-1}}{1+\mathsf{\beta }}-\frac{{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{E}}_{\mathrm{k}-2}}{1+\mathsf{\beta }}$

  Therefore, for the sliding surface Equation (17), it can be derived Equation (23),

  (23)${\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}}=\frac{{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}-1}}{1+\mathsf{\beta }}$

  Since $\mathsf{\beta }>0$, $\left|\left|{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}}\right|\right|<\left|\left|{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}-1}\right|\right|$.

  Then $\underset{\mathrm{k}\to \infty }{\lim }\left|\left|{\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}}\right|\right|$=0. Since ${\mathsf{\phi }}_{\mathrm{k}}$ is bounded, ${\mathsf{\phi }}_{\mathrm{k}}{}^{\mathrm{T}}$ is bounded. Therefore, the following Equation (24) can be derived.

  (24)$\underset{\mathrm{k}\to \infty }{\lim }\left|\left|{\mathrm{S}}_{\mathrm{k}}\right|\right|=0$

  Therefore, the sliding surface will converge to 0.

  Substituting Equation (17) into Equation (24), $\underset{\mathrm{k}\to \infty }{\lim }\left|\left|\left({\mathrm{E}}_{\mathrm{k}}-{\mathrm{E}}_{\mathrm{k}-1}\right)+{\mathsf{\alpha }\mathrm{E}}_{\mathrm{k}}\right|\right|$=0. Then, Equation (25) can be deduced.

  (25)$\underset{\mathrm{k}\to \infty }{\lim }\left|\left|{\mathrm{E}}_{\mathrm{k}}\right|\right|=\underset{\mathrm{k}\to \infty }{\lim }\left|\left|\frac{{\mathrm{E}}_{\mathrm{k}-1}}{1+\mathsf{\alpha }}\right|\right|=\underset{\mathrm{k}\to \infty }{\lim }\left|\left|\frac{{\mathrm{E}}_0}{{\left(1+\mathsf{\alpha }\right)}^{\mathrm{k}}}\right|\right|$

  Since $\mathsf{\alpha }>0$, the following Equation (26) can be derived.

  (26)$\underset{\mathrm{k}\to \infty }{\lim }\left|\left|{\mathrm{E}}_{\mathrm{k}}\right|\right|=0$

  Therefore, the terminal error will converge to 0.

  As mentioned above, the convergence analysis is made.

• 4.. As it is mentioned in Section 2, the system input is constant of every cycle at an arbitrary sampling time for one batch and it is executed by X axis. The compensation is executed by the compensation module if it is necessary; the compensation input for the ${\left(\mathrm{k}+1\right)}_{\mathrm{th}}$ workpiece is computed as follows.

  (27)${\mathrm{com}}_{\mathrm{k}+1=}{\mathrm{com}}_{\mathrm{k}}+\frac{{\mathsf{\phi }}_{\mathrm{k}+1}{\mathrm{E}}_{\mathrm{k}}}{{\left|\left|{\mathsf{\phi }}_{\mathrm{k}+1}\right|\right|}^2}-\frac{\left(2+\mathsf{\alpha }+\mathsf{\beta }\right){\mathsf{\phi }}_{\mathrm{k}+1}{\mathrm{E}}_{\mathrm{k}}-{\mathsf{\phi }}_{\mathrm{k}+1}{\mathrm{E}}_{\mathrm{k}-1}}{\left(1+\mathsf{\alpha }\right)\left(1+\mathsf{\beta }\right){\left|\left|{\mathsf{\phi }}_{\mathrm{k}+1}\right|\right|}^2}$

Here, the compensation is computed by the previous two-dimensional errors, and ${\mathsf{\phi }}_{\mathrm{k}}$ is computed in Equation (14).

The flowchart of the compensation method is shown in Figure 4. Note that in the judgment module, the compensation direction needs to be defined. If ${\mathrm{E}}_{\mathrm{k}}$ is less than $\left(-{\mathrm{E}}_{\lim }\right)$, that is to say, the dimension of this machined workpiece is less than the desired dimension, the compensation input will be positive. If ${\mathrm{E}}_{\mathrm{k}}$ is greater than positive error limit $\left(+{\mathrm{E}}_{\lim }\right)$, that is to say, the dimension of this machined workpiece is greater than the desired dimension, the compensation will be negative.

For the machine tool, stability and safe operation are of great importance. Under this premise, if the compensation value is very large, the compensation axis may result in irreversible fault. In order to ensure the safety of the machine tool, a compensation limit ${\mathrm{com}}_{\max }$ is introduced to avoid excessive fluctuation. If the computed ${\mathrm{com}}_{\mathrm{k}}$ is in an accepted range ($-{\mathrm{com}}_{\max }$, ${\mathrm{com}}_{\max }$), ${\mathrm{com}}_{\mathrm{k}}$ is executed; otherwise, the compensation value is set to ${\mathrm{com}}_{\max }$.

4. Simulation and Verification

4.1. Simulation

To verify the effectiveness of the upper error compensation algorithm, the step signal and impulse signal simulating dimensional errors is respectively introduced into the error compensation module. The desired value can be any value in the simulation. Here, it is defined as 9.435 mm, which is consistent with the desired dimension of experiment. Since the error acceptance of A level indexable insert is 0.005 mm, in order to increase dimensional precision, the error limit ${\mathrm{E}}_{\lim }$ is set to 0.003 mm for the simulation. Moreover, in order to operate steadily, the maximum compensation value ${\mathrm{C}}_{\max }$ is chosen as 0.015 mm. For the compensation module, the initial value of updated index ${\mathsf{\phi }}_1=1$. Values of weighting factors are defined according to the established conditions of Equations (14), (17), (18) and historical data analyzing, which have great influence on simulation results. Here, values of weighting factors are: $\mathsf{\rho }=1, \mathsf{\mu }=0.5, \mathsf{\beta }=10, \mathsf{\alpha }=10$.

The step error is 0.01 mm, which is imposed on the 15th and subsequent workpieces. The simulation result is shown in Figure 5. It is evident that the compensation method works effectively. The dimension recovers to the desired value only in two operation cycles. That is, the dimensions from the 16th workpiece are all converged to the desired value. For the second simulation, the desired dimension and the impulse error are the same to the first simulation. The impulse error is imposed on the 15th workpiece. The simulation result is shown in Figure 6. The dimension returns to the desired value in three operation cycles. The simulation results show that this algorithm can quickly converge to the desired value.

4.2. Experimental Investigation

A group of indexable inserts is grinded on Agathon 400combi shown in Figure 7. Figure 8 shows the geometry of this batch ${35}^{°}$ diamond indexable insert, and the material is YT15, a kind of tungsten–cobalt–titanium hard alloy. The diameter of the inscribed circle for the indexable insert is the controlled dimension (D in Figure 8), which is measured by a grating displacement sensor. The displacement sensor used in experiment is HEIDENHAIN MT60K. Its measurement range is 60 mm, and the system accuracy is up to $\pm 0.5 \mathrm{um}$. The measurement diagram is shown in Figure 9. After the indexable insert is grinded, the displacement sensor measured the diameter of inscribed circle and transmitted it to the controller to decide whether the compensation is needed or not. The grinding conditions for the experiments are shown in Table 1. For the experiment, first, the machine tool runs without grinding any workpiece for half an hour. Next, it begins to grind indexable inserts.

For this group experiment, the desired dimension of indexable insert is 9.435 mm. For the compensation module, the error limit ${\mathrm{E}}_{\lim }$ is set to 0.003 mm also in consideration of the error acceptance of A level. The maximum compensation value ${\mathrm{C}}_{\max }$ is chosen as 0.01 mm in order to ensure the safety. Furthermore, the initial value of updated index ${\mathsf{\phi }}_1=1$. For the experiment, under the grinding conditions shown in Table 1, a group of indexable inserts have been grinded on the machine tool without compensation. The dimension data are listed in columns 2 of Table 2, which are the baseline data.

Figure 10 shows the comparison of dimensional variation with different sets of compensation parameters, which is a group of compensation experiment. The first set compensation parameters are: $\mathsf{\rho }=1, \mathsf{\mu }=0.5, \mathsf{\beta }=1, \mathsf{\alpha }=1$. The second set compensation parameters are: $\mathsf{\rho }=1, \mathsf{\mu }=0.5, \mathsf{\beta }=10, \mathsf{\alpha }=10$. The third set compensation parameters are: $\mathsf{\rho }=1, \mathsf{\mu }=0.5, \mathsf{\beta }=10, \mathsf{\alpha }=20$. It is evident that under the first set compensation parameters, the dimensions fluctuated largely compared with the other two sets. The compensation effectiveness of the second set of compensation parameters is much better.

Figure 11 shows the comparison of dimensional variation before and after compensation. The blue line represents the baseline data without compensation. It is shown that the dimensional error is non monotonic. At the beginning of batch grinding, the dimensional error is relatively large. Then the dimensional error tends to within error tolerance. After a certain time grinding, the dimensional error increases largely again. Column 3 of Table 2 shows the data of the dimensions using this compensation method, which is indicated by the black line. The compensation parameters for this group experiment are: $\mathsf{\rho }=1, \mathsf{\mu }=0.5, \mathsf{\beta }=10, \mathsf{\alpha }=10$. It can be seen from Figure 11 that the dimensions fluctuate around the desired value after error compensation, and the compensation effectiveness is really good. From Figure 11 and Table 2, it is evident that there are only 33 workpieces out of 68 within the error tolerance (−0.005, 0.005) before error compensation, while after error compensation, the amount increases to 65. The qualified rate increases from 48.5% to 95.6% and the maximum error decreases from 0.019 to 0.007 mm. From the experiment results, it can be seen that this sliding mode terminal iterative learning compensation method performs well.

5. Conclusions

This paper focuses on the dimensional error compensation in CNC batch grinding. A data-driven sliding mode terminal iterative learning error compensation method is proposed to solve the problem. Through the above research, the following conclusions can be achieved:

• (1). Based on the theory of sliding mode terminal iterative learning control, the compensation method is presented. This method only needs to measure the terminal dimension using a touch probe, without using other sensors to detect all kinds of error sources. Moreover, there is no need of model information, which is data-driven;

• (2). Step and impulse input were introduced into the compensation module to simulate errors. Through many simulations, the values of compensation parameters were defined, and the simulation results showed that the compensation effectiveness was obvious;

• (3). In order to verify the compensation performance in an actual industrial environment, two groups of experiment were carried out on the machine tool. Based on the comparisons, the compensation parameters defined by simulation performed well and the performance effectiveness of this proposed compensation method was obvious. In the experiment, for 68 pieces of indexable inserts, the qualified rate increased from 48.5% to 95%.

Using this compensation method, the machined indexable inserts will remain within tolerance in successive machining. This method is appropriate for batch manufacturing, which is much more economical and universality.

Author Contributions

Conceptualization, T.C. and X.T.; methodology, T.C.; software, T.C.; validation, T.C.; formal analysis, T.C.; investigation, T.C. and X.T.; resources, X.T.; data curation, T.C.; writing—original draft preparation, T.C.; writing—review and editing, T.C.; visualization, T.C.; supervision, X.T.; project administration, X.T. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Structure of CNC batch grinding machine tool.

Figure 2. Diagram of grinding model.

Figure 3. Diagram of influence of grinding wheel diameter decrease on dimensional error.

Figure 4. Flowchart of dimensional error compensation method.

Figure 5. Step response of dimension compensation system.

Figure 6. Impulse response of dimension compensation system.

Figure 7. Agathon 400combi grinding machine tool.

Figure 8. Geometry diagram of the indexable insert.

Figure 9. Dimension measurement diagram. (a) The displacement sensor. (b) Actual measurement diagram.

Figure 10. Dimensional variation with three sets of compensation parameters for the compensation experiment.

Figure 11. Dimensional variation before and after error compensation.

Table 1

Grinding conditions for workpiece.

Parameters	Value
diamond wheel	resin bonded, mesh #220, concentration 100%
wheel diameter (mm)	350
wheel rotation (mm/s)	18
feed rate (mm/s)	0.03–0.05
no spark time(s)	3
coolant	water soluble
temperature of coolant (°C)	12–13

Table 2

Dimensions of the machined workpieces before and after compensation.

No.	Dimensions of the Group Workpieces without Compensation (mm)	Dimensions of the Group Workpieces with Compensation (mm)
1	9.422	9.422
2	9.424	9.434
3	9.428	9.438
4	9.428	9.434
5	9.43	9.436
6	9.43	9.436
7	9.428	9.434
8	9.428	9.434
9	9.428	9.434
10	9.428	9.434
11	9.432	9.438
12	9.432	9.434
13	9.436	9.438
14	9.434	9.432
15	9.432	9.431
16	9.432	9.436
17	9.428	9.432
18	9.428	9.436
19	9.43	9.438
20	9.434	9.442
21	9.436	9.436
22	9.434	9.434
23	9.438	9.438
24	9.434	9.434
25	9.434	9.434
26	9.434	9.434
27	9.434	9.434
28	9.434	9.434
29	9.434	9.434
30	9.434	9.434
31	9.434	9.434
32	9.434	9.434
33	9.438	9.438
34	9.436	9.436
35	9.436	9.436
36	9.436	9.436
37	9.436	9.436
38	9.436	9.436
39	9.436	9.436
40	9.434	9.434
41	9.436	9.436
42	9.438	9.438
43	9.44	9.439
44	9.446	9.440
45	9.448	9.436
46	9.446	9.434
47	9.446	9.434
48	9.448	9.436
49	9.45	9.438
50	9.452	9.436
51	9.448	9.432
52	9.45	9.434
53	9.454	9.438
54	9.452	9.432
55	9.45	9.431
56	9.452	9.438
57	9.452	9.438
58	9.45	9.436
59	9.45	9.436
60	9.446	9.432
61	9.446	9.436
62	9.444	9.434
63	9.442	9.432
64	9.444	9.438
65	9.444	9.438
66	9.444	9.438
67	9.444	9.438
68	9.44	9.434
Number of workpieces within the dimensional error tolerance 0.005 mm	33	65

References

1. Nojedeh, M.V.; Habibi, M.; Arezoo, B. Tool path accuracy enhancement through geometric error compensation. Int. J. Mach. Tools Manu.; 2011; 51, pp. 471-482. [DOI: https://dx.doi.org/10.1016/j.ijmachtools.2011.02.005]

2. Schmitz, T.L.; Ziegert, J.C.; Canning, J.S.; Raul, Z. Case study: A comparison of error sources in high-speed milling. Prec. Eng.; 2008; 32, pp. 126-133. [DOI: https://dx.doi.org/10.1016/j.precisioneng.2007.06.001]

3. Yin, S.H.; Ohmori, H.; Dai, Y.T.; Uehara, Y.; Chen, F.J.; Tang, H.N. ELID grinding characteristics of glass—Ceramic materials. Int. J. Mach. Tools Manu.; 2009; 49, pp. 333-338. [DOI: https://dx.doi.org/10.1016/j.ijmachtools.2008.10.013]

4. Khan, A.W.; Chen, W.Y. A methodology for systematic geometric error compensation in five-axis machine tools. Int. J. Adv. Manuf. Technol.; 2011; 53, pp. 615-628. [DOI: https://dx.doi.org/10.1007/s00170-010-2848-3]

5. Li, J.; Xie, F.G.; Liu, X.J.; Li, W.D.; Zhu, S.W. Geometric error identification and compensation of linear axes based on a novel 13-line method. Int. J. Adv. Manuf. Technol.; 2016; 87, pp. 2269-2283. [DOI: https://dx.doi.org/10.1007/s00170-016-8580-x]

6. Chen, J.; Lin, S.; Zhou, X.; Gu, T. A ballbar test for measurement and identification the comprehensive error of tilt table. Int. J. Mach. Tools Manu.; 2016; 103, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.ijmachtools.2015.12.002]

7. Wu, C.J.; Fan, J.W.; Wang, Q.H.; Pan, R.; Tang, Y.H.; Li, Z.S. Prediction and compensation of geometric error for translational axes in multi-axis machine tools. Int. J. Adv. Manuf. Technol.; 2018; 95, pp. 3413-3435. [DOI: https://dx.doi.org/10.1007/s00170-017-1385-8]

8. Li, Y.; Zhao, W.; Lan, S.; Ni, J.; Wu, W.; Lu, B. A review on spindle thermal error compensation in machine tools. Int. J. Mach. Tools Manu.; 2015; 95, pp. 25-38. [DOI: https://dx.doi.org/10.1016/j.ijmachtools.2015.04.008]

9. Wu, C.W.; Tang, C.H.; Chang, C.F.; Shao, Y.S. Thermal error compensation method for machine center. Int. J. Adv. Manuf. Technol.; 2012; 59, pp. 681-689. [DOI: https://dx.doi.org/10.1007/s00170-011-3533-x]

10. Wang, H.T.; Wang, L.P.; Li, T.M.; Han, J. Thermal sensor selection for the thermal error modeling of machine tool based on the fuzzy clustering method. Int. J. Adv. Manuf. Technol.; 2013; 69, pp. 121-126. [DOI: https://dx.doi.org/10.1007/s00170-013-4998-6]

11. Yang, J.; Shi, H.; Feng, B.; Zhao, L.; Ma, C.; Mei, X.S. Thermal error modeling and compensation for a high-speed motorized spindle. Int. J. Adv. Manuf. Technol.; 2015; 77, pp. 1005-1017. [DOI: https://dx.doi.org/10.1007/s00170-014-6535-7]

12. Liu, K.; Liu, Y.; Sun, M.J.; Wu, Y.L.; Zhu, T.J. Comprehensive thermal growth compensation method of spindle and servo axis error on a vertical drilling center. Int. J. Adv. Manuf. Technol.; 2017; 88, pp. 2507-2516. [DOI: https://dx.doi.org/10.1007/s00170-016-8972-y]

13. Zhang, Z.L.; Cai, L.G.; Cheng, Q.; Liu, Z.F.; Gu, P.H. A geometric error budget method to improve machining accuracy reliability of multi-axis machine tools. J. Intell. Manuf.; 2019; 30, pp. 495-519. [DOI: https://dx.doi.org/10.1007/s10845-016-1260-8]

14. He, G.; Ma, W.; Yu, G.; Lang, A. Modeling and experimental validation of cutting forces in five-axis ball-end milling based on true tooth trajectory. Int. J. Adv. Manuf. Technol.; 2015; 78, pp. 189-197. [DOI: https://dx.doi.org/10.1007/s00170-014-6548-2]

15. Ma, W.K.; He, G.Y.; Zhu, L.M.; Guo, L.Z. Tool deflection error compensation in five-axis ball-end milling of sculptured surface. Int. J. Adv. Manuf. Technol.; 2016; 84, pp. 1421-1430. [DOI: https://dx.doi.org/10.1007/s00170-015-7793-8]

16. Manna, M.A.; Ramesh, R.; Poo, A.N. Error compensation in machine tools-a review part I: Geometric, cutting-force induced and fixture-dependent errors. Int. J. Mach. Tools Manu.; 2007; 40, pp. 1235-1256.

17. Chen, T.T.; Tian, X.C.; Yan, L. Intelligent dimensional error pre-compensation in CNC grinding using iterative learning approach. Int. J. Adv. Manuf. Technol.; 2013; 67, pp. 1825-1835. [DOI: https://dx.doi.org/10.1007/s00170-012-4611-4]

18. Chen, T.T.; Tian, X.C. An intelligent self-learning method for dimensional error pre-compensation in CNC grinding. Int. J. Adv. Manuf. Technol.; 2014; 75, pp. 1349-1356. [DOI: https://dx.doi.org/10.1007/s00170-014-6249-x]

19. Xu, J.X.; Chen, Y.Q.; Heng, L.T.; Yamamoto, S. Terminal iterative learning control with an application to RTPCVD thickness control. Automatic; 1999; 35, pp. 1535-1542. [DOI: https://dx.doi.org/10.1016/S0005-1098(99)00076-X]

20. Chi, R.H.; Liang, H.; Lin, N.; Zhang, R.K.; Huang, B. Constraint data driven optimal terminal ILC of end product quality for a class of I/O constrained batch processes. Can. J. Chem. Eng.; 2018; 96, pp. 215-222. [DOI: https://dx.doi.org/10.1002/cjce.22934]

21. Knight, J.; Shirsavar, S.; Holderbaum, W. An improved reliability cuk based solar inverter with sliding mode control. IEEE Trans. Power Electron.; 2017; 21, pp. 1107-1115. [DOI: https://dx.doi.org/10.1109/TPEL.2006.876786]

22. Li, S.; Du, H.; Yu, X. Discrete-time terminal sliding mode control systems based on Euler’s discretization. IEEE Trans. Autom. Control; 2014; 59, pp. 546-552. [DOI: https://dx.doi.org/10.1109/TAC.2013.2273267]

23. Wang, W.H.; Hou, Z.S. New adaptive quasi-sliding mode control for nonlinear discrete-time systems. J. Syst. Eng. Electron.; 2008; 19, pp. 154-160.

24. Hou, Z.S.; Jin, S.T. A novel data-driven control approach for a class of discrete-time nonlinear systems. IEEE T. Contr. Syst. T.; 2011; 19, pp. 1549-1558. [DOI: https://dx.doi.org/10.1109/TCST.2010.2093136]

25. Chi, R.H.; Wang, D.W.; Hou, Z.S.; Jin, S.T. Data-driven optimal terminal iterative learning control. J. Process Control; 2012; 22, pp. 2026-2037. [DOI: https://dx.doi.org/10.1016/j.jprocont.2012.08.001]

26. Chi, R.H.; Wei, Y.C.; Wang, R.R. SMC-based data-driven terminal ILC for nonlinear discrete-time systems with unavailable intermediate measurements. Proceedings of the 31st Chinese Control and Decision Conference; Nanchang, China, 3–5 June 2019.