Implementation of Iterative Learning Control on a Pneumatic Actuator

Abstract

Pneumatic actuators demonstrate various nonlinear and uncertain behavior, and as a result, precise control of such actuators with model-based control schemes is challenging. The Iterative Learning Control (ILC) algorithm is a model-free control method usually used for repetitive processes. The ILC uses information from previous repetitions to learn about a system’s dynamics for generating a more suitable control signal. In this paper, an ILC method to overcome the nonlinearities and uncertainties in a pneumatic cylinder-piston actuator is suggested. The actuator is modeled using MATLAB SimScape blocks, and the ILC scheme has been expanded for controlling nonlinear, non-repetitive systems so that it can be used to control the considered pneumatic system. The simulation results show that the designed ILC controller is capable of tracking a non-repetitive reference signal and can overcome the internal and payload uncertainties with the precision of 0.002 m. Therefore, the ILC can be considered as an approach for controlling the pneumatic actuators, which is challenging to obtain their mathematical modeling.

1. Introduction

Pneumatic actuators have a wide range of applications in industrial automation and robotics, such as in load positioning [1], vehicles’ active suspension [2], air-brake systems [3], conveyor belt systems [4], pneumatic muscle actuators (PMA) [5,6] and soft robotics [7,8]. Nevertheless, the precise control of pneumatic systems requires further consideration as it cannot be satisfactorily achieved through many model-based control schemes. This emanates from the nonlinear and uncertain characteristics of pneumatic systems, which are caused by various factors such as dead-zone, air compressibility, frictional forces effects, and changes in airflow rate parameters, making obtaining an accurate model for such systems challenging [9,10,11,12].

Among the early attempts to control pneumatic systems, proportional-integral-derivative (PID) controllers [13,14] and their combination with acceleration feedback and nonlinear compensators [15,16] received more attention. Despite PID controllers’ popularity in industrial applications, they could not obtain the desired level of robustness, accuracy, and speed in controlling pneumatic systems. To improve robustness and overcome uncertainties, different robust control methods, including the H_∞ [17,18,19], quantitative feedback theory (QFT) [20,21] and sliding mode control (SMC) theory [22,23,24], have been proposed. The SMC controllers have demonstrated satisfactory performance in overcoming uncertainties in pneumatic systems, provided that the uncertainties’ boundaries are available to the controller. Adaptive controllers have also shown promising performance in controlling pneumatic actuators [25,26,27,28], and although they do not need any prior information on the uncertainties’ boundaries, an accurate model of each plant variation must be available to the controller for adaptation. However, obtaining accurate uncertainty boundaries or system models for most pneumatic actuators is challenging, and therefore, model-free control approaches, such as intelligent control methods, have attracted much interest. This includes fuzzy logic [6,29,30,31,32] and neural network-based controllers [33,34]. The intelligent control methods have also been combined with other controllers, such as the adaptive method, to achieve better performance in controlling pneumatic systems [35,36]. Nevertheless, employing a fuzzy logic controller requires expert knowledge or field test results, and using neural network controllers needs a training phase.

The iterative learning control (ILC) algorithm is an alternative intelligent method usually used for repetitive processes [37]. In this method, information from previous repetitions is used to learn about a system’s dynamics for generating a more suitable control signal. This learning process is performed in an iterative manner to improve the controller’s performance from one iteration to the other for achieving a zero-error convergence. ILC algorithms are particularly useful in real-time control systems, given their relatively quick response to the changes of the input signal. Many industrial processes are repetitive, which means the same control action should be performed repeatedly. Therefore, it is reasonable to make use of previously acquired data for improving a controller’s convergence and robustness in such processes. The difference between ILC and other learning-type control methods, such as the neural network and adaptive control, is that the ILC only modifies the control signal according to predefined control law. In contrast, other learning-type controllers monitor the system’s performance and update their control law during the process, accordingly [37]. Unlike an adaptive controller, ILC methods do not need any information on the system’s model and only operate based on the historical input and output. Moreover, contrasting to other intelligent controllers, no training is required, and a well-selected ILC method should be able to converge to the expected state within a few iterations [38]. The ILC has been used for controlling a pneumatic actuated X-Y table [39] and was combined with PID to control a simplified model of a pneumatic servo system [40]. Recently, the ILC method was proposed for accurate tracking of PMA [41,42] and controlling pneumatic valves [43].

This paper suggests an ILC method to overcome the nonlinearities and uncertainties resulting from air characteristics, pressure loss, leakage and load variations in a pneumatic system consisting of a valve, a mechanical actuator and the connecting pipes. The considered pneumatic system for this study is a cylinder-piston actuator, which is the most commonly used type of pneumatic system in the industry. Although ILC is a model-free control method, in Section 2, the mathematical model and physical properties of the considered pneumatic actuator have been discussed to demonstrate the system’s nonlinearities and uncertainties. MATLAB SimScape blocks have been used for simulating the system to provide a more detailed model than those used in previous studies on the ILC-controlled pneumatic systems. The asymptotic stability, monotonic convergent and zero steady-state error are the three aspects that should be considered in designing an ILC algorithm, which originally limits the ILC application to linear, repetitive systems. In this paper, the application of ILC has been further expanded to control the considered pneumatic system containing nonlinear behavior and responding to non-repetitive signals. The design of the proposed ILC method is done through a detailed mathematical analysis based on the physical properties of the considered actuator that is provided in Section 3. The simulation results, provided in Section 4, show that the designed ILC controller is capable of tracking a non-repetitive reference signal and can overcome the internal and payload uncertainties.

2. Mathematical Modeling and Simulation of the System

The considered pneumatic actuator in this paper follows the system given in [44], which is a double-acting pneumatic cylinder controlled by a 4-port-3-position (4/3) electro-pneumatic valve, as is depicted in Figure 1. According to the ideal gas law

$$\mathrm{P} = \mathrm{z}\mathsf{\rho }\mathrm{RT}$$

(1)

$\mathrm{P}, \mathsf{\rho }, \mathrm{z},$ and $\mathrm{T}$ are the gas’s pressure, density, compressibility factor, and temperature, respectively, and $\mathrm{R}$ is the gas constant. The ideal gas model is sufficiently accurate for modeling dry air under standard conditions.

In the control volume approach, each pneumatic component is considered as an internal node enclosed by a control surface. The mass flow rate in the control surface can be expressed as

$${\dot{\mathrm{m}}}_{\mathrm{in}}-{\dot{\mathrm{m}}}_{\mathrm{out}}=\frac{\partial \mathrm{M}}{\partial \mathrm{P}}\frac{{\mathrm{dP}}_{\mathrm{l}}}{\mathrm{dt}}+\frac{\partial \mathrm{M}}{\partial \mathrm{T}}\frac{{\mathrm{dT}}_{\mathrm{l}}}{\mathrm{dt}}$$

(2)

${\dot{\mathrm{m}}}_{\mathrm{in}}$ and ${\dot{\mathrm{m}}}_{\mathrm{out}}$ are the mass flow rates of gas entering and leaving the control surface. The gas volume properties are denoted by subscript $\mathrm{l}$, such as in ${\mathrm{P}}_{\mathrm{l}}$ and ${\mathrm{T}}_{\mathrm{l}}$, representing the pressure and temperature of the gas volume in the internal node, respectively. $\partial \mathrm{M}/\partial \mathrm{P}$ is the mass flow rate of the gas volume with respect to pressure at constant temperature and volume; and $\partial \mathrm{M}/\partial \mathrm{T}$ is the mass flow rate of the gas volume with respect to temperature at constant pressure and volume. For an ideal gas,

$$\begin{array}{c}\frac{\partial \mathrm{M}}{\partial \mathrm{P}}= \mathrm{V}\frac{{\mathsf{\rho }}_{\mathrm{l}}}{{\mathrm{P}}_{\mathrm{l}}} \\ \frac{\partial \mathrm{M}}{\partial \mathrm{T}}=-\mathrm{V}\frac{{\mathsf{\rho }}_{\mathrm{l}}}{{\mathrm{T}}_{\mathrm{l}}}\end{array}$$

(3)

Similarly, the control surface heat flow rate can be expressed as

$${\mathsf{\Phi }}_{\mathrm{in}}-{ \mathsf{\Phi }}_{\mathrm{out}}+ \mathrm{Q}=\frac{\partial \mathrm{U}}{\partial \mathrm{P}}\frac{{\mathrm{dP}}_{\mathrm{l}}}{\mathrm{dt}}+\frac{\partial \mathrm{U}}{\partial \mathrm{T}}\frac{{\mathrm{dT}}_{\mathrm{l}}}{\mathrm{dt}}$$

(4)

where ${\mathsf{\Phi }}_{\mathrm{in}}$ and ${\mathsf{\Phi }}_{\mathrm{out}}$ are the energy flow rates due to the gas entering and leaving the control surface. $\mathrm{Q}$ is the heat flow rate as a result of heat transferring between the control system and its surrounding, and $\mathrm{U}$ is the internal energy of the gas volume in the internal node. For an ideal gas,

$$\begin{array}{c}\frac{\partial \mathrm{U}}{\partial \mathrm{P}}= \mathrm{V}\left(\frac{{\mathrm{h}}_{\mathrm{l}}}{{\mathrm{zRT}}_{\mathrm{l}}}-1\right) \\ \frac{\partial \mathrm{U}}{\partial \mathrm{T}}={ \mathrm{V}\mathsf{\rho }}_{\mathrm{l}}\left({\mathrm{C}}_{{\mathrm{p}}_{\mathrm{l}}}-\frac{{\mathrm{h}}_{\mathrm{l}}}{{\mathrm{T}}_{\mathrm{l}}}\right)\end{array}$$

(5)

where ${\mathrm{h}}_{\mathrm{l}}$ and ${\mathrm{C}}_{{\mathrm{p}}_{\mathrm{l}}}$ are the specific enthalpy and specific heat capacity of the gas volume in the internal node, respectively.

