An Innovative MIMO Iterative Learning Control Approach for the Position Control of a Hydraulic Press

To improve the performance of hydraulic press position control and eliminate the need to manually define control signals, this paper proposes a multi-input-multi-output (MIMO) Iterative Learning Control (ILC) algorithm. The MIMO ILC algorithm design is based on the inversion of the known low frequency dynamics of the hydraulic press, whereas the unknown and uncertain high frequency dynamics are discarded due to their low influence in the learning transient. Moreover, for the MIMO ILC convergence condition, a graphical method is proposed, in which the ILC learning filter eigenvalues are analyzed. This method allows studying the stability and convergence rate of the algorithm intuitively. Theoretical analysis and results prove that with the MIMO ILC algorithm the position control is automated and that high precision in the position tracking is gained. A comparison with other model inverse ILC approaches is carried out and it is shown that the proposed MIMO ILC algorithm outperforms the existing algorithms, reducing the number of iterations required to converge while guaranteeing system stability. Furthermore, experimental results in a hydraulic test rig are presented and compared to those obtained with a conventional PI controller.


I. INTRODUCTION
Hydraulic presses have traditionally been used for high-force applications due to their easy operation and adaptability to suit a wide range of forming conditions. One of the advantages of hydraulic presses is that the force can be applied throughout the stroke, unlike mechanical presses where the force is maximum at the bottom.
During the hydraulic press operation, shown in Fig. 1, position and force control needs to be done to guarantee the correct forming of the workpiece. The position control ensures the correct trajectory tracking during the Free fall and Drawing phases, whereas the force control maintains a specific force level during the Making Force phase.
Nowadays, in the working operation of a hydraulic press the proportional valve control signals are defined manually, The associate editor coordinating the review of this manuscript and approving it for publication was Wonhee Kim . with a predefined closing ramp, so the position is correctly tracked during the Free fall and Drawing phases. This is a tedious and time-consuming process that typically includes several days and can involve high costs. Furthermore, the success of such a control depends on the valve pre-defined signals design, which usually depends on the operator's ability and experience.
We seek to automate the rather tedious manual position control, with no need to manually define any signal. Nowadays, the desired falling velocity during the Free Fall phase is obtained by manually setting a constant spool position to the auxiliary chamber proportional valve, see Fig. 2. Furthermore, during the Free Fall and Drawing phases transition, a closing ramp must be defined to this valve which is modified every time the position reference changes.
We propose to control the auxiliary chamber pressure with the variable axial piston pump, while controlling the slide position with the auxiliary chamber proportional valve.  This changes the press position control from a singleinput-single-output (SISO) problem to a multiple-inputmultiple-output (MIMO) problem, as the pump and the valve control the cylinder pressure and position, respectively. The MIMO position controller will not differentiate between phases, as the control will be carried out continuously from the start of the Free Fall phase to the end of the Drawing phase, regardless of the position reference.
It is necessary to guarantee a certain pressure level in the auxiliary chamber during the Drawing phase, as the forming of the workpiece is carried out in this phase. In the event that more force is required for the workpiece forming, if the auxiliary chamber is pressurized at a certain level it will take less effort for the pump to supply said force. Else, if no pressure existed in the auxiliary chamber, the main chamber would have to be pressurized completely before more force could be applied.
Therefore, in order to be prepared for a higher force demand during the workpiece forming, maintaining a specific pressure level in the auxiliary chamber is also of great importance.
However, a hydraulic circuit is highly nonlinear regarding flow and pressure relationship, friction, oil temperature variation, component coupling, etc. These nonlinear behaviors make the process of designing a combined position and force control algorithm a challenge.
In [1], fuzzy laws are designed to overcome the coupling effects and to improve the position and force control performance of an electro-hydraulic system. Although successful results were obtained, we aim to design a self-learning controller that adapts automatically to every force and position scenario, eliminating the need to manually redefine control signals as is the case with fuzzy logic.
A decoupling compensation method is proposed in [2] for an electro-hydraulic circuit. The study focuses on the pressure coordinate control with the pump and valve, rather than decoupling the valve and pump for position and force control.
In [3], a single-input-multiple-output (SIMO) PI control system is developed for the force applied by 12 actuators in a stamping machine. The SIMO PI controller gains are designed based on the estimated perturbation model, however, as the authors point out, the fine-tuning of the controller could result in time-consuming and expensive.
In [4] the pressure and velocity of a hydraulic actuator are decoupled by means of velocity feedback. In [5], a feed-forward control scheme is presented for decoupled pressure control in an actuator. In [6], a MIMO inversion-based feed-forward controller is derived by input-output linearization for a hydraulic system. However, these works focus on the energy efficiency of the hydraulic system rather than on the control performance.
Considering the limitations of the existing MIMO controllers in hydraulic systems and taking advantage of the repetitive process of a hydraulic press, we try Iterative Learning Control (ILC) as a solution to the MIMO position control problem. With the MIMO ILC algorithm, the valve and pump control loops coupling will be counteracted and, as a result, a self learning MIMO position and force control will be obtained.
Several studies have been realized regarding the design of MIMO ILC algorithms. In [7], a conjugate-gradient ILC was proposed to guarantee fast monotonic convergence without the need to calculate the inverse of the system matrix. This method converged fast at the first iterations, but the algorithm performance suffered from the control loop couplings, penalizing the convergence rate later on. This was improved in [8], with a Quasi-Newton optimized ILC. However, a gain was included in the learning function to counteract for model mismatches, which penalized the convergence rate of the algorithm.
To improve the convergence rate, optimization based MIMO algorithms have been proposed, in which the learning function is based on an error minimization function [9]- [11]. These approaches rely on the tuning of the weighting matrix in order to obtain a global minimum in the minimization problem. It is hard to explicitly define the ILC multiple objectives in a cost function. In such a cost function we would have to include a condition for fast convergence, robustness towards uncertainties, long-term stability and satisfy the input constraints. This approach requires a too complex optimization problem.
To reduce the model dependent design, model-free MIMO ILC approaches were proposed in [12], [13], to reduce the modeling requirement. However, these approaches directly affect the convergence rate of the ILC algorithm, and the learning process requires more iterations to converge.
The advantages of using model-based approaches were shown in [14], where an overview of which MIMO ILC design approach to use, illustrated on an industrial printer, is given. It concludes that the centralized, modelbased, approach achieves the better performance in terms on convergence, however it requires more designing and computation effort, in comparison to the decentralized, model-free, approach. In the position control problem of a hydraulic press, we can manage the modeling designing effort better than a slow convergence. Investing in an accurate plant model design can optimize the convergence of the algorithm, which results in commissioning cost savings. Else, the slower the convergence, the more iterations will be required to obtain an appropriate workpiece forming that fulfills all the design requisites.
So far, no studies have been done regarding the implementation of a MIMO ILC algorithm in hydraulic systems. Previous works have just focused on the application of SISO ILC algorithms either for the position control or for the force control of hydraulic circuits.
An ILC algorithm for the position control of an electro-hydraulic system is developed in [15]. The ILC algorithm improves the PID controller performance, however no analysis of the robustness and the convergence rate of the algorithm is provided.
In [16] the tuning of the learning gain and the delay time is analyzed. Increasing the learning gain value the convergence rate is improved, however the system response oscillates harshly. Increasing the delay time the tracking error is decreased, but as iterations go on the system becomes unstable.
Most of the ILC studies in hydraulic systems do not consider the convergence rate of the algorithm, the focus has been on the applicability of ILC in hydraulic systems, no matter the number of iterations required to converge.
One should aim to obtain the fastest convergence rate possible to reduce the number of iterations needed. The fewer iterations required, the shorter the commissioning time of the hydraulic press control will be, which translates into a reduction of the production costs.
The convergence rate of the ILC is improved in [17], to track the displacement curve of the hydraulic press slide. A fuzzy ILC is developed, to adapt the ILC gains based on the error and change in the error. The requirement of defining fuzzy strategies manually opposes the objective of this study to automate the control process.
A PD-type ILC is implemented for the trajectory tracking in [18]. The number of iterations required to converge is considerably reduced compared to a conventional proportional type ILC, as it is also shown in [19]. However, no analysis of how to design the PD-type ILC gains is provided. The only procedure followed is that the higher the gain the faster the convergence will be. However, this design method can lead to stability problems as if the gains are set too high the learning transients will not converge.
One can conclude that there exists a gap in the literature regarding the design of MIMO ILC algorithms in hydraulic systems. No analysis of the convergence rate improvement of the MIMO ILC algorithm to optimize the press working operation has been done. Furthermore, this improvement should be enhanced with a stability analysis of the algorithm, as stability is an essential property for industrial practice.
The main contributions of this paper are summarized as follows. A MIMO ILC is proposed for the position control of a hydraulic press, to automate the press operation and eliminate the position control reliance on pre-defined control inputs. With the MIMO ILC controller, we can control simultaneously the pump outflow and the valve spool position, to control the cylinder pressure and piston position, respectively. The MIMO ILC algorithm design is based on the inversion of the low frequency hydraulic press dynamics, discarding the high frequency dynamics. This design method obtains a better convergence rate compared to other model inverse ILC algorithms in the literature. Furthermore, a graphical method is proposed which provides the user the possibility to analyze the convergence rate and the stability of the MIMO ILC algorithm intuitively.
In Section II, the MIMO ILC algorithm is introduced for the position control. The results of the MIMO ILC algorithm are shown in Section III. A comparison of the proposed ILC design with other ILC methods is carried out in Section IV. The proposed MIMO ILC is implemented in a hydraulic test rig in Section V. Lastly, the conclusions are provided in Section VI.