For simulating the gas behavior in SimScape, a supply unit, as shown in Figure 2, has been set up. The gas properties, together with the thermal conductivity and dynamic viscosity, which are used in modeling the gas transport behavior, are defined in the SimScape Gas properties block. The system’s reference temperature and pressure, which are taken as atmospheric temperature and pressure in this model, are defined by the SimScape Reservoir block, and the air compressor is modeled by a Pressure source block. The SimScape Pressure source block has two ports and is capable of maintaining a constant pressure difference between its ports. Therefore, by connecting the input port to the reservoir, a set pressure gas can be injected into the system.

The pipes in the system can be modeled using SimScape’s Pipe (G), thermal mass, convective heat transfer and temperature source blocks, as shown in Figure 3.

The parameters used in calculating the gas pressure drop and convective heat transfer between the gas flowing in the pipe and the pipe’s wall are set in the Pipe (G) block. The pressure drop at each end of the pipe is given as

$$\begin{array}{c}{\mathrm{P}}_{\mathrm{in}}-{ \mathrm{P}}_{\mathrm{l}}={\left(\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathrm{A}}_{\mathrm{t}}}\right)}^2\times \left(\frac{1}{{\mathsf{\rho }}_{\mathrm{l}}}-\frac{1}{{\mathsf{\rho }}_{\mathrm{in}}}\right)+{ \Delta \mathrm{P}}_{\mathrm{inl}}\\ {\mathrm{P}}_{\mathrm{out}}-{ \mathrm{P}}_{\mathrm{l}}={\left(\frac{{\dot{\mathrm{m}}}_{\mathrm{out}}}{{\mathrm{A}}_{\mathrm{t}}}\right)}^2\times \left(\frac{1}{{\mathsf{\rho }}_{\mathrm{l}}}-\frac{1}{{\mathsf{\rho }}_{\mathrm{out}}}\right)+{ \Delta \mathrm{P}}_{\mathrm{outl}}\end{array}$$

(6)

where ${\mathrm{A}}_{\mathrm{t}}$ is the pipe’s cross-sectional area, and ${\Delta \mathrm{P}}_{\{·\}\mathrm{l}}$ is the pressure losses due to viscous friction, $\mathsf{\mu }$. The $\{·\}$ notation is used to abbreviate the equations and can be replaced by $“\mathrm{in}”$ or $“\mathrm{out}”$. This notation will be followed in the rest of this section. The pressure loss in a pipe depends on Reynolds numbers. The Reynolds number at each end of the pipe is equal to

$${\mathrm{Re}}_{\{·\}}=\frac{\left|{\dot{\mathrm{m}}}_{\{·\}}\right|\mathrm{D}}{{\mathrm{A}}_{\mathrm{t}}{\mathsf{\mu }}_{\mathrm{l}}}$$

(7)

where $\mathrm{D}$ denotes the pipe’s diameter. If the calculated Reynolds number is less than 2000, the gas flow follows the laminar regime, and the pressure drop is equal to

$${\Delta \mathrm{P}}_{\{·\}\mathrm{l}}={\mathrm{f}}_{\mathrm{shape}}\frac{{\dot{\mathrm{m}}}_{\{·\}}{ \mathsf{\mu }}_{\mathrm{l}}}{2{\mathsf{\rho }}_{\mathrm{l}}{\mathrm{D}}^2{\mathrm{A}}_{\mathrm{t}}}\times \frac{{\mathrm{L}}_{\mathrm{t}}+{\mathrm{L}}_{\mathrm{eqv}}}{2}$$

(8)

For Reynolds numbers greater than 4000, the gas flow follows the turbulent flow regime, with a pressure drop of

$${\Delta \mathrm{P}}_{\{·\}\mathrm{l}}={\mathrm{f}}_{{\mathrm{Darcy}}_{\{·\}}}\frac{{\dot{\mathrm{m}}}_{\{·\}}{ \mathsf{\mu }}_{\mathrm{l}}}{2{\mathsf{\rho }}_{\mathrm{l}}{\mathrm{DA}}_{\mathrm{t}}^2}\times \frac{{\mathrm{L}}_{\mathrm{t}}+{\mathrm{L}}_{\mathrm{eqv}}}{2}$$

(9)

${\mathrm{L}}_{\mathrm{t}}$ denotes the pipe’s length and ${\mathrm{L}}_{\mathrm{eqv}}$ is the aggregate equivalent length of local resistances, which has been taken as 10% of the pipe’s length. The shape factor, ${\mathrm{f}}_{\mathrm{shape}}$, is considered as a constant number equal to 2.59, but the Darcy friction factor, ${\mathrm{f}}_{\mathrm{Darcy}}$, is a function of the Reynolds number and is equal to

$${\mathrm{f}}_{{\mathrm{Darcy}}_{\{·\}}}={\left[-1.8\log \left(\frac{6.9}{{\mathrm{Re}}_{\{·\}}}+{\left(\frac{{\mathsf{\epsilon }}_{\mathrm{rough}}}{3.7\mathrm{D}}\right)}^{1.11}\right)\right]}^{-2}$$

(10)

where ${\mathsf{\epsilon }}_{\mathrm{rough}}$ is the internal surface absolute roughness. For the Reynolds numbers between 2000 and 4000, a transition between laminar and turbulence regimes has been assumed.

The convective heat transfer between the gas flowing in the pipe and the pipe’s wall are modeled as

$$\begin{array}{cc}\hfill \mathrm{Q} =\left|\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}-{\dot{\mathrm{m}}}_{\mathrm{out}}}{2}\right|& {{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{avg}}\left({\mathrm{T}}_{\mathrm{H}}-{\mathrm{T}}_{\mathrm{in}}\right){\displaystyle (}1- \exp {\displaystyle (}-\frac{{\mathrm{N}}_{\mathrm{u}}{\mathsf{\kappa }}_{\mathrm{avg}}{\mathsf{\pi }\mathrm{L}}_{\mathrm{t}}}{\left|\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}-{\dot{\mathrm{m}}}_{\mathrm{out}}}{2}\right|{{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{avg}}}{\displaystyle )}{\displaystyle )}\hfill \\ & +{ \mathsf{\kappa }}_{\mathrm{l}}{\mathsf{\pi }\mathrm{L}}_{\mathrm{t}}\left({\mathrm{T}}_{\mathrm{H}}-{ \mathrm{T}}_{\mathrm{l}}\right)\hfill \end{array}$$

(11)

where ${{\mathrm{C}}_{\mathrm{p}}}_{\mathrm{avg}}$ and ${\mathsf{\kappa }}_{\mathrm{avg}}$ are the specific heat capacity and thermal conductivity calculated at the average temperature. ${\mathrm{N}}_{\mathrm{u}}$ is the Nusselt number, ${\mathsf{\kappa }}_{\mathrm{l}}$ is the gas volume thermal conductivity, and ${\mathrm{T}}_{\mathrm{in}}$ and ${\mathrm{T}}_{\mathrm{H}}$ are the inlet and pipe internal wall temperatures, respectively. The value of the Nusselt number depends on the flow regime. For the laminar flow, the Nusselt number is taken as a constant, 3.66, and for turbulence flow is equal to

$${\mathrm{N}}_{\mathrm{u}\_\mathrm{tur}}=\frac{\frac{{\mathrm{f}}_{\mathrm{Darcy}}}{8} \left({\mathrm{Re}}_{\mathrm{avg}}-1000\right){ \mathrm{Pr}}_{\mathrm{avg}}}{1+12.7\sqrt{\frac{{\mathrm{f}}_{\mathrm{Darcy}}}{8}}+{ \mathrm{Pr}}_{\mathrm{avg}}^{2/3}-1}$$

(12)

${\mathrm{Re}}_{\mathrm{avg}}$ and ${\mathrm{Pr}}_{\mathrm{avg}}$ are the Reynolds and Prandtl numbers, which are obtained at the average temperature.

The pipe properties for modeling the conduction heat transfer in the pipe’s wall can be adjusted in the thermal mass block, and the parameters for modeling the convective heat transfer between the pipe’s wall and the surrounding environment can be set in the convective heat transfer and temperature source blocks. The temperature source can maintain the atmospheric temperature regardless of the amount of heat flow into the system. The conduction heat transfer in the pipe’s wall can be mathematically modeled by

$$\mathrm{Q} ={\mathsf{\rho }}_{\mathrm{t}}{\mathsf{\pi }\mathrm{DL}}_{\mathrm{t}}{\mathrm{C}}_{\mathrm{t}}\left({\mathrm{T}}_{\mathrm{t}}-{ \mathrm{T}}_{\mathrm{H}}\right)$$

(13)

where ${\mathsf{\rho }}_{\mathrm{t}}$, ${\mathrm{C}}_{\mathrm{t}},$ and ${\mathrm{T}}_{\mathrm{t}}$ are the pipe’s wall density, specific heat capacity, and temperature. Finally, the convective heat transfer between the pipe’s wall and the surrounding environment can be expressed as

$$\mathrm{Q} ={\mathsf{\kappa }}_{\mathrm{air}}{\mathsf{\pi }\mathrm{DL}}_{\mathrm{t}}\left({\mathrm{T}}_{\mathrm{atm}}-{ \mathrm{T}}_{\mathrm{t}}\right)$$

(14)

where ${\mathsf{\kappa }}_{\mathrm{air}}$ and ${\mathrm{T}}_{\mathrm{atm}}$ are the surrounding air thermal conductivity and temperature.

The next component of the pneumatic actuator is the valve, which can be modeled as a set of restrictions capable of controlling the gas flow according to an input control signal. The restriction causes contraction of the gas at its input port, followed by the gas expansion at its output port, which results in a pressure drop across the ports. Considering the process as adiabatic, this pressure difference can be modeled as

$$\Delta \mathrm{P} =\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathrm{C}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{R}}}.\left|\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathsf{\rho }}_{\mathrm{R}}{\mathrm{C}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{R}}}\right|\left(\frac{1+ \mathrm{r}}{2}\left(1- \mathrm{r}\frac{{\mathsf{\rho }}_{\mathrm{R}}}{{\mathsf{\rho }}_{\mathrm{in}}}\right)- \mathrm{r}\left(1- \mathrm{r}\frac{{\mathsf{\rho }}_{\mathrm{R}}}{{\mathsf{\rho }}_{\mathrm{out}}}\right)\right)$$

(15)

where ${\mathrm{A}}_{\mathrm{R}}$ is the restriction’s cross-sectional area, and ${\mathsf{\rho }}_{\mathrm{R}}$ is the gas volume density at it. ${\mathrm{C}}_{\mathrm{d}}$ denotes the discharge coefficient, and $\mathrm{r}={ \mathrm{A}}_{\mathrm{R}}/{\mathrm{A}}_{\mathrm{port}}$ is the ratio of the restriction cross-sectional area to the ports cross-sectional area. It is assumed that both ports have the same dimensions. As is seen from Equation (15), the pressure difference is proportional to the square of the gas flow rate, ${\dot{\mathrm{m}}}_{\mathrm{in}}$, which is typical in the turbulence regime. However, in the laminar regime, the pressure difference is linearly proportional to the gas flow rate, and thus $\Delta \mathrm{P}$ can be approximated as

$$\Delta {\mathrm{p}}_{\mathrm{lam}}=\sqrt{\frac{{\mathsf{\rho }}_{\mathrm{R}}.\frac{{\mathrm{p}}_{\mathrm{in}}+{ \mathrm{p}}_{\mathrm{out}}}{2}\left(1-{ \mathrm{B}}_{\mathrm{lam}}\right)}{2}}\left(1- \mathrm{r}\right)$$

(16)

where ${\mathrm{p}}_{\mathrm{in}}$ and ${\mathrm{p}}_{\mathrm{out}}$ are the pressures of the inflow and outflow gases into the restriction, and ${\mathrm{B}}_{\mathrm{lam}}$ is the laminar flow pressure ratio, which is taken as a constant, 0.999. The amount of pressure at the restriction is equal to

$${\mathrm{P}}_{\mathrm{R}}={\mathrm{p}}_{\mathrm{in}}-\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathrm{C}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{R}}}.\left|\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathsf{\rho }}_{\mathrm{R}}{\mathrm{C}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{R}}}\right|\left(\frac{1+ \mathrm{r}}{2}\left(1- \mathrm{r}\frac{{\mathsf{\rho }}_{\mathrm{R}}}{{\mathsf{\rho }}_{\mathrm{in}}}\right)\right)$$

(17)

For the laminar regime, this can be approximated as

$${\mathrm{P}}_{\mathrm{R}}=\frac{{\mathrm{p}}_{\mathrm{in}}+{ \mathrm{p}}_{\mathrm{out}}}{2}-\frac{1}{{\mathsf{\rho }}_{\mathrm{R}}}{\left(\frac{{\dot{\mathrm{m}}}_{\mathrm{in}}}{{\mathrm{C}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{R}}}\right)}^2\frac{1-{ \mathrm{r}}^2}{2}$$

(18)

The local restriction (G) block has been used for modeling the valve in the SimScape environment. Each active state of the valve can be modeled by two local restriction blocks connecting ports P and T to ports A and B in Figure 1. The block allows modeling the valve leakage by defining a non-zero minimum restriction area, ${\mathrm{A}}_{\mathrm{leakage}}$. Moreover, the valve spool displacement, $\mathrm{x}$, has been used to adjust the orifice area. The orifice cross-sectional area is linearly proportional to the spool displacement as

$${\mathrm{A}}_{\mathrm{R}}=\frac{{\mathrm{A}}_{\max }-{\mathrm{A}}_{\mathrm{leakage}}}{{\mathrm{x}}_{\max }-{\mathrm{x}}_{\mathrm{leakage}}}\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{leakage}}\right)+{\mathrm{A}}_{\mathrm{leakage}}$$

(19)

where ${\mathrm{x}}_{\max }$ is the maximum spool displacement causing the maximum cross-sectional area in the restriction, ${\mathrm{A}}_{\max }$, and ${\mathrm{x}}_{\mathrm{leakage}}$ is the minimum value for the spool displacement. Figure 4 shows the model of the valve using SimScape components.

The complete model for the cylinder using SimScape blocks is depicted in Figure 5. The cylinder consists of two chambers, where each can be modeled using the translational mechanical converter (G) block in the SimScape to model the relation between the gas pressure inside a chamber and the applied mechanical force to the interface. Each cylinder chamber can be considered as an internal node, with a mass flow rate of

$$\dot{\mathrm{m}}=\frac{\partial \mathrm{M}}{\partial \mathrm{p}}.\frac{{\mathrm{dp}}_{\mathrm{l}}}{\mathrm{dt}}+\frac{\partial \mathrm{M}}{\partial \mathrm{T}}.\frac{{\mathrm{dT}}_{\mathrm{l}}}{\mathrm{dt}}+{ \mathsf{\rho }}_{\mathrm{l}}\frac{\mathrm{dV}}{\mathrm{dt}}$$

(20)

A heat flow rate of

$$\varnothing + \mathrm{Q} =\frac{\partial \mathrm{U}}{\partial \mathrm{p}}.\frac{{\mathrm{dp}}_{\mathrm{l}}}{\mathrm{dt}}+\frac{\partial \mathrm{U}}{\partial \mathrm{T}}.\frac{{\mathrm{dT}}_{\mathrm{l}}}{\mathrm{dt}}+{ \mathsf{\rho }}_{\mathrm{l}}{\mathrm{h}}_{\mathrm{l}}\frac{\mathrm{dV}}{\mathrm{dt}}$$

(21)

where $\varnothing$ is the energy flow rate as a result of the gas transportation into/out of the chamber, and $\mathrm{Q}$ is due to the convective heat transfer between the gas in the chamber and the cylinder’s body. The partial derivative terms in Equations (20) and (21) can be obtained from Equations (3) and (5), respectively. Equations (11)–(14) are also applicable here to model the heat transfer between the gas in the chamber and the surrounding environment.

The gas volume in the chamber depends on the displacement of the moving interface, and is equal to

$$\mathrm{V} ={ \mathrm{V}}_{\mathrm{d}} \pm { \mathrm{A}}_{\mathrm{p}}{\mathrm{x}}_{\mathrm{c}}$$

(22)

where ${\mathrm{V}}_{\mathrm{d}}$ denotes the dead volume, and ${\mathrm{A}}_{\mathrm{p}}$ and ${\mathrm{x}}_{\mathrm{c}}$ represent the interface cross-sectional area and displacement. The sign of the displacement value depends on the movement direction of the interface. The applied force to the interface can be expressed as

$$\mathrm{F} =\left({\mathrm{p}}_{\mathrm{atm}}-{ \mathrm{p}}_{\mathrm{l}}\right){\mathrm{A}}_{\mathrm{p}}$$

(23)

where ${\mathrm{p}}_{\mathrm{atm}}$ is the atmospheric pressure. Therefore, the total force caused by the gas pressure in chambers 1 and 2, considering the forces’ directions, is equal to

$${\mathrm{F}}_{\mathrm{p}}=\left({\mathrm{p}}_{\mathrm{l}1}-{ \mathrm{p}}_{\mathrm{l}2}\right){\mathrm{A}}_{\mathrm{p}}$$

(24)

Moreover, a translational hard stop block is used to restrict the cylinder interface movement within the length of the cylinder. The viscous friction coefficient for the piston movement can be modelled by a translational damper block. Two mass blocks are used for modelling the load and piston masses. The convective heat transfer, thermal mass, and temperature source blocks can model the heat transfer between the chamber gases and the surrounding environment. The motion equation for the load connected to the piston rod is equal to

$$\left({\mathrm{M}}_{\mathrm{L}}+{ \mathrm{M}}_{\mathrm{P}}\right){\dot{\mathrm{x}}}_{\mathrm{c}}={ \mathrm{F}}_{\mathrm{p}}- \mathsf{\beta }\frac{{\mathrm{dx}}_{\mathrm{c}}}{\mathrm{dt}}$$

(25)

where ${\mathrm{x}}_{\mathrm{c}}$ denotes the rod’s displacement, ${\mathrm{M}}_{\mathrm{P}}$ is the piston mass, and ${\mathrm{M}}_{\mathrm{L}}$ is the load mass, which causes the payload force of ${\mathrm{F}}_{\mathrm{L}}={ \mathrm{M}}_{\mathrm{L}}{\dot{\mathrm{x}}}_{\mathrm{c}}$ on the actuator. $\mathsf{\beta }$ represents the viscous friction coefficient for the piston movement.

3. Controller Design

This section covers the procedure for designing an ILC controller for the considered pneumatic actuator covered in Section 2. The ILC method can generally be applied to repetitive processes. In such processes, every iteration starts and ends at the same conditions and lasts for a fixed period of time, $\mathrm{T}$. The state-space equation for a repetitive linear time-invariant system (with no feedthrough) can be represented as

$$\begin{array}{c}{\dot{x}}_{\mathrm{k}}\left(\mathrm{t}\right)=A_c{x}_{\mathrm{k}}\left(\mathrm{t}\right)+B_c{\mathrm{u}}_{\mathrm{k}}\left(\mathrm{t}\right)\\ {\mathrm{y}}_{\mathrm{k}}\left(\mathrm{t}\right)=C_c{x}_{\mathrm{k}}\left(\mathrm{t}\right) \end{array}; \mathrm{t}\in \left[0,\mathrm{T}\right] x_{\mathrm{k}}\left(0\right)=x_0 \forall \mathrm{k}$$

(26)