II. MIMO ITERATIVE LEARNING CONTROL
ILC is ideal for a system that performs the same process repeatedly and where uncertainties abound, as it is capable of learning the necessary control action via the repeated process. ILC was first proposed for improving the reference tracking of a system that follows a repeated trajectory by Uchiyama in [20]. It was further extended by Arimoto et al. [21], for a mechanical robot operation. The learning control scheme proposed by Arimoto was: which is a past-error feed-forward law. U j (s) is the Laplace transform of the entire input vector to the system G(s) at the j-th learning iteration. The iteration error E j (s) is given by the difference between the position reference and the iteration system output, filter and Q(s) is a gain that enlarges the stability region of the algorithm.

A. MIMO ILC ALGORITHM ANALYSIS
In the position control loop, the slide hydraulic circuit can be regarded as a general time invariant MIMO linear continuous state-variable equation: where x is the four-dimensional state vector, consisting of the piston position and velocity, and the pressure inside the two cylinder chambers. u is the two-dimensional learning control input vector, consisting of the pump swash angle and the valve spool position. y is the two-dimensional output vector, consisting of the piston position and the cylinder auxiliary chamber pressure. The Laplace transform of the system equation is: G(s) = C s (sI − A s ) −1 B s . Note that G(s) is a twoby-two transfer matrix. In [22] the ILC error propagation equation for a SISO system is shown. Transforming it into a MIMO problem results in: where S(s) is the sensitivity transfer matrix S(s) = (I + G(s)C(s)) −1 , and L(s) the learning transfer matrix. C(s) is the controller transfer matrix containing the pump PI controller, C P (s), and the valve PI controller, C V (s): From equation (3), if the term (s) = I − S(s)G(s)L(s) is less than one, then the error vector at the next iteration will be smaller than the error at the current iteration. Otherwise, the error would increase causing instability in the system. This brings a sufficient condition for the stability and the convergence of the error vector for the design of L(s), which is given for the frequency domain by [23]: The tracking error will converge to zero as iterations go on if (s) has all its eigenvalues, λ i , less than one.
From equation (5), the learning transfer matrix L(s) must be designed so the eigenvalues values of (s) are less than one. The two-by-two L(s) learning transfer matrix will have the following structure: For the L(s) design the plant dynamics inversion is a benchmark method to achieve convergence in one iteration. However, the uncertainty present in plant models can lead to sub-optimal tracking, stability issues and convergence rate penalization. Therefore, one must accept that without an exact plant model perfect tracking of the reference is unattainable.
Several studies have addressed the incorporation of model dynamics in the L(s) design, to increase the convergence rate of the ILC algorithm. The convergence condition and the robustness of an inverse model-based ILC are analyzed in [24]. A quadratic criterion is introduced to analyze the convergence on a model-based ILC in [25]. The optimization problem for the ILC controller is also developed in [26], [27]. As it has been explained in Section I, it would be too laborious to solve the optimization problem of a cost function that considers the algorithm convergence rate, stability and input constraints.
For systems with unstable zeros, pseudo-inverse and stable inversion methods have been proposed in [28] and [29], respectively. A thorough analysis of which inversion approach to use depending on the system characteristics is carried out in [30]. The hydraulic slide circuit in Fig. 2 is a minimum-phase system, therefore, direct inversion can be applied.
In [31] direct inversion was applied in an electronic printer, where the learning filter was designed as L(s) = (G(s)S(s)) −1 . Instead of including a low-pass filter in the L(s) design, the Q(s) filter was used to guarantee robustness to modeling errors. However, Q(s) should be designed properly, as a Q(s) with a module different from one penalizes the algorithm performance [32].
Direct inversion was also applied in [33], for a linear SISO motor of a wafer system. In [33], it is proposed to include the controller in the model inverse ILC design based on the series and parallel ILC designs. The series, P s (s), and parallel, P p (s), ILC architectures have the following model transfer functions: P s (s) and P p (s) are the process sensitivity functions, which are used as a parametric model of the system to design the learning filter. In [33], unlike in [31], the low-pass filter was applied to the inverse of the process sensitivity functions to deal with the high frequency differences. This is the most common procedure that one can find in the literature when applying direct inversion [33]- [36]. VOLUME 9, 2021 Although fast convergence tracking is achieved, the algorithm performance is penalized as the inclusion of the low-pass filter to counteract model mismatches affects both G(s) and C(s).
We claim that it is not necessary that the inclusion of the low-pass filter in the learning filter design affects C(s). One can design the learning filter so that the model differences only affect the inverse of G(s). For MIMO systems with a parallel structure consisting of an ILC algorithm and feedback controllers, we propose the following model inverse L(s) design: In (8), only G(s) is inverted and C(s) is completely known, see (4). The low-pass filters to deal with high frequency differences will only apply to G −1 (s), not affecting C(s). With this L(s) design, the convergence rate is improved with respect to the traditional model inverse approach, as the model differences are reduced.
Furthermore, in the G(s) inversion we only include the plant low frequencies that we are interested in, discarding the high frequencies. The low frequencies provide us with information about the pressure change due to the velocity, which is necessary for the position and force control. The hydraulic press model uncertainties abound at high frequencies, such as pump and valve dynamics and friction forces, these uncertain high frequencies will be attenuated with zero-phase filtering.