$\mathrm{k} >0$ denotes the trial number, and $x_{\mathrm{k}}\left(\mathrm{t}\right)\in {ℝ}^{\mathrm{n}}$, ${\mathrm{y}}_{\mathrm{k}}\left(\mathrm{t}\right)\in ℝ,$ and ${\mathrm{u}}_{\mathrm{k}}\left(\mathrm{t}\right)\in ℝ$ respectively represent the state variable vector, output, and input of the system at the ${\mathrm{k}}^{\mathrm{th}}$ iteration. $A_c\in {ℝ}^{\mathrm{n}\times \mathrm{n}},B_c\in {ℝ}^{\mathrm{n}\times 1}, C_c\in {ℝ}^{1\times \mathrm{n}}$ are the state space system matrices. Suppose the system’s desired output for every iteration is ${\mathrm{y}}_{\mathrm{d}}\left(\mathrm{t}\right)$, which makes the error at the ${\mathrm{k}}^{\mathrm{th}}$ iteration equal to

$${\mathrm{e}}_{\mathrm{k}}\left(\mathrm{t}\right)={ \mathrm{y}}_{\mathrm{d}}\left(\mathrm{t}\right)-{ \mathrm{y}}_{\mathrm{k}}\left(\mathrm{t}\right)$$

(27)

The idea of the ILC is to define a control law using previous trials’ information such that the error monotonically decreases in every new iteration until it reaches zero ($\underset{\mathrm{k}\to \infty }{\lim }{\mathrm{e}}_{\mathrm{k}}\left(\mathrm{t}\right)=0$). The control signal is calculated according to a recursive law as

$${\mathrm{u}}_{\mathrm{k}+1}\left(\mathrm{t}\right)=\mathrm{F}\left({\mathrm{u}}_0\left({\mathrm{t}}^{\prime }\right),\dots ,{\mathrm{u}}_{\mathrm{k}}\left({\mathrm{t}}^{\prime }\right),{\mathrm{y}}_0\left({\mathrm{t}}^{\prime }\right),\dots ,{\mathrm{y}}_{\mathrm{k}}\left({\mathrm{t}}^{\prime }\right),{\mathrm{y}}_{\mathrm{d}}\left({\mathrm{t}}^{\prime }\right)\right); 0 \le {\mathrm{t}}^{\prime }\le \mathrm{T}$$

(28)

If $\mathrm{F}$ is designed in a way that ${\mathrm{t}}^{\prime }> \mathrm{t}$, then the learning law is known as noncausal. Moreover, the control signal generated at the ${\left(\mathrm{k}+1\right)}^{\mathrm{th}}$ iteration can be calculated based on the information collected from all previous iterations. This method is called the high-order ILC (HOILC). However, a simplified law is preferred as long as it can attain the convergence with satisfactory speed.

The ILC algorithm can be presented in discrete-time, which is more suitable for being implemented by a microcontroller. For a sampling time of ${\mathrm{T}}_{\mathrm{s}}$, where $\mathrm{T} ={ \mathrm{NT}}_{\mathrm{s}}$, the system can be considered as

$$\begin{array}{c}{x}_{\mathrm{k}}\left[\mathrm{i}+1\right]=A{x}_{\mathrm{k}}\left[\mathrm{i}\right]+B{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]\\ {\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]=C{x}_{\mathrm{k}}\left[\mathrm{i}\right] \end{array} ; \mathrm{i} \in \left[0,\mathrm{N}\right] x_{\mathrm{k}}\left[0\right]=x_0 \forall \mathrm{k}$$

(29)

$A={ \mathrm{e}}^{A_c{\mathrm{T}}_{\mathrm{s}}}$, $B={{\displaystyle \int }}_0^{{\mathrm{T}}_{\mathrm{s}}}{\mathrm{e}}^{A_c\mathsf{\alpha }}\mathrm{d}\mathsf{\alpha }{B}_c$ and $C=C_c$. The system’s output can then be calculated as

$${\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]=C{\left(\mathrm{q}I-A\right)}^{-1} B{ \mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+C{A}^{\mathrm{k}} x_0= \mathrm{P}\left(\mathrm{q}\right){\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+{ \mathrm{d}}_{\mathrm{k}}$$

(30)

where $\mathrm{q}$ denotes the forward time-shift operator as $\mathrm{qx}\left[\mathrm{i}\right]\equiv \mathrm{x}\left[\mathrm{i}+1\right]$. For a rational LTI system, $\mathrm{P}\left(\mathrm{q}\right)$ can be expanded into

$$\mathrm{P}\left(\mathrm{q}\right)=CB{\mathrm{q}}^{-1}+CAB{\mathrm{q}}^{-2}+C{A}^{2}B{\mathrm{q}}^{-3}+\dots$$

(31)

The error at the ${\mathrm{k}}^{\mathrm{th}}$ iteration is equal to

$${\mathrm{e}}_{\mathrm{k}}\left[\mathrm{i}\right]={\mathrm{y}}_{\mathrm{d}}\left[\mathrm{i}\right]-{\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]$$

(32)

The ILC control law in the discrete-time domain can then be presented as

$${\mathrm{u}}_{\mathrm{k}+1}\left[\mathrm{i}\right]=\mathrm{F}\left({\mathrm{u}}_{\mathrm{k}}[·],\dots ,{ \mathrm{u}}_0[·],{ \mathrm{y}}_{\mathrm{k}}[·],\dots ,{ \mathrm{y}}_0[·],{\mathrm{y}}_{\mathrm{d}}[·]\right)$$

(33)

where $[·]$ represents any sample in the range of $\left[0, \mathrm{N}\right].$ In such cases, it is common to implement the ILC method in the form of digital filters.

We begin our design by considering a general form for the ILC law as

$${\mathrm{u}}_{\mathrm{k}+1}\left[\mathrm{i}\right]= \mathrm{Q}\left(\mathrm{q}\right){\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+ \mathrm{L}\left(\mathrm{q}\right){\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]+ \mathrm{M}\left(\mathrm{q}\right){\mathrm{y}}_{\mathrm{d}}\left[\mathrm{i}\right]$$

(34)

To represent (34) in matrix format, let us consider $\mathrm{N}$ sample sequences for the input, output and desired signals as

$$\begin{array}{c}{u}_{\mathrm{k}}={{\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{u}}_{\mathrm{k}}\left[1\right]& {\mathrm{u}}_{\mathrm{k}}\left[2\right] \end{array}& \begin{array}{cc}\cdots & {\mathrm{u}}_{\mathrm{k}}\left[\mathrm{N}\right]\end{array}\end{array}\right]}^{\mathrm{T}}}^{\mathrm{T}}\\ y_{\mathrm{k}}={\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{y}}_{\mathrm{k}}\left[1\right]& {\mathrm{y}}_{\mathrm{k}}\left[2\right] \end{array}& \begin{array}{cc}\cdots & {\mathrm{y}}_{\mathrm{k}}\left[\mathrm{N}\right]\end{array}\end{array}\right]}^{\mathrm{T}}\\ y_{\mathrm{d}}={\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{y}}_{\mathrm{d}}\left[1\right]& {\mathrm{y}}_{\mathrm{d}}\left[2\right]\end{array}& \begin{array}{cc}\cdots & { \mathrm{y}}_{\mathrm{d}}\left[\mathrm{N}\right]\end{array}\end{array}\right]}^{\mathrm{T}}\\ d_{\mathrm{k}}={\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{d}}_{\mathrm{k}} & {\mathrm{d}}_{\mathrm{k}} \end{array}& \begin{array}{cc}\cdots & {\mathrm{d}}_{\mathrm{k}} \end{array}\end{array}\right]}^{\mathrm{T}}\end{array}$$

(35)

And $M$, $L$, $Q,$ and $P$ as matrices equal to

$$\begin{array}{cc}M=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\mathrm{m}}_0\\ {\mathrm{m}}_1\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{m}}_{\mathrm{N}-1}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{m}}_{-1}\\ {\mathrm{m}}_0\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{m}}_{\mathrm{N}-2}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{m}}_{-\left(\mathrm{N}-1\right)}\\ {\mathrm{m}}_{-\left(\mathrm{N}-2\right)}\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{m}}_0\end{array}\end{array}\end{array}\end{array}\right]& L=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\mathrm{l}}_0\\ {\mathrm{l}}_1\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{l}}_{\mathrm{N}-1}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{l}}_{-1}\\ {\mathrm{l}}_0\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{l}}_{\mathrm{N}-2}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{l}}_{-\left(\mathrm{N}-1\right)}\\ {\mathrm{l}}_{-\left(\mathrm{N}-2\right)}\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{l}}_0\end{array}\end{array}\end{array}\end{array}\right]\\ Q=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\mathrm{q}}_0\\ {\mathrm{q}}_1\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{q}}_{\mathrm{N}-1}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{q}}_{-1}\\ {\mathrm{q}}_0\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{q}}_{\mathrm{N}-2}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathrm{q}}_{-\left(\mathrm{N}-1\right)}\\ {\mathrm{q}}_{-\left(\mathrm{N}-2\right)}\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{q}}_0\end{array}\end{array}\end{array}\end{array}\right]& P=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\mathrm{P}}_1\\ {\mathrm{P}}_2\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{P}}_{\mathrm{N}}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ {\mathrm{P}}_1\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{P}}_{\mathrm{N}-1}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}⋮\\ {\mathrm{P}}_1\end{array}\end{array}\end{array}\end{array}\right]\end{array}$$

(36)

Therefore, (34) can be represented as

$$\begin{array}{cc}\hfill u_{\mathrm{k}+1}& =Q{u}_{\mathrm{k}}+L{y}_{\mathrm{k}}+M{y}_{\mathrm{d}}\hfill \\ & =\left(Q+LP\right)u_{\mathrm{k}}+\left(M{y}_{\mathrm{d}}+L{d}_{\mathrm{k}}\right)\hfill \end{array}$$

(37)

An ILC method is regarded to be asymptotically stable (AS) if

$$\exists \overline{\mathrm{u}}\in ℝ : \left|{\mathrm{u}}_{\mathrm{k}+1}\left[\mathrm{i}\right]\right| \le \overline{\mathrm{u}} \forall \mathrm{k} >0,\forall \mathrm{i}\in \left[1,\mathrm{N} -1\right]$$

(38)

The converged control signal can be defined as ${\mathrm{u}}_{\infty }\left[\mathrm{i}\right]=\underset{\mathrm{k}\to \infty }{\lim }{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]$, and in order for (37) to be AS,

$$\mathsf{\rho }\left(Q+LP\right)<1$$

(39)

where $\mathsf{\rho }\left(A\right)=\underset{\mathrm{j}}{\max }\left|{\mathsf{\lambda }}_{\mathrm{j}}\left(A\right)\right|$ and ${\mathsf{\lambda }}_{\mathrm{j}}\left(A\right)$ is the ${\mathrm{j}}^{\mathrm{th}}$ eigenvalue of matrix $A$.

Theorem 1.

The asymptotic error of the controlled system is equal to

$$\begin{array}{c}\\ e_{\infty }=\left(I-P{\left(I-Q-LP\right)}^{-1}M\right)y_{\mathrm{d}}-\left(P{\left(I-Q-LP\right)}^{-1}L+I\right)d_{\infty }\end{array}$$

(40)

where$\mathbf{I}$is an$N\times N$identity matrix.

Proof of Theorem 1.

Based on (32), $e_{\mathrm{k}}=y_{\mathrm{d}}-y_{\mathrm{k}}$ and following (30) $y_{\mathrm{k}}=P{u}_{\mathrm{k}}+d_{\mathrm{k}}$, which makes $e_{\mathrm{k}}=y_{\mathrm{d}}-P{u}_{\mathrm{k}}-d_{\mathrm{k}}$. Therefore, the asymptotic error, $k\to \infty ,$ can be obtained as $e_{\infty }=y_{\mathrm{d}}-P{u}_{\infty }-d_{\infty }$ and by using (37) when $k\to \infty$

$$\begin{array}{c} u_{\infty }=\left(Q+LP\right)u_{\infty }+\left(M{y}_{\mathrm{d}}+L{d}_{\infty }\right)\\ \therefore u_{\infty }={\left(I-Q-LP\right)}^{-1}\left(M{y}_{\mathrm{d}}+L{d}_{\infty }\right) \end{array}$$

□

Theorem 2.

In order to have ${\mathit{e}}_{\infty }=\mathbf{0}$, $\mathit{M}$, $\mathit{L},$ and $\mathit{Q}$ should be selected as

$$\begin{array}{c}Q=I \\ M=-L\end{array}$$

(41)

Proof of Theorem 2.

From (40), for $e_{\infty }=0$, we should have

$$\begin{gathered} \begin{array}{c}P{\left(I-Q-LP\right)}^{-1}L=-I\\ \begin{array}{c}I-Q-LP=-LP\\ \therefore Q=I\end{array}\end{array} \\ \begin{array}{c}P{\left(I-Q-LP\right)}^{-1}M=I\\ \begin{array}{c}I-Q-LP=MP\\ \begin{array}{c}I-I-LP=MP\\ \therefore -L=M\end{array}\end{array}\end{array} \end{gathered}$$

□

Although by selecting matrices as (41), a zero asymptotic error can be achieved, and the transient error also needs to be analyzed to prevent having significant transient errors in the system response. The controlled system is called monotonically convergent if

$$\begin{array}{c}\\ \Vert e_{\infty }-e_{\mathrm{k}+1}\Vert \le \mathsf{\gamma } \Vert e_{\infty }-e_{\mathrm{k}}\Vert \forall \mathrm{k}, 0\le \mathsf{\gamma }<1\end{array}$$

(42)

where $\Vert · \Vert$ is the Euclidean norm, and $\gamma$ is the convergence rate.

Theorem 3.

For the ILC law given in (34) we have,

$$\mathsf{\gamma }=\mathsf{\sigma }\left(\frac{I+LP}{2}\right)$$

(43)

where$\sigma (·)$denotes the maximum singular value operator.

Proof of Theorem 3.

From (40), we have