B. HYDRAULIC MODELING
In order to carry out the model inverse L(s) design, the hydraulic press modeling is derived in this section.
The axial-piston pump outflow is proportional to the swash plate angle: where q A is the volumetric flow out of the pump (m 3 /s), q N is the nominal flow (m 3 /s), ω the shaft rotational speed (rad/s), ω N the nominal shaft rotational speed (rad/s) and α is the swash plate angle. The terms q N and ω N are used to normalize α within [0-1] interval. Equation (10) shows the relationship between the pressure into the cylinder main chamber and the volumetric flow out of the pump [37]: where A A is the cylinder main chamber area (m 2 ), V A is the main cylinder chamber dead volume (m 3 ), β is the hydraulic compressibility (1/bar), x is the piston position (m) and P A is the main chamber pressure (bar). The relationship between the pressure in the cylinder auxiliary chamber and the volumetric flow out of the chamber is as follows [37]: where q B is the volumetric flow out of auxiliary chamber (m 3 /s), A B is the auxiliary chamber area (m 2 ), V B is the auxiliary chamber dead volume (m 3 ), l is the cylinder piston length (m) and P B is the auxiliary chamber pressure (bar). The relationship between the pressure in the auxiliary chamber and the flow through the valve is shown [38]: (12) where K v (y v ) is the hydraulic conductance (m 3 /(bar·s)), which is function of the valve spool position y v . q ref and P ref are included to normalize K v (y v ) with respect to the nominal flow and the nominal pressure of the valve, respectively.
The force balance equation of the cylinder is: where m is the cylinder moving mass (kg) and F is the output force of the rod (N). During the Free Fall and Drawing phases, the oil flow rate into the cylinder main chamber will be equal to the pump outflow. Likewise, the oil flow rate out of the auxiliary chamber will equal the oil flow through the proportional valve.
Rearranging (9) with (10), (11) with (12), and using (13) to relate both cylinder chamber pressures, we can obtain the slide hydraulic circuit system equations, from which the linearized state-space system is derived in equation (14): The cylinder piston position and velocity, and the pressure in both chambers are chosen as state variables. The swash plate angle and the valve spool position are the system inputs; the piston position and the main chamber pressure are the system outputs.
are the small-signal deviation from an operating point obtained in steady-state conditions,P x ,x,ȳ v andᾱ respectively. The definition of every term A xx is shown in Appendix A.

C. HYDRAULIC PRESS MIMO ILC DESIGN
By Laplace transforming equation (14), the system transfer matrix G(s) is obtained, from which the system inverse will be calculated to carry out the L(s) = G −1 (s) + C(s) design. G(s) is a two-by-two transfer matrix with four transfer functions with a common denominator, from which the system poles can be obtained. To analyze how the system parameters affect the poles of the system, the G(s) denominator based on the system parameters is obtained. The denominator is a fourth-order polynomial, each coefficient is shown in Table 1 We can put the polynomial in Table 1 into numbers, with the values provided by the data-sheets, shown in Appendix B, and obtain the poles of our system: From Table 2 we can see that we have a pole in zero, a low frequency pole and a high frequency complex pole. The high frequency poles have really small damping, as the hydraulic slide energy does not dissipate anywhere as no friction has been included in the system model, e.g. friction between the slide and the rails or between the cylinder piston and the chambers. Although some authors include the friction forces in the force balance equation of the cylinder, see [39], we will discard them as they will only affect our plant at high frequencies. The low frequency pole describes the cylinder pressure change due to the velocity.
In the G(s) inverse design, we will not include the high frequencies, as these frequencies are where the uncertainties of our plant abound. The effects of the oil compressibility appear at high frequencies, so we could discard every term containing β in the coefficients shown in Table 1. However, the terms containing β and P B provide us with information of the rate at which the pressure is increasing, necessary for the pressure control.
As shown in Appendix B, β is a small number. In Table 1, in the fourth-order coefficient β appears squared, so the value of the coefficient will be insignificant in comparison to the others. In the third-order coefficient, β does not multiply any pressure variable, therefore we can discard it as well, it will not provide us with any valuable information at low frequencies. After these considerations the resulting polynomial, with the low frequency poles is as follows: The numerators, b g xx (s), corresponding to each transfer function of G(s) = [g 11 (s) g 12 (s); g 21 (s) g 22 (s)] are shown: In the same way as for the poles, the zeros at high frequencies are undesired. In b g 22 (s), a high frequency zero appears in the terms multiplied by β and no pressure variable appears. This term is discarded resulting in: The Bode diagram of the system shown in (14), G(s), and the system without the high frequency poles and zeros,Ĝ(s), is shown in Fig. 4. At low frequencies, both systems have the same response and as frequency increases, both systems' responses deviate.
The new learning transfer matrix design will be L(s) = G −1 (s) + C(s), with the simplified system. Not incorporating the high frequencies in L(s) will affect our learning rate, however, it will have little impact as the low frequencies have the greatest influence on the convergence rate. Furthermore, we will obtain an equivalent lower-order system from which the inverse calculation will be simpler.
We include a low-pass filter, to attenuate the high frequencies of the system that are not included in the plant design. Furthermore, with this filter, we will attenuate the non-modeled physical dynamics of the proportional valve and the pump plate, as they are unknown.
Incorporating a fourth-order low-pass filter to the plant inverse modeling we obtain the following L(s) design: where ω c is the filter cutoff frequency in rad/s. ω c is set to 2 rad/s to ensure the robustness of the learning law to the high frequency dynamics that we have not modeled. However, producing a cutoff we introduce a considerable phase lag in our system, see the response ofĜ −1 f (s) in Fig. 5.
To avoid the additional phase lag, zero-phase low-pass filters (ZPF) are widely used in ILC design [22], [23], as they allow to filter the error signal forward and then, apply backward filtering to the filtered signal. Backward filtering generates phase lead to compensate for the phase loss of the forward filtering to achieve a zero-phase effect. The fourth-order filter is divided into two second-order low-pass filters, by which the iteration is filtered twice, once in the positive direction of time and once in the negative direction of time. The resulting plant inverse transfer matrix is: To perform the filtering the entire iteration error vector must be available for processing, therefore the ZPF cannot be done in real time, but one can perform the filtering between iterations. In Fig. 5, the Bode diagram of the resulting plant inverse with ZPF is shown,Ĝ −1 zpf (s), where the unknown high frequencies are attenuated without losing phase.
In order to show the advantages of not affecting C(s) with the low-pass filter, we compare the plant inverse design proposed in (18), L MIC (s), and the traditional direct model inverse design used in literature [33]- [36], L MI (s), with the inverse design without system simplifications nor low-pass filter introduction L D (s). The learning filters are defined as follows: The Bode diagrams of the three different model inverse approaches are shown in Fig. 6. It can be seen that L D (s) contains the high frequency dynamics of the slide hydraulic system, as no simplifications have been done and no low-pass filter has been applied.
We can see, in This is a consequence of applying a really low cutoff frequency in the low-pass filter introduced, as we want to get rid of most of the high frequencies and just kept those very low frequencies that we are interested in. In case we applied a higher cutoff frequency, the differences between L MIC (s) and L MI (s) would reduce at low frequencies, there would only be differences at the high frequency regions.