$$\begin{array}{cc}\hfill e_{\infty }-e_{\mathrm{k}+\mathrm{1}}=(I& -P{\left(I-Q-LP\right)}^{-1}M)y_{\mathrm{d}}-\left(P{\left(I-Q-LP\right)}^{-1}L+I\right)d_{\infty }-{\mathbf{y}}_{\mathrm{d}}\hfill \\ & +\mathbf{P}({\mathbf{Q}\mathbf{u}}_{\mathbf{k}}+\mathbf{L}({\mathbf{P}\mathbf{u}}_{\mathbf{k}}+{\mathbf{d}}_{\mathbf{k}})+{\mathbf{M}\mathbf{y}}_{\mathbf{d}}+{\mathbf{d}}_{\mathrm{k}+1}\hfill \\ & =(-P{\left(I-Q-LP\right)}^{-1}M+PM){\mathbf{y}}_{\mathrm{d}}+P(Q+LP){\mathbf{u}}_{\mathrm{k}}\hfill \\ & -\left(P{\left(I-Q-LP\right)}^{-1}L+I\right)d_{\infty }+{PLd}_{\mathrm{k}}+d_{\mathrm{k}+1}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill e_{\infty }-e_{\mathrm{k}}& =\left(-P{\left(I-Q-LP\right)}^{-1}M\right)y_{\mathrm{d}}+P{u}_{\mathrm{k}}-\left(P{\left(I-Q-LP\right)}^{-1}L+I\right)d_{\infty }+d_{\mathrm{k}}-{\mathbf{y}}_d+\mathbf{P}{\mathbf{u}}_k+{\mathbf{d}}_k\hfill \\ & =\left(-P{\left(I-Q-LP\right)}^{-1}M-I\right)y_{\mathrm{d}}+2P{u}_{\mathrm{k}}-\left(P{\left(I-Q-LP\right)}^{-1}L+I\right)d_{\infty }+2d_{\mathrm{k}}\hfill \end{array}$$

Using values obtained in (41)

$$e_{\infty }-e_{\mathrm{k}+1}=P\left(I+LP\right)u_{\mathrm{k}}-\left(I+PL\right)y_{\mathrm{d}}+PL{d}_{\mathrm{k}}+d_{\mathrm{k}+1}$$

and

$$e_{\infty }-e_{\mathrm{k}}=-2I{y}_{\mathrm{d}}+2P{u}_{\mathrm{k}}+2d_{\mathrm{k}}$$

Therefore,

$$\frac{\Vert e_{\infty }-e_{\mathrm{k}+1}\Vert }{\Vert e_{\infty }-e_{\mathrm{k}}\Vert }=\frac{\Vert -\left(I+PL\right)y_{\mathrm{d}}+P\left(I+LP\right)u_{\mathrm{k}}+PL{d}_{\mathrm{k}}+d_{\mathrm{k}+1}\Vert }{\Vert -2I{y}_{\mathrm{d}}+2P{u}_{\mathrm{k}}+2d_{\mathrm{k}}\Vert }$$

Considering that in ILC $\left\{x_k\left[0\right]=x_0 \forall k\right\}$, we have $d_{\mathrm{k}}=d_{\mathrm{k}+1}=C{A}^{\mathrm{k}} x_0 1^{\mathrm{T}}=d$ where $1$ denotes an all one $\mathbf{1}\times \mathrm{N}$ vector. Therefore,

$$\frac{\Vert e_{\infty }-e_{\mathrm{k}+1}\Vert }{\Vert e_{\infty }-e_{\mathrm{k}}\Vert }=\frac{\Vert -\left(I+PL\right)y_{\mathrm{d}}+P\left(I+LP\right)u_{\mathrm{k}}+\left(I+PL\right)d\Vert }{\Vert -2I{y}_{\mathrm{d}}+2P{u}_{\mathrm{k}}+2d\Vert }$$

Since for the matrices $P$ and $L$ given in (36) it can be proved that $\Vert \left(I+PL\right)\Vert =\Vert \left(I+LP\right)\Vert$, we have

$$\begin{array}{c}\frac{\Vert e_{\infty }-e_{\mathrm{k}+1}\Vert }{\Vert e_{\infty }-e_{\mathrm{k}}\Vert }\le \frac{\Vert \left(I+LP\right)\Vert }{2}\frac{\Vert -I{y}_{\mathrm{d}}+P{u}_{\mathrm{k}}+d\Vert }{\Vert -I{y}_{\mathrm{d}}+P{u}_{\mathrm{k}}+d\Vert }\\ \therefore \frac{\Vert e_{\infty }-e_{\mathrm{k}+1}\Vert }{\Vert e_{\infty }-e_{\mathrm{k}}\Vert }\le \frac{\Vert \left(I+LP\right)\Vert }{2}=\mathsf{\sigma }\left(\frac{I+LP}{2}\right)\end{array}$$

□

Using the results obtained from (39), (41), and (43), the ILC law given in (37) can achieve asymptotic stability as well as monotonic convergent and zero steady-state error when it is in the form of

$$\begin{array}{c}{u}_{\mathrm{k}+1}=u_{\mathrm{k}}+M\left(y_{\mathrm{d}}-y_{\mathrm{k}}\right)\\ \begin{array}{c}\mathsf{\rho }\left(I-MP\right)<1\\ 0 \le \mathsf{\sigma }\left(\frac{I-MP}{2}\right)<1\end{array}\end{array}$$

(44)

However, further consideration has to be taken into account before using (44) to design a controller for the pneumatic actuator. First, the considered pneumatic actuator is a nonlinear system with the state space equation of

$$\{\begin{array}{c}{\dot{x}}_{\mathrm{k}}\left(\mathrm{t}\right)=\mathrm{f}\left(x_{\mathrm{k}}\left(\mathrm{t}\right),{\mathrm{u}}_{\mathrm{k}}\left(\mathrm{t}\right)\right)\\ {\mathrm{y}}_{\mathrm{k}}\left(\mathrm{t}\right)=\mathrm{g}\left(x_{\mathrm{k}}\left(\mathrm{t}\right)\right) \end{array}; \mathrm{t}\in \left[0,\mathrm{T}\right] x_{\mathrm{k}}\left(0\right)=x_0 \forall \mathrm{k}$$

(45)

We assume that $\mathrm{f}:{ℝ}^{\mathrm{n}+1\to \mathrm{n}}$ and $\mathrm{g}:{ℝ}^{\mathrm{n}\to 1}$ are global Lipschitz continuous (GLC) functions. Therefore,

$$\begin{array}{c}\Vert \mathrm{f}\left(x_1,{\mathrm{u}}_1,\mathrm{t}\right)- \mathrm{f}\left(x_2,{\mathrm{u}}_2,\mathrm{t}\right)\Vert \le { \mathrm{L}}_{\mathrm{f}}\left(\Vert x_1-x_2\Vert +\Vert { \mathrm{u}}_1-{\mathrm{u}}_2\Vert \right)\\ \Vert \mathrm{g}\left(x_1,\mathrm{t}\right)- \mathrm{g}\left(x_2,\mathrm{t}\right)\Vert \le { \mathrm{L}}_{\mathrm{g}}\Vert x_1-x_2\Vert \end{array}$$

(46)

By taking into account that the considered pneumatic system is asymptotically stable and for a sampling time of ${\mathrm{T}}_{\mathrm{s}}$, where $\mathrm{T}={\mathrm{NT}}_{\mathrm{s}}$, the system can be considered as

$$\begin{array}{cc}\hfill x_{\mathrm{k}}\left[\mathrm{i}+1\right]& =A{x}_{\mathrm{k}}\left[\mathrm{i}\right]+B{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+ \mathsf{\xi }\left(x_{\mathrm{k}}\left[\mathrm{i}\right],{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]\right)\hfill \\ \hfill {\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]& =C{x}_{\mathrm{k}}\left[\mathrm{i}\right]+ \mathsf{\eta }\left(x_{\mathrm{k}}\left[\mathrm{i}\right]\right) ; x_{\mathrm{k}}\left[0\right]=x_0 \forall \mathrm{k},\mathrm{i}\in \left[0,\mathrm{N}\right]\hfill \end{array}$$

(47)

where $\underset{\Vert x\Vert \to 0}{\lim }\frac{\Vert \mathsf{\xi }\Vert }{\Vert x\Vert }=0$ and $\underset{\Vert x\Vert \to 0}{\lim }\frac{\Vert \mathsf{\eta }\Vert }{\Vert x\Vert }=0$. This follows the Lyapunov stability analysis, where it is assumed that for a physically stable system all the linear terms are in $A{x}_{\mathrm{k}}\left[\mathrm{i}\right]+B{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]$ and $C{x}_{\mathrm{k}}\left[\mathrm{i}\right]$ and higher-order terms are in $\mathsf{\xi }$ and $\mathsf{\eta }$, in the sense that when $x$ gets small, $\mathsf{\xi }$ and $\mathsf{\eta }$ get small faster. This means that the nonlinearities’ effects in the system’s dynamic would eventually vanish, and that the system has a dominating linear characteristic at its steady-state condition. The remaining nonlinearities is as a result of input effect, ${\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]$, which can be denoted as $\mathsf{\psi }\left({\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]\right)$. Using (47), the system’s output can be presented as

$$\begin{array}{cc}\hfill {\mathrm{y}}_{\mathrm{k}}\left[\mathrm{i}\right]& =C{\left(\mathrm{q}I-A\right)}^{-1} B{ \mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+C{A}^{\mathrm{k}} x_0+ \mathsf{\psi }\left({\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]\right)\hfill \\ & = \mathrm{P}\left(\mathrm{q}\right){\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]+{\mathrm{d}}_{\mathrm{k}}+ \mathsf{\psi }\left({\mathrm{u}}_{\mathrm{k}}\left[\mathrm{i}\right]\right)\hfill \end{array}$$

(48)

where $\mathsf{\psi }(·)$ term is as a result of the nonlinearities’ residue in the model. By applying the proposed ILC law given in (37), we have

$$\begin{array}{cc}\hfill u_{\mathrm{k}+1}& =Q{u}_{\mathrm{k}}+L{y}_{\mathrm{k}}+M{y}_{\mathrm{d}}\hfill \\ & =\left(Q+LP\right)u_{\mathrm{k}}+\left(M{y}_{\mathrm{d}}+L{d}_{\mathrm{k}}\right)+L\mathsf{\psi }\left(u_{\mathrm{k}}\right) \hfill \end{array}$$

(49)

In order to make sure that the ILC method is AS, an additional term has been added to the ILC law as

$$u_{\mathrm{k}+1}=Q{u}_{\mathrm{k}}+L{y}_{\mathrm{k}}+M{y}_{\mathrm{d}}-\mathsf{\alpha }{u}_{\mathrm{k}}$$

(50)

where $\Vert L\frac{\partial \mathsf{\psi }}{\partial u_{\mathrm{k}}}\Vert \le \mathsf{\alpha } \forall \mathrm{k}$. The existence of $\mathsf{\alpha }$ is guaranteed as $\mathrm{f}(·)$ and $\mathrm{g}(·)$ are assumed to be GLC functions. As the nonlinearities in the system can be controlled (removed) by adding the proportional control term, $-\mathsf{\alpha }{u}_{\mathrm{k}}$, the rest of the system can be considered linear, and (39), (41), and (43) will be held.

The ILC, in theory, is effective for repetitive processes. However, we would like to expand the ILC application to control the pneumatic system responding to non-repetitive inputs and disturbances. For this purpose, the proof given for Theorem 3 has to be revisited. As part of this proof, we had

$$\frac{\Vert e_{\infty }-e_{\mathrm{k}+1}\Vert }{\Vert e_{\infty }-e_{\mathrm{k}}\Vert }=\frac{\Vert -\left(I+PL\right)y_{\mathrm{d}}+P\left(I+LP\right)u_{\mathrm{k}}+PL{d}_{\mathrm{k}}+d_{\mathrm{k}+1}\Vert }{\Vert -2I{y}_{\mathrm{d}}+2P{u}_{\mathrm{k}}+2d_{\mathrm{k}}\Vert }$$

(51)

However, as the process is not repetitive, we cannot consider $d_{\mathrm{k}}=d_{\mathrm{k}+1}$. Instead, we should select the $\Vert PL\Vert$ value as its maximum, since this choice can help the condition to be asymptotically held as $\underset{\Vert PL\Vert \to \infty }{\lim }\frac{\Vert PL{d}_{\mathrm{k}}+d_{\mathrm{k}+1}\mathbf{\Vert }}{\Vert PL{d}_{\mathrm{k}}+d_{\mathrm{k}}\mathbf{\Vert }}=1$. Therefore, the ILC law for controlling the pneumatic actuator can be summarized as:

$$\begin{array}{c}{u}_{\mathrm{k}+1}=\left(1-\mathsf{\alpha }\right)u_{\mathrm{k}}+M\left(y_{\mathrm{d}}-y_{\mathrm{k}}\right)\\ \begin{array}{c}\mathsf{\rho }\left(I-MP\right)<1\\ \begin{array}{c}0 \le \mathsf{\sigma }\left(\frac{I-MP}{2}\right)<1\\ 1 \ll \Vert PL\Vert \end{array}\end{array}\end{array}$$

(52)

4. Results and Discussion

Table 1. Parameters values used in simulation.

	Parameter Name	Parameter Value
Gas Properties	Supply pressure (Ps)	7 × 105 (Pa)
	Supply temperature (Ts)	293.15 (K)
	Gas constant (R)	287 (J/(kg·K))
	Specific enthalpy (h)	293.6 (kJ/kg)
	Compressibility factor (z)	0.999
	Specific heat capacity (Cp)	1.01 (kJ/(kg·K))
	Thermal conductivity ($\mathsf{\kappa }$)	25.7 (mW/(K.m))
	Dynamic viscosity ($\mathsf{\mu }$)	18.2 ($\mathsf{\mu }$Pa·s)
	Specific heat ratio ($\mathsf{\gamma }$)	1.4
	Reference temperature (T0)	293.15 (K)
	Atmosphere pressure (Patm)	1 × 105 (Pa)
Valve Properties	Discharge coefficient (Cd)	0.82
	Max orifice area (Amax)	4 × 10−6 (m2)
	Leakage area (Aleakage)	1 × 10−10 (m2)
	Displacement for leakage area (xleakage)	2 × 10−4 (m)
	Displacement limit (xmax)	5 × 10−3 (m)
	Input voltage range	[–12, 12] (V)
Pipe Properties	Length (Lt)	1 (m)
	Pipe cross-sectional area (At)	5 × 10−6 (m)
	Internal surface absolute roughness (${ϵ}_{\mathrm{rough}}$)	15 × 10−6 (m)
	Pipe wall density (${\mathsf{\rho }}_{\mathrm{t}}$)	1500 (kg/m3)
	Pipe wall specific heat capacity (Ct)	1250 (J/(kg·K))
	Wall-air heat transfer coefficient (${\kappa }_{\mathrm{air}}$)	20 (W/(m2·K))
Cylinder Properties	Initial interface displacement (Linit)	0 (m)
	Max piston stroke (Lp)	0.2 (m)
	Interface cross-sectional area (Ap)	0.002 (m2)
	Dead volume (Vd)	4 × 10−5 (m3)
	Gas-wall heat transfer coefficient (${\mathsf{\kappa }}_{\mathrm{p}}$)	100 (W/(m2·K))
	Actuator wall specific heat (CCylinder)	870 (J/(kg·K))
	Actuator mass (MC)	3 (kg)
	Piston mass (MP)	1 (kg)
	Hard stop stiffness (khs)	1 × 107 (N/m)
	Hard stop damping (${\mathsf{\beta }}_{\mathrm{hs}}$)	1500 (N/(m·s))
	Mechanical damping (${\mathsf{\beta }}_{\mathrm{Mech}}$)	200 (N/(m·s))

presents the parameters’ values used in simulating the piston-cylinder pneumatic actuator.

The effect of applied voltage on the valve’s cross-sectional area, where the valve supplies port B for the negative voltage values and connects the supply port P to port A for the positive values, is shown in Figure 6. As is seen, a dead-band behavior appears in the vicinity of 0 V. This follows the pattern of the same in the physical actuator studied in [44], shown in Figure 6 of this reference.

The gas flow rate, ${\dot{\mathrm{m}}}_{\mathrm{in}}$, in the valve with respect to the applied voltage values is presented in Figure 7. As a result of leakage, the gas flow rate for the applied voltage less than ~0.5 V remains at zero. By applying more voltage, the gas flow rate increases by following a laminar regime until it reaches a point that the flow rate enters into a turbulence regime, in which a decrease in the flow rate slope can be observed. This follows the pattern of the same in the physical actuator studied in [44] shown in the Figure 7 of this reference.

We then apply an input voltage to the valve and measure the rod’s displacement. This is depicted in Figure 8. The input signal is selected in such a way as to take the valve’s spool to its extreme positions. An applied 12 V signal charges chamber 1 and allows gas to be exhausted from chamber 2, which causes the rod to extend to its maximum displacement, 0.2 m. The signal has been applied long enough so that the system can be settled. The valve is then moved to its neutral state by applying 0 V, and as is shown, the rod’s displacement remains unchanged. By applying −12 V to the valve, the gas in chamber 1 exhausts, and chamber 2 fills with gas, moving the piston inward until it reaches zero. Again, a 0 V signal has been applied to move the valve to its neutral state, where holding the rod’s position unchanged. The response transition time is around 2 s.

The frequency response of the system has also been studied. As was discussed, the air is compressible and has a low damping characteristic, which causes a nonlinear response and increases the pneumatic system’s dynamic order. Moreover, before the system can apply any force to a load, the pipes and cylinders have to be filled with air. This results in further nonlinearities in the form of dead-band and transmission attenuations. Moreover, a pneumatic system is affected by frictional forces caused by the mechanical parts’ movements. All these uncertainties and nonlinearities make linearizing a pneumatic system a complicated practice that generates inaccurate results. Therefore, a system identification approach has been implemented to estimate the frequency response of the pneumatic actuator. In this approach, a set of sinusoidal signals with different frequencies are applied to the system and the piston’s rod displacement is measured. The collected results are then used to draw the system’s bode plot. The estimated bode plot for the pneumatic actuator obtained by applying sinusoidal voltage signals in the range of 0.5 to 100 Hz is given in Figure 9, showing around 1.6 rad/s frequency bandwidth for the considered pneumatic actuator.

According to [45], servo pneumatic actuators in approximately 70% of industrial applications should move 1–10 kg payloads with ±2 to ±0.02 mm precision. Therefore, in this paper, we consider the same criteria in studying the performance of the designed controllers. One major advantage of ILC is that it does not need accurate knowledge of the system model, and instead, it can learn from the system’s historical input and output. However, some degree of estimation can help to obtain a better design with a faster convergence. For this purpose, the estimated bode plot, given in Figure 9, will be used. As is seen, the system is a lowpass with the bandwidth of 1.6 rad/s. Therefore, we use

$$\mathrm{P}\left(\mathrm{s}\right)=\frac{0.01}{\mathrm{s}+1}$$

(53)

to estimate the linear part of the system. Since the system’s bandwidth is 1.6 rad/s, it can be assumed that the system does not have much of fluctuations over a period of one second. Therefore, by taking $\mathrm{T} =0.01 \mathrm{s}$, we can consider that the system would be seen as a repetitive process from the controller’s perspective. The $\mathrm{P}\left(\mathrm{q}\right)$ for the system given by (53) based on $\mathrm{T}=0.01 \mathrm{s}$ is equal to

$$\mathrm{P}\left(\mathrm{q}\right)=9.95\times {10}^{-5}\times \left({\mathrm{q}}^{-1}+\left(0.99\right){\mathrm{q}}^{-2}+{\left(0.99\right)}^2{\mathrm{q}}^{-3}+\dots \right)$$

(54)

and

$$P=9.95\times {10}^{-5}\times \left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}1\\ \left(0.99\right)\end{array}\\ \begin{array}{c}⋮\\ {\left(0.99\right)}^{\mathrm{N}-1}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 1\end{array}\\ \begin{array}{c}⋮\\ {\left(0.99\right)}^{\mathrm{N}-2}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}⋮\\ 1\end{array}\end{array}\end{array}\end{array}\right]$$

(55)

$M$ should be selected such that $\mathsf{\rho }\left(I-MP\right)<1$ and $0 \le \mathsf{\sigma }\left(\frac{I-MP}{2}\right)<1$. Theoretically, $M=P^{-1}$ will perfectly satisfy both conditions. However, as

$$P^{-1}=\frac{1}{9.95\times {10}^{-5}}\times \left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}1\\ -0.99\end{array}\\ \begin{array}{c}\begin{array}{c}0\\ ⋮\end{array}\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 1\end{array}\\ \begin{array}{c}\begin{array}{c}-0.99\\ ⋮\end{array}\\ 0\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\dots \\ \dots \end{array}\\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}\begin{array}{c}⋮\\ 0\end{array}\\ 1\end{array}\end{array}\end{array}\end{array}\right]$$