D. ZERO-PHASE FILTER WITH ANTI-WINDUP
The ZPF is applied to obtain the control inputs for the pump swash plate and the proportional valve position. Both actuators have a physical limitation, the input obtained from the ZPF must not exceed the range [-1, 0] and [0, 1] for the valve and pump, respectively. To limit the signal obtained from the ZPF, an anti-windup has been included to the filtering process.
The general recursive formulation for discrete-time filters in z-transform holds as follows: The forward and backward filters have been designed as second-order filters, therefore, we set M = J = 2. The resulting linear differential equation is: where n is a nonnegative integer. The static gain of the filters is given as: The filter output is given by (21). However, the obtained output value of the filter must not exceed the specified ranges, else the valve and pump control inputs will saturate.
If the output of the filter exceeds the ranges of interest, the anti-windup is applied to the recursive filter for the last J outputs and M inputs. In this way, we avoid the integration of the error for the recursive outputs and inputs, increasing the performance of the filter.
The resulting filter algorithm including the anti-windup is shown in Algorithm 1.

E. MIMO ILC STABILITY ANALYSIS
The convergence condition of the ILC problem has been widely considered in the literature to guarantee algorithm stability. The Hurwitz stability is used in [40], the Schur stability is applied in [41] and vertex Markov parameters are used in [42].
These methods require a large number of computations and cannot be applied directly to ILC. A more intuitive method is to analyze the stability graphically, from a frequency point of view, as it is done for SISO systems in [22], [24], [42], [43].

K s
Saturate previous inputs by multiplying the maximum value by the inverse of the static gain 17: end for 18: end if 19: end for 20: return: y The entire filtered signal is returned In these studies, the frequency response for all the frequencies up to the Nyquist frequency is plotted to analyze the convergence condition:

|1 − G(jw)S(jw)L(jw)| ≤ 1/Q(jw) ∀ω ∈ [−∞, ∞].(23)
Equation (23) can be interpreted on a Nyquist plot; if the frequency response remains inside the unit circle, the ILC algorithm will converge to zero error. Else, if the response gets out of the unit circle, the ILC algorithm will be unstable.
For the convergence condition of a MIMO ILC problem, the eigenvalues should be analyzed at all frequencies, instead of the frequency response as in the SISO case. We introduce again the MIMO convergence condition with the new L(s) design: So far, no studies have been done regarding the graphical analysis of the convergence condition of a MIMO ILC algorithm. We propose to apply the same procedure as in the SISO ILC, however, instead of looking to a Nyquist plot of the frequency response, we analyze the eigenvalues of (24) at all frequencies. VOLUME 9, 2021 From equation (24), if every eigenvalue of the term 2 (s) = I −S(s)G(s)(Ĝ −1 zpf (s)+C(s)) is less than one, then the convergence condition of the MIMO ILC will be guaranteed. 2 (s) is a two-by-two transfer matrix from which two eigenvalues can be obtained at each frequency point. Therefore, in order to guarantee the stability for all frequencies we evaluate the frequency response of the term 2 (s) up to the Nyquist frequency.
Ideally, if an accurate model inverse is included in the L(s) design, 2 (s) would be a zero matrix with every eigenvalue at zero for all frequencies. However, our L(s) design is far from being ideal as some simplifications have been considered in the inverse model and low-pass filters have been added to attenuate the uncertain high frequencies.
Hence, one should aim to obtain every eigenvalue as close as possible to the origin. The closer the eigenvalues to the origin, the faster those error frequencies will be corrected by the ILC algorithm. Figure 7 shows each eigenvalue obtained from 2 (s) at each frequency. We are not interested in the phase of the eigenvalues, but in their module, else the error would not asymptotically converge. It would in general oscillate around a value and eventually decay to it.
At low frequencies, the eigenvalues are close to zero, as it has been pointed out in Fig. 7 for a frequency of ω = 2 rad/s. At high frequencies the eigenvalues value increases, therefore, it will take more iterations to correct the high frequencies of the error. Although the eigenvalues increase, none of them get out of the unit circle, which guarantees the stability of the algorithm at high frequencies.
If the eigenvalues got out of the unit circle at a specific frequency, we could make use of Q(s) to enlarge the stability circle so the eigenvalues at all frequencies remain inside the new stability circle. By setting a |Q(s)| < |I| the stability circle increases accordingly, see in (24) that the smaller the value of Q(s) the bigger the value of the right-hand side. However, one should be cautious in the Q(s) design, as a Q(s) value different from one affects the convergence performance of the ILC algorithm, as it was shown in [44].
Some authors design Q(s) as a low-pass filter, to ensure the condition for stability is met [44], [45]. Using this Q(s) design approach, one can determine which frequencies are emphasized in the learning process. However, perfect tracking will not be achieved.

III. SIMULATION RESULTS
The designed MIMO control ILC algorithm has been implemented in a nonlinear model of a hydraulic press in Matlab/ Simulink. The press model has been developed in a novel library made by Ikerlan [46]. With this library, components are easily parametrized with data-sheet information in contrast to Simscape, which requires acquiring constructive parameters obtainable only from laboratory experiments. This novel approach allows reproducing the physical behavior of industrial components with high precision without losing Real-Time capabilities. The ILC algorithm is compared to the performance of a PI controller, which is the most prevalent controller applied in hydraulic circuits. Two different PI controllers have been designed, for the pump swash plate angle and the auxiliary chamber proportional valve spool position. The proportional valve PI controller gains have been set to K Pv = 3 and K Iv = 3. The pump PI controller gains have been set to K Pp = 3.86e-05 and K Ip = 0.82. Both PI controllers have been discretized applying Tustin approximation as the hydraulic press models run in discrete time with a fixed step-size of 0.002s.

A. POSITION CONTROL RESULTS
As explained in Sec. I, the control objective is to track a predefined slide trajectory by controlling the pump swash plate angle and the proportional valve spool position. Figure 8 shows both Free Fall and Drawing phases, where the position control is carried out. The ILC signal is introduced at t ≈ 0.5s, before the transition between the two phases.
At the first iteration, with the PI controller, there exists a manifest rebounding which results in the deceleration of the slide velocity. After the rebounding, the proportional valve spool position oscillates around -0.2, see Fig. 9, and it is not until t ≈ 3.5s that the position reference is reached.
As iterations go on, the ILC corrects the proportional valve spool position during both phases transition, introducing an overshoot before oscillating around -0.2 input signal value. This overshoot avoids the slide rebounding, resulting in a perfect tracking of the Free Fall and Drawing transition. After this transition, the input signal reaches the value -0.2 and then, oscillates smoothly around this value to maintain the constant falling velocity of the slide.
This effect can be seen in the slide velocity progress through iterations in Fig. 10. At the first iteration, there exists a notorious overshoot in the velocity signal which is attenuated as iterations go on. With the introduction of the ILC signal, once the slide reaches -100 mm/s speed it is maintained at this level. The pump swash plate angle normalized between [0-1], and the auxiliary chamber pressure are shown in Fig. 11 and Fig. 12, respectively. During the Free Fall phase the auxiliary chamber pressure is given by the slide weight and no pressure control is done. This can be seen in Fig. 11 where the pump plate is closed until t ≈ 1.0s.
Once the Drawing phase starts, from t ≈ 1.0s onwards, the pump swash plate is opened completely to reach the pressure reference of 180 bar. At the first iteration, with the PI controller, the pressure signal oscillates harshly around the reference as the swash plate angle total input has large oscillations. These swash plate oscillations are reduced as iterations go on, and the total input settles around 0.7 value. With this improvement the pressure reference signal is reached faster and the oscillations are considerably reduced.
It should be pointed out that perfect pressure tracking is not needed, it is sufficient to keep the auxiliary chamber pressure close to the reference so, when the slide strikes the workpiece during the Drawing phase, the extra pressure needed can be supplied by the auxiliary chamber pressure.
As a performance index for the proposed position control ILC algorithm, the RMSE between the desired trajectory and the slide trajectory over iterations is shown in Fig. 13. The RMSE is reduced by a factor of 7 in the position tracking, and a fast convergence is obtained as at the seventh iteration the error is considerably reduced, with respect to the first iteration with a PI controller.

IV. ILC PERFORMANCE COMPARISON
In this study, we propose an ILC algorithm design based on the simplified model inverse and the feedback controller (MIC-ILC). With this design, we optimize the convergence rate of the ILC algorithm to reduce the required number of iterations.
To analyze the performance of the MIC-ILC, we compare it with the proportional-type ILC (P-ILC) [47] and the traditional direct model-inverse ILC (MI-ILC) algorithm approach with direct inversion and the low-pass filter influencing the controller [34]. The comparison will be done for  the MIMO position control problem, and the RMSE between the position reference and the slide position will be used as the performance index.

A. P-ILC DESIGN
The P-ILC is one of the most used ILC algorithms in the literature due to its simplicity, as it calculates an input signal proportional to the error. In this section we will extend the P-ILC algorithm proposed by Arimoto in [47] to the MIMO case. The MIMO P-ILC learning control law is as follows: where K P is the proportional learning matrix, which is used to maximize the convergence rate of the algorithm. The following K P design is proposed: The eigenvalues of the P-ILC algorithm at each frequency are shown in Fig. 14, which we compare with the response obtained with MIC-ILC in Section II-E. With the P-ILC, at zero frequency, the eigenvalues lie at (+1,0), whereas with VOLUME 9, 2021 FIGURE 11. Total input to the pump swash plate and ILC input to the pump swash plate over iterations. MIC-ILC, at zero frequency, the eigenvalues lie at the origin. This is due to the cancellation of the model dynamics at low frequencies, that optimizes the convergence rate of the ILC algorithm.
The low frequencies will be corrected faster with MIC-ILC algorithm than with the P-ILC, however, we can still obtain fast convergence with the P-ILC if K P is increased so that the eigenvalues are as close as possible to the origin at low frequencies. However, the very low frequencies will always be close to (+1,0), penalizing its performance. As it is pointed out in Fig. 14, in both designs at 2 rad/s frequency the eigenvalues are close to the origin, remarkably closer with the MIC-ILC.
The RMSE between the position reference and the slide position with the P-ILC and the MIC-ILC algorithm is shown in Fig. 15. A similar convergence rate is obtained with both algorithms, although the P-ILC converges to a bigger steady-state error.
As it was shown in [19], the P-ILC can become unstable in the presence of measurement noise and state disturbance, and one would expect to obtain a noisy signal from the hydraulic press position sensor. Therefore, we add a varying white noise to the position sensor to analyze the P-ILC performance under noisy conditions. The noisy slide position error during the Free Fall and Drawing phases is shown in Fig. 16.
The RMSE over iterations is shown in Fig. 17, the P-ILC algorithm converges fast at the first iterations, however, as iterations go on it gets unstable and the error increases. We could decrease the gains in (26), so the amplification of the error is not that severe. However, decreasing the P-ILC gains the convergence rate is also decreased, see Fig. 17. The MIC-ILC algorithm remains stable due to the low-pass filters included in its design.