(56)

this design would take the system into saturation due to its high gain value of $1/9.95\times {10}^{-5}$. Considering the maximum rod’s displacement as 0.2 m and the maximum input voltage as 12 V, the maximum gain value should be limited to $12/0.2=60.$ Another observation from (56) is that the maximum number of non-zero elements in each column is two, which means that regardless of the selected value for $\mathrm{N}$, only two consecutive samples are used in the control law. As a result, $\mathrm{N} =2$ would be sufficient for implementing the ILC to control this pneumatic actuator. Therefore, we consider $M$ as a $2\times 2$ lower triangular matrix and calculate its arrays’ values.

$$\begin{array}{l} M=\left[\begin{array}{cc}\begin{array}{c}{\mathrm{m}}_0\\ {\mathrm{m}}_1\end{array}& \begin{array}{c}0\\ {\mathrm{m}}_0\end{array}\end{array}\right]P=\left[\begin{array}{cc}\begin{array}{c}{\mathrm{P}}_1\\ {\mathrm{P}}_2\end{array}& \begin{array}{c}0\\ {\mathrm{P}}_1\end{array}\end{array}\right]\\ \det {\displaystyle (}\mathsf{\lambda }I-{\displaystyle (}I-MP{\displaystyle )}{\displaystyle )}=\det \left(\begin{array}{cc}\mathsf{\lambda }-1+{\mathrm{m}}_0{\mathrm{P}}_1& 0\\ {\mathrm{m}}_0{\mathrm{P}}_2+{\mathrm{m}}_1{\mathrm{P}}_1& \mathsf{\lambda }-1+{\mathrm{m}}_0{\mathrm{P}}_1\end{array}\right)=0\hfill \\ \mathsf{\lambda }=1-{\mathrm{m}}_0{\mathrm{P}}_1; \mathsf{\rho }\left(I-MP\right)=1-{\mathrm{m}}_0{\mathrm{P}}_1\hfill \\ \therefore 0\le {\mathrm{m}}_0\le \frac{1}{{\mathrm{P}}_1}=10050\hfill \\ \mathsf{\sigma }\left(\frac{I-MP}{2}\right)=\frac{I-MP}{2}=\frac{\sqrt{2{\left(1-{\mathrm{m}}_0{\mathrm{P}}_1\right)}^2+{\left({\mathrm{m}}_0{\mathrm{P}}_2+{\mathrm{m}}_1{\mathrm{P}}_1\right)}^2}}{2}\hfill \\ \therefore \frac{-{\mathrm{m}}_0{\mathrm{P}}_2-\sqrt{4-2{\left(1-{\mathrm{m}}_0{\mathrm{P}}_1\right)}^2}}{{\mathrm{P}}_1}\le {\mathrm{m}}_1\le \frac{\sqrt{4-2{\left(1-{\mathrm{m}}_0{\mathrm{P}}_1\right)}^2}-{\mathrm{m}}_0{\mathrm{P}}_2}{{\mathrm{P}}_1}\hfill \end{array}$$

(57)

Although from the theoretical perspective, a larger value of ${\mathrm{m}}_0$ improves the AS condition of the ILC method, from the practical aspect, it should be limited to 60 to prevent the system from saturation. Figure 10 shows the relation between $\mathsf{\sigma }\left(\frac{I-MP}{2}\right)$ and ${\mathrm{m}}_1\in \left[-60, 60\right]$ to avoid saturation. From a repetitive process perspective, the best value would be ${\mathrm{m}}_1=-60$ to minimize the transient error as is given by (43). However, the inputs to the considered pneumatic system are not repetitive, and following (52), the values in (57) are selected as ${\mathrm{m}}_0={ \mathrm{m}}_1=60$, which makes $\mathsf{\rho } =0.994$ and $\mathsf{\sigma } =0.5$, satisfying the conditions in (52). The value for $\mathsf{\alpha }$ should be in the range of $\left[0, 1\right)$ to satisfy the AS condition of the ILC method, which we adjusted to 0.25 in this design. However, choosing the optimal value for $\mathsf{\alpha }$ requires the system’s knowledge, which is not available.