B. TRADITIONAL MODEL INVERSE DESIGN
In this section the traditional model inverse ILC (MI-ILC) will be designed, in which the low-pass filter is applied to the entire sensitivity function. The approach proposed by [34] is followed and the learning filter is designed as follows: In order to attenuate the model high frequency uncertainties and the noise added to the position sensor, a fourth-order low-pass filter is added to L MI (s) and zero-phase filtering is carried out to avoid the phase loss. Following [33] design, the filter is added to the entire L MI (s) transfer matrix: with ω c = 2 rad/s. Note that the placement of the low-pass filter changes with respect to the MIC-ILC algorithm, recall equation (18), where the low-pass filter is only applied to G −1 (s). This variation directly affects the convergence, as the model differences are bigger and the eigenvalues lie further from the origin, see Fig. 18.
In both MI-ILC and MIC-ILC designs, at zero frequency the eigenvalues lie at the origin as a result of the model dynamics cancellation. Therefore, the MI-ILC will perform  better than the P-ILC as the very low frequencies lie close to the origin. However, due to the application of the low-pass filter to the entire L MI (s), the eigenvalues grow apart from the origin faster in the MI-ILC than in the MIC-ILC algorithm. This can be seen in particular at 2 rad/s frequency point.
In Fig. 19, a faster convergence rate is obtained with the MIC-ILC in comparison to the MI-ILC. The former converges in 10 iterations and the latter in 20 iterations. Although the MI-ILC obtains fast convergence, from iteration 20 onwards the convergence rate gets reduced and the error converges gradually. This is a consequence of the eigenvalues distance from the origin as the frequency increases. With the MIC-ILC, the eigenvalues are closer to the origin so the low frequency components of the error are corrected faster.
Apart from the RMSE as a performance indicator, the Power Spectral Density (PSD) also provides useful insights for the comparison between ILC algorithms. In Fig. 20, the PSD of the error signal over iterations is shown. It can be seen that the energy is concentrated on the frequency band [0 ∼ +30] Hz.  The simulation has been carried out without noise, so the performance of the three algorithms can be analyzed in an ideal scenario. The main enhancement is at low frequencies, where the MIC-ILC outperforms the other ILC algorithms. At the first 10 iterations, the MIC-ILC reduces faster the low frequency components of the error. Especially the frequency band from 0 Hz to 5 Hz is decreased considerably, as those frequencies are where the eigenvalues lie close to the origin.
At higher frequencies, the three algorithms perform similarly as, at those frequencies, we only seek stability regardless of the convergence rate. It can be seen that at high frequencies the error reduction is slower, however, the MIC-ILC still achieves a faster correction rate. In line with the results obtained with the RMSE, at iteration 20, the three algorithms learning process has already converged and the difference between iterations is minimal, the frequency component of the error remains similar.

V. EXPERIMENTAL VALIDATION IN HYDRAULIC TEST RIG
The proposed MIC-ILC position control algorithm has been implemented in the hydraulic test rig shown in Fig 21. VOLUME 9, 2021 It consists of two identical hydraulic circuits and each circuit comprises of: a double-acting cylinder, a 4-way proportional valve, a hydraulic pump, a pressure relief valve, an accumulator, a pressure sensor at each cylinder chamber and a flow sensor.
The hydraulic pumps are fixed displacement pumps, which flow displacement varies depending on the rotational speed of the pump. With a fixed displacement pump, instead of controlling the swash angle, as in the pump shown in Section II-B, the pump's angular velocity is controlled. The relationship between the fixed displacement axial piston pump outflow and the shaft speed is as follows: where ω s is the shaft rational speed (rad/s). Every system actuator and sensor are connected via Ethernet communication to a Beckhoff 6930-0050. The Beckhoff Industrial PC (IPC) acts as a control cabinet and enables real-time communications and high-performance to control the hydraulic test rig. The Beckhoff IPC receives and processes the sensor signals, and delivers back the corresponding control inputs to the actuators. The position MIC-ILC algorithm is implemented in TwinCAT 3, which is a platform developed by Beckhoff, where PLC programming is carried out.
The position reference the press must follow during the Free Fall and Drawing phases is shown in Fig. 22. From the figure, the difference between the two phases falling velocity can be seen, this transition takes place at t ≈ 2s. For the pressure reference tracking, a 60 bar pressure level is specified, so it can be maintained during the Free Fall and Drawing phases.
As a first approach we design two PI controllers, to control the pump and the valve. The pump PI controller gains values are: K Pp = 4 and K I p = 3.8. The valve PI controller gains values are: K Pv = 1.5 and K I v = 0.25. The resulting position and force control is shown in Fig. 23 where the coupling between the two control loops can be seen. The pressure signal oscillates around the pressure reference of 60 bar. When the pressure is higher than the reference is due to an excessive pump velocity, which yields faster falling velocity of the slide. On the contrary, when the pressure is lower than the reference, is due to insufficient pump velocity, which translates into a slide velocity slowing down.
To improve the PI controller performance the MIC-ILC algorithm is implemented. The MIC-ILC input signals are introduced from the start of the Free Fall phase to the end of the Drawing phase.
We design the learning gain based on (18). The two PI controllers' gains are depicted above, and the system is modeled following (14) state space design with the test rig parameters depicted in Appendix C. The resulting eigenvalues of the frequency response matrix are shown in Fig. 24.
A similar response to that shown in Fig. 7 for the simulations is obtained. At low frequencies, the eigenvalues are close to zero, as it has been pointed out for a frequency of ω = 2 rad/s. At high frequencies the eigenvalues value increases, diverging from the origin. However, the responses get out of the unit circle, therefore we make use of the Q matrix to enlarge the stability circle. A Q = 0.93I is set, which will penalize the convergence of the algorithm but will keep it stable.
The result is shown in Fig. 25. The position tracking has already converged at iteration 10, where the oscillations existing with the PI controllers are reduced. In Fig. 26, how the MIC-ILC has adjusted the valve spool position to the slide velocity change can be seen. At iteration 10, from t ≈ 0s to t ≈ 2.2s, the valve spool position is reduced progressively from -0.2 to -0.1, to achieve fast velocity during the Free Fall phase. Once in the Drawing phase, the valve spool position is maintained at a -0.1 value, so a constant falling velocity is achieved.
The oscillations in the pressure reference tracking are also reduced, as shown in Fig. 27. At the beginning of the step the MIC-ILC introduces a positive velocity, see Fig. 28, to reach  the reference pressure as fast as possible, and then, a negative signal, to reduce the overshoot. From t ≈ 2s on, the MIC-ILC introduces a constant velocity signal to maintain the pressure level at the reference.
As a performance index for the proposed MIMO ILC algorithm, the RMSE between the position reference and the slide position signal is shown in Fig. 29. The RMSE is reduced by a factor of 4 in the position tracking, and a fast convergence is obtained as at the tenth iteration the error is VOLUME 9, 2021   considerably reduced, with respect to the first iteration with the PI controller.   The RMSE between the 60 bar pressure reference and auxiliary chamber pressure signal is shown in Fig. 30. The error is reduced by a factor of 2.6 regarding the first iteration with the PI controller.

VI. CONCLUSION
A new hydraulic press multiple-input-multiple-output (MIMO) position control based on Iterative Learning Control (ILC) is proposed. The MIMO ILC position control eliminates the need to manually define control signals and reduces the number of iterations required to converge with respect to other ILC algorithms.
The proposed ILC algorithm design is based on the simplification of the model inverse and the existing plant feedback controller (MIC-ILC) to increase the error convergence rate. A graphical evaluation is provided to analyze the stability and convergence condition. With this method, the stability of the algorithm is fulfilled if the learning filter eigenvalues lie inside the unit circle. The closer the eigenvalues lie to the origin, the faster the MIMO ILC algorithm will converge.
The MIMO ILC position control has been implemented in a nonlinear hydraulic press model in Matlab/Simulink for which satisfactory results have been obtained. The proposed MIC-ILC algorithm has been compared to the most common ILC algorithms in the literature. By means of the proposed graphical stability design, as the MIC-ILC eigenvalues lie closer to the origin than the other algorithms, faster convergence rate has been achieved. This is one of the main contributions of the MIC-ILC, the ability of correcting the low frequency components faster than the existing ILC algorithms in the literature. The enhancement of the MIC-ILC has been verified in the frequency domain, with the Power Spectral Density, analyzing the error signal at each iteration for each ILC algorithm. The evolution of the error signal in the frequency domain also shows that the MIC-ILC outperforms other ILC methods in terms of stability and convergence rate.
The MIC-ILC has been implemented in a hydraulic test rig available at Ikerlan Research Center. The implementation results correspond to the results obtained through simulations. The MIC-ILC not only achieves good tracking in the position and pressure signals but also gets fast convergence.
The proposed algorithm has shown to perform correctly in hydraulic press tests, in which the experiments are mainly carried out at low frequencies. It remains to analyze the performance of the algorithm at higher frequency response experiments, such as in fatigue testing machines or servo drives.
Although ILC yields satisfactory results and the convergence rate that we obtain is high with respect to other works, it requires starting the learning from scratch every time the reference is modified. This could be a drawback when implementing it in a real press operation, as every time the design of a workpiece is modified the ILC algorithm would need to perform the entire learning task to obtain the optimal input.
Techniques such as Reinforcement Learning (RL) could be of use, to obtain a nonlinear policy from a specific cost function. Some works have used PILCO probabilistic model [48], for learning a throttle valve control in a combustion engine with successful results [49]. It is to be analyzed whether via RL techniques the performance of a PI controller could be further improved.

A. LINEARIZED SYSTEM PARAMETERS
The complete expressions of the terms introduced in the linearized state-space model, introduced in (14), are as follows: The system parameters appearing in the terms introduced in the linearized state-space model, shown in (14), are as follows: • Cylinder moving mass: m = 26500 kg.