Based on the above discussions, the designed ILC method for the pneumatic actuator system is implemented as

$${\mathrm{u}}_{\mathrm{k}+1}\left[\mathrm{n}\right]=0.75{\mathrm{u}}_{\mathrm{k}}\left[\mathrm{n}\right]+\left(60+60{\mathrm{q}}^{-1}\right)\left({\mathrm{y}}_{\mathrm{d}}\left[\mathrm{n}\right]-{\mathrm{y}}_{\mathrm{k}}\left[\mathrm{n}\right]\right)$$

(58)

We further used the obtained frequency response of the system, given in Figure 9, to design a PID controller in order to compare its performance to that of the proposed ILC method. The frequency domain techniques, namely the Nichols chart and Inverse Nichols chart, have been used to achieve the required performance in reference tracking and disturbance rejection. To obtain the precision of ±0.002 m, the tracking boundaries are designed to be in the range of

$$\begin{array}{c}20\log \left(\frac{0.2-0.002}{0.2}\right) \le { \mathrm{TB}}_{\mathrm{dB}} \le 20\log \left(\frac{0.2+0.002}{0.2}\right)\\ -0.09 \mathrm{dB} \le { \mathrm{TB}}_{\mathrm{dB}} \le 0.08 \mathrm{dB}\end{array}$$

(59)

and the sensitivity bound is designed to be less than 3 dB for the lower frequencies ($\mathsf{\omega } <8 \mathrm{rad}/\mathrm{s}$) and less than 6 dB for all frequencies. The obtained PID controller is calculated as

$$\mathrm{C}\left(\mathrm{s}\right)=100+\frac{10}{\mathrm{s}}+\frac{8}{1+8\frac{1}{\mathrm{s}}}$$

(60)

Moreover, two recently proposed intelligent model-free control approaches have been selected to compare the performance of the designed ILC method with theirs. The first approach, proposed in [46], uses type-2 Takagi–Sugeno (T-S) fuzzy systems with the memory state feedback control. The controller is expressed using the state variable $\mathrm{x}\left[\mathrm{n}\right]$ as

$$\begin{array}{c}\mathrm{u}\left[\mathrm{n}\right]={\displaystyle {\displaystyle \sum }_{\mathsf{\mu }=1}^{\mathrm{q}}}{\mathsf{\eta }}_{\mathsf{\mu }}\left(\mathrm{x}\left[\mathrm{n}\right]\right)\left\{{\mathcal{K}}_{\mathsf{\mu }}\mathrm{x}\left[\mathrm{n}\right]+{\mathcal{G}}_{\mathsf{\mu }}\mathrm{x}\left[\mathrm{n} -1\right]\right\}\\ \begin{array}{c}\mathrm{Controller} \mathrm{rule} \mathsf{\mu }:\mathrm{IF} {\mathcal{W}}_1\left(\mathrm{x}\left[\mathrm{n}\right]\right) \mathrm{is} {\mathbb{W}}_1^{\mathsf{\mu }} \mathrm{and} {\mathcal{W}}_2\left(\mathrm{x}\left[\mathrm{n}\right]\right) \mathrm{is} {\mathbb{W}}_2^{\mathsf{\mu }} \\ \begin{array}{c}\mathrm{and}\dots \mathrm{and} {\mathcal{W}}_{\mathsf{\upsilon }}\left(\mathrm{x}\left[\mathrm{n}\right]\right) \mathrm{is} {\mathbb{W}}_{\mathsf{\upsilon }}^{\mathsf{\mu }} \mathrm{THEN}\\ \mathrm{u}\left[\mathrm{n}\right]={\mathcal{K}}_{\mathsf{\mu }}\mathrm{x}\left[\mathrm{n}\right]+{\mathcal{G}}_{\mathsf{\mu }}\mathrm{x}\left[\mathrm{n}-1\right]\end{array}\end{array}\end{array}$$

(61)

where ${\mathbb{W}}_{\mathsf{\upsilon }}^{\mathsf{\mu }}$ denotes the interval type-2 (IT2) fuzzy set of the premise variable, ${\mathcal{W}}_{\mathsf{\upsilon }}\left(\mathrm{x}\left[\mathrm{n}\right]\right)$, and $\mathsf{\mu }$ and $\mathsf{\upsilon }$ are the number of fuzzy IF-THEN rules and fuzzy sets of controller, respectively. ${\mathsf{\eta }}_{\mathsf{\mu }}\left(\mathrm{x}\left[\mathrm{n}\right]\right)$ are nonlinear weight coefficient functions that are calculated based on upper and lower membership functions. The footprint of uncertainty (FOU) for the employed T-S fuzzy system is given in Figure 11. Two fuzzy rules are defined taking the system error as the state variable and the coefficients are selected according to a procedure explained in the reference as ${\mathsf{\eta }}_1\left(\mathrm{x}\left[\mathrm{n}\right]\right)={\cos }^2\left(\mathrm{x}\left[\mathrm{n}\right]\right)$ and ${\mathsf{\eta }}_2\left(\mathrm{x}\left[\mathrm{n}\right]\right)=1-{\mathsf{\eta }}_1\left(\mathrm{x}\left[\mathrm{n}\right]\right)$. The feedback control gains are then calculated as $\left[{\mathcal{K}}_1,{\mathcal{K}}_2\right]=\left[5.8123, -5.8122\right]$ and $\left[{\mathcal{G}}_1,{\mathcal{G}}_2\right]=\left[9.1279, -9.1278\right]$.

In the other approach proposed in [47], a type-1 T–S fuzzy observer is used to develop a piecewise control strategy for dealing with the nonlinear terms in the system. The membership functions employed for constructing the observer is depicted in Figure 12. The controller gain is calculated as [5.37, −5.37].

The performance of the controllers in tracking a reference signal under a constant load of 10 kg is shown in Figure 13. As is seen, the designed ILC is capable of tracking a non-repetitive signal. Next, we compare the performance of the ILC with the other controllers, which is given in

Table 2. Comparison between the controllers’ performance in tracking a signal.

Controller	Rise Time (s)	Overshoot	Settling Time	Error
ILC	0.2	13.5%	0.52	0.0044
PID	0.26	0	0.46	0.0049
IT2 T-S fuzzy	0.24	0	0.35	0.0048
T-S fuzzy obs.	0.38	0	--	0.0081

. It can be seen that the ILC demonstrates a faster response, however, it demonstrates overshoots in the transient response. The maximum overshoot percentage happens during the forward motion when the actuator moves from the resting to reach 0.1 m, which equals 13.5%. This can be considered as the major drawback of the ILC, which is caused by the way that ILC calculates the control signal. The ILC only modifies the control signal according to predefined control law and based on the historical input and output. Therefore, a sudden change in the input signal will cause overshoot in the control signal. However, in a well-designed ILC that holds monotonic convergent, the gap between the input and output is reduced after sufficient iterations. The ILC can obtain ±0.002 m precision (2% settling time) in 0.52 s compared to 0.46 s for the PID and 0.35 s for the IT2 T_S fuzzy controller. The controller based on T-S fuzzy observer could not achieve the required precision. The difference between the reference and output of the system controlled by each of these controllers is shown in Figure 14, and the summation of the square of the error (${\displaystyle \int }{\mathrm{e}}^2\left(\mathrm{t}\right)\mathrm{dt}$) is given in Table 2. The ILC demonstrates superior performance than other controllers in terms of both speed (rise time) and tracking accuracy.

Next, we compare the performance of the controllers in overcoming the payload uncertainty, as is shown in Figure 15. The top plot in this figure shows the payload variation applied to the piston while the reference signal remains constant at 0.1 m. As is seen, the ILC-controlled system can maintain the required precession (±0.002 m) regardless of changes in the load. However, this cannot be adequately achieved with the other controllers. The plot of the disturbance rejection error is given in Figure 16. The summation of the square of the error for the controllers is given in

Table 3. Comparison between the controllers’ performance in disturbance rejection.

Controller	Error
ILC	$2.6\times {10}^{-5}$
PID	$3.3\times {10}^{-4}$
IT2 T-S fuzzy	$3.8\times {10}^{-4}$
T-S fuzzy obs.	$9.2\times {10}^{-4}$

. Again, the ILC shows a better performance than the other controllers in overcoming the uncertainties.

5. Conclusions

This paper suggested an effective model-free control scheme to overcome the nonlinearities and uncertainties resulting from air characteristics, pressure loss, leakage, and load variations in a pneumatic system. In the majority of the control algorithms used in controlling the position of a pneumatic actuator, the model of the system has to be achieved prior to the design of the controller. However, due to the physical behavior of the gases, modelling a pneumatic system is usually done based on many assumptions that might result in an inaccurate model during system operation. As a result, such controllers may not achieve the required performance. Moreover, intelligence-based control algorithms cannot easily be implemented in real-time applications due to their extensive calculation requirements. This leads to a need for developing real-time control systems capable of controlling pneumatic systems without needing to obtain the mathematical model of the system. An ILC algorithm uses information from previous repetitions to learn about the system’s dynamics for generating a more suitable control signal. ILC algorithms are particularly useful in real-time control systems, given their relatively quick response to changes in the input signal. However, their application is only limited to repetitive processes, where the same control action should be performed repeatedly. In this paper, the application of the ILC algorithm has been expanded for controlling nonlinear, non-repetitive systems. Moreover, a pneumatic cylinder-piston actuator has been simulated using MATLAB SimScape blocks and the simulation results showed that the model can successfully demonstrate the behavior of a pneumatic system. The procedure for designing an ILC controller without needing any information on the system model and by only using the system’s input and output measurements has been covered. The simulation results showed that the designed ILC controller is capable of tracking a non-repetitive reference signal and can overcome the internal and payload uncertainties with the precision of 0.002 m. Therefore, the ILC can be considered as an approach for controlling the pneumatic actuators, which is challenging to obtain their mathematical modeling.