Marco Faroni
email: m.faroni003@unibs.it

Manuel Beschi
email: manuel.beschi@itia.cnr.it

Lorenzo Molinari Tosatti

Antonio Visioli
email: antonio.visioli@unibs.it

Lorenzo Molinari Tosatti
email: lorenzo.molinaritosatti@itia.cnr.it

A Predictive Approach to Redundancy Resolution for Robot Manipulators

Keywords: robot manipulators, redundant robots, continuous model predictive control, inverse kinematics, hard joint limits, optimal redundancy handling

In this paper, we propose a new method for the online redundancy resolution of robot manipulators, which implements a predictive strategy to calculate the optimal control action. In this way, it is possible to obtain a more efficient handling of the constraints, which represents one of the main issues in online resolution methods. The predictive model has been obtained by considering every joint as a kth-order integral system, and the predictive equations are derived from a continuous-time formulation. This allows the use of an irregular distribution of the prediction and control time instants and, as a consequence, longer prediction and control horizons can be obtained, without increasing the computational complexity of the algorithm. Finally, joint hard bounds are easily included in a linear-model-predictive-like framework, and the optimal control action is calculated by solving a linear quadratic problem. Simulation results for a 4-degree-of-freedom planar arm show the effectiveness of the method compared to purely local resolution techniques.

INTRODUCTION

Redundant manipulators have been broadly exploited in industry and research, since they allow greater versatility in the task execution. As a matter of fact, the extra degrees of freedom permit to satisfy (completely or partially) secondary tasks, or to optimise a given objective function such as manipulability maximisation, joint availability or minimum energy consumption. Several methods for the optimal resolution of the redundancy has been proposed during the years. They are classified mainly in local and global methods [START_REF] Nakamura | Advanced Robotics: Redundancy and Optimization[END_REF]: the former class solves the redundancy using information related to the actual configuration of the robot, whereas the latter aims at obtaining a globally optimal solution, taking into account the whole desired task in the optimisation. Of course, local methods generally lead to poorer global performance [START_REF] Faroni | A global approach to manipulability optimisation for a dual-arm manipulator[END_REF], but they allow an online resolution thanks to their lighter computational burden. On the contrary, the computational complexity of global methods is too high for real-time control applications and, therefore, they need to be performed offline. As regards local methods, they are based on the resolution of the inverse kinematics of the robot, in order to obtain joint velocities and/or accelerations. The typical approach is based on the pseudo-inverse of the Jacobian of the manipulator [START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF]. One of the main issues in the redundancy handling is the satisfaction of joint position, velocity and acceleration bounds. Pseudo-inverse methods do not explicitly take them into account. Typically, hard position bounds are converted into soft ones by means of a cost function that keeps the joints in the middle of their ranges [START_REF] Liegeois | Automatic supervisory control of the configuration and behavior of multibody mechanisms[END_REF]. However, the impossibility of handling velocity and acceleration bounds would generate the saturation of joint commands, leading to the deformation of the primary task. In order to consider also these constraints, the local redundancy problem may be considered as a Quadratic Program (QP), with linear equality and inequality constraints [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF]. In this way, joint constraints are considered explicitly and an optimal solution is found in case the task is locally feasible. However, the solution of these method is only locally optimal, and it does not take into account the future evolution of the task, the constraints, and the kinematic states. This means that the iterative local resolution of redundancy might move the robot to a disadvantageous configuration with respect to the constraints. This would cause undesired velocity/acceleration peaks or even the impossibility of satisfying the constraints with bounded control actions. In case no solution that satisfies all the constraints exists, the saturation of the actuators leads to the deformation of the task. Recent methods, for instance, implement task scaling techniques to preserve at least the geometrical task when no local feasible solution exists [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF]. Anyway, considering that constraints and bounds along the task are usually known a priori, those situations could often be avoided, and this represents one of the main drawbacks of purely local approaches, as pointed out in [START_REF] Schuetz | Predictive online inverse kinematics for redundant manipulators[END_REF]. This paper proposes a method to solve the inverse kinematic of a redundant manipulator by means of a model predictive approach. One great advantage of MPC techniques is the efficient handling of constraints. Moreover, the predictive approach permits to take into account the constraints in the successive time instants as well. For this reason, such approach seems to be a good candidate for overcoming the problems of punctual resolution of the redundancy in a robotic manipulator. Following this principle, [START_REF] Schuetz | Predictive online inverse kinematics for redundant manipulators[END_REF] exploited Pontryagin's Minimum Principle to set up a nonlinear optimisation problem over a finite horizon. In [START_REF] Zube | Cartesian nonlinear model predictive control of redundant manipulators considering obstacles[END_REF], by contrast, a typical nonlinear MPC (NMPC) is applied to solve the inverse kinematics of a redundant manipulator. Although the results are certainly satisfactory in terms of kinematic performance, real-time applications for longenough horizons seem to be still far away. In this work, we propose a different approach, where every joint of the robot is treated as an kth-order integral system, while the desired trajectory is considered as a nonlinear constraint in the optimisation. As shown in Section 3.2, this constraint can be easily linearised using the prediction of the state variables. In this way it is possible to reduce the optimisation of some typical cost functions (e.g. minimum control effort and joint availability) to a linear quadratic programming problem, which can be solved much more rapidly. In order to obtain control and prediction horizons sufficiently long for the purposes of the paper, control and prediction time instants are not equally distributed along the horizons. As a matter of fact, the time instants are more frequent at the beginning of the horizons and less frequent in the end. This strategy permits to reduce the number of variables involved in the prediction and in the optimisation and, therefore, it dramatically decreases the computational burden of the algorithm. Another advantage of the proposed method is that first and second (or even third) order resolution levels may be obtained just by using either a single or a double (or triple) integral system respectively, without any modification to the algorithm. Moreover, the intrinsic structure of linear MPC permits to handle position, velocity and acceleration boundaries in an elegant and efficient manner.

BACKGROUND

Consider an n-degree-of-freedom manipulator and an mdimensional task with m < n. The associated direct kinematics is given by: x = f (q) (1) where x ∈ R m is the task position vector and q ∈ R n is the joint position vector. Differentiating with respect to time gives: ẋ = J q (2) where J = ∂f /∂q is the m × n Jacobian matrix, ẋ ∈ R m is the task velocity vector and q ∈ R n is the joint velocity vector. Typical inverse kinematics are based on the solution of the linear system (2) by means of the pseudoinverse of the Jacobian. This can be performed at velocity or acceleration level [START_REF] De Luca | Issues in acceleration resolution of robot redundancy[END_REF]. However, pseudoinverse methods do not provide an intrinsic way of handling boundaries and constraints. As a matter of fact, their mere application may cause the violation of position, velocity and acceleration boundaries if no strategies are implemented to avoid it. Moreover, the violation of velocity and acceleration bounds would saturate the actuators, with consequent deformation of the primary task. Alternatively, it is possible to represent the redundancy problem as a quadratic optimisation, as follows:

q = argmin 1 2 e( q) T e( q) s.t. J(q) q = ẋ ; Qmin ≤ q ≤ Qmax (3)
where Qmin and Qmax represent the kinematic constraints at velocity level.

When the function e( q) is linear with respect to q, the minimisation becomes a QP.

The extension of such method to the acceleration level is easily derived and gives: q = argmin 1 2 e(q) T e(q)

s.t. J(q)q = ẍ -J q ; Qmin ≤ q ≤ Qmax (4)
where Qmin and Qmax represent the kinematic constraints at acceleration level.

The kinematic constraints typically include position, velocity and/or acceleration limits, as extensively explained in [START_REF] Flacco | Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF].

PROPOSED METHOD

Formulation of the predictive model

Consider an n-degree-of-freedom manipulator and an mdimensional task, wit m < n. From a control perspective, the aim of redundancy resolution consists in finding the optimal control actions that have to be applied to the joints, taking into account the constraints introduced by the Cartesian task. To this end, every joint is modelled as a kth-order integral system:

q i (s) = 1 s k u i (s) ∀ i = 1 . . . n (5)
Equivalently, choosing the state vector as

ξ i := [q i qi . . . q (k-1) i ] T ∀ i = 1 . . . n (6)
the following state-space formulation ( Ã, B, C) results:

ξi (t) = Ãξ i (t) + Bu i (t) (7) q i (t) = Cξ i (t) (8)
where the matrices à ∈ R k×k , B ∈ R k×1 , and C ∈ R 1×k are given by:

à = 0 1 0 0 0 1 0 0                 (9) B = [ 0 • • • 0 1 ] T (10) C = [ 1 0 • • • 0 ]
(11) (Note that in the remarkable case of k = 1, the matrices reduce to the scalar case, in which à = 0, B = 1, C = 1.)

In this way the control action u i results to be the kth derivative of the corresponding joint position. Choosing k = 1, 2, or 3 the joint will be controlled at velocity, acceleration, or jerk level respectively. It is now possible to obtain the kinematic model of the whole manipulator from the single-joint models by coupling them in state-space formulation (A, B, C), where:

A := blkdiag Ã, . . . , Ã ∈ R nk×nk (12) B := blkdiag B, . . . , B ∈ R nk×n (13)
C := blkdiag C, . . . , C ∈ R n×nk (14) where blkdiag(•) generates a block-diagonal matrix from the given matrices. The complete state vector results:

ξ = [q 1 , q1 , . . . q (k-1) 1 , . . . , q n , qn , . . . q (k-1) n ] T
(15) whereas the output vector is:

q = [q 1 , q 2 , . . . q n ] T (16)
Once the state-space model of the whole manipulator has been obtained, the application of classic discrete-time MPC would result straightforward. On the other hand, in this work, a continuous-time approach has been preferred. Such choice allows the irregular distribution of prediction and control time instants along the horizons, as will be explained later. Therefore, consider a time instant t 0 and the corresponding state ξ 0 = ξ(t 0 ). The predicted output at a generic time t = t 0 + τ , τ > 0, is given by:

q(t 0 + τ ) = C e Aτ ξ 0 + C e Aτ τ 0 e -Aγ Bu(γ) dγ (17)
Now assume that B is constant and u is a staircase function. Let ti be the time instants, with respect to t 0 , in which the steps of the function are applied and let u(t 0 + ti ) be the input value in the time interval between ti-1 and ti . These intervals are defined, then, as ∆ t = [ t1 -t 0 , t2 -t1 , , . . . , tctc-1 ], c > 0. Note that the lengths of the intervals does not need to be the same. Thus, (17) becomes:

q(t 0 + τ ) = C e Aτ ξ 0 + C c i=1 Γ i u(t 0 + ti ) (18)
where: 18) with respect to τ gives:

Γ i :=      e A(τ -ti) ∆ ti 0 e Aγ dγ B, if τ ≥ ti τ -ti-1 0 e Aγ dγ B, if ti-1 < τ < ti 0, if τ < ti ∧ τ < ti-1 (19) assuming t0 = 0. Differentiating (
q(t 0 + τ ) = CA e Aτ ξ 0 + C c i=1 Γi u(t 0 + ti ) ( 20
)
where:

Γi =      A e A(τ -ti) ∆ ti 0 e Aγ dγ B, if τ ≥ ti e A(τ -ti-1) B, if ti-1 < τ < ti 0, if τ < ti ∧ τ < ti-1 (21)
Repeating for a dth-order derivative, the following equation results:

q (d) (t 0 + τ ) = CA d e Aτ ξ 0 + C c i=1 Γ (d) i u(t 0 + ti ) (22)
where d = 1, . . . , k -1, and

Γ (d) i =    A d Γ, if τ ≥ ti A d-1 e A(τ -ti-1) B, if ti-1 < τ < ti 0, if τ < ti ∧ τ < ti-1 (23) assuming A 0 = I in case A = 0.
Having defined the control time instants t, a series of constant matrices Φ d,τ can be defined for a given prediction time τ , such that:

Φ d,τ := C Γ (d) 1 , C Γ (d) 2 , . . . , C Γ (d) c ( 24
)
for d = 0, . . . , k -1, having assumed Γ 24) in ( 18) and ( 22) gives respectively:

(d) i = Γ i if d = 0. Substituting (
q(t 0 + τ ) = C e Aτ ξ 0 + Φ 0,τ u ( t) (25)
q (d) (t 0 + τ ) = CA d e Aτ ξ 0 + Φ d,τ u ( t) (26)
where u ( t) := [u(t 0 + t1 ), . . . , u(t 0 + tc )] T .

Extending this reasoning to a generic number of prediction times [τ 1 , . . . , τ p ], p > 0, leads to the following systems of equations:

       q(t 0 + τ 1 ) = C e Aτ1 ξ 0 + Φ 0,τ1 u ( t)
. . .

q(t 0 + τ p ) = C e Aτp ξ 0 + Φ 0,τp u ( t) (27)
       q (d) (t 0 + τ 1 ) = CA d e Aτ1 ξ 0 + Φ d,τ1 u ( t) . . . q (d) (t 0 + τ p ) = CA d e Aτp ξ 0 + Φ d,τp u ( t) (28)
As for the control time instants, the prediction instants [τ 1 , . . . , τ p ] are not necessarily equally spaced.

The systems ( 27) and ( 28) can be rewritten in a single matrix form as:

q(d) = L d ξ 0 + F d u ( t) d = 0, . . . , k -1 ( 29
)
where q ∈ R np×1 is the predicted output vector at times (t 0 + τ 1 , . . . , t 0 + τ p ), u ( t) ∈ R nc×1 is the vector of the control inputs, and:

L d =   
CA d e Aτ1 . . .

CA d e Aτp    ∈ R np×nk , F d =    Φ d,τ1 . . . Φ d,τp    ∈ R np×nc .
The formulation presented in ( 29) is convenient since the output and its derivatives have clear physical meaning.

In this way, the series of matrices L d and F d permits to easily express joint positions, velocities or accelerations (depending on the order k of the model) with respect to the control action. For instance, choosing k = 2, the joint positions and velocities are given by: q

= L 0 ξ 0 + F 0 u ( t) (30) q = L 1 ξ 0 + F 1 u ( t)
(31) whereas the vector u ( t) represents the future accelerations at time instants (t 0 + t1 , . . . , t 0 + tc ).

In classic MPC, prediction and control instants are equally distributed along the respective horizons and correspond to the sampling times. However, robotic systems generally have sampling periods of the order of few milliseconds, not allowing long-enough horizons and low computational burdens at the same time. Although the derivation of the predictive equations in continuous time might seems less intuitive, their great advantage consists in the possibility of using a free distribution of prediction and control instants along their horizons, disregarding the real sample time of the system. As a matter of fact, it is not necessary for the time instants ti and τ i to be neither equally distributed nor coincident with a sample instant. Adopting an irregular partition of the predictive and control horizons, thus, it is possible to obtain much longer horizons, without increasing the computational burden of the algorithm.

Linearisation of the task constraint

Up to this point, a predictive model of the robot in the joint space has been obtained, with longer horizons thanks to the continuous-time formulation of the model. On the other hand, nothing has been said about the satisfaction of the Cartesian task. In order to perform the task, the robot has to satisfy (2) for all time instants. Imposing this constraint to every control instant ti of the control horizon gives:

J * q(t 0 + ti ) q(t 0 + ti ) = ẋ(t 0 + ti ) ∀ i = 1, . . . , c (32) where:

J * q(t 0 + ti ) = blkdiag J q(t 0 + t1 ) , . . . , J q(t 0 + tc ) .

This constraint is nonlinear with respect to the variables q(t 0 + ti ), due to the Jacobian. However, using (30) it is possible to calculate the prediction of the joint position. Therefore, approximating q(t 0 + ti ) with its prediction at the previous cycle, (32) becomes linear and depends only on the joint velocities q(t 0 + ti ). Moreover, considering that the relationship between the future joint velocities and the control action is given by (31), the task constraint becomes linear and dependent on the sole control action u ( t) :

J * (q)F 1c u ( t) = ẋ(t 0 + ti ) -J * (q)L 1c ξ 0 ∀ i = 1, . . . , c (33)
where F 1c and L 1c are the velocity prediction matrices as defined in ( 29), but related only to the control time instants. In other words, they correspond to the first nc rows of the matrices F 1 and L 1 used in (31). (Note that, in case of first-order model, the joint velocity is given simply by the control action u ( t) . So, in that case, the matrix F 1c would result equal to the identity matrix I c , while the matrix L 1c would be null.) Obviously, the constraint must be updated at each iteration with the new prediction q in the Jacobians, the new initial state ξ 0 , and the Cartesian velocities ẋ related to the corresponding time instants.

Handling of constraints

The formulation of the predictive model ( 29) allows a straightforward inclusion of constraints on linear combinations of the state variables. As a matter of fact, the use of ( 29) as prediction of the kinematic variables permits to express such constraints with respect to the control variables in a easy and intuitive manner.

In the resolution of the inverse kinematics, the main constraints are represented by the mechanical limits of the joints and the boundaries of the kinematic variables, such as velocity and acceleration. Therefore, chosen the order k of the kinematic model, these constraints are given by: q min ≤ q ≤ q max . . .

q (k) min ≤ q (k) ≤ q (k) max (34)
By using ( 29) it is possible to impose these boundaries with respect to the control actions u ( t) as:

q min ≤ L 0 ξ 0 + F 0 u ( t) ≤ q max . . . q (k) min ≤ u ( t) ≤ q (k) max ( 35
)
which can be easily rewritten in the typical form of linear inequality constraints with respect to u ( t) as:

F d u ( t) ≤ q (d) max -L d ξ 0 ∀ d = 0, . . . , k -1 -F d u ( t) ≤ -q (d) min + L d ξ 0 ∀ d = 0, . . . , k -1 u ( t) ≤ q (k) max -u ( t) ≤ -q (k) min (36)

Optimal control action

In classic linear MPC, the control action to apply to the system is obtained by minimising a cost function in the following form:

η(u) = τp 0 e 2 (t)dt + λ tc 0 u 2 (t)dt ( 37
)
where e is a linear function of the output and, due to the linearity of the model, of the control variable as well. This cost function may be approximated as follows:

η(u) ≈ e(τ ) T ∆τ e(τ ) + λu T ( t) ∆t u ( t) (38)
where ∆τ = diag(τ 1 , τ 2 -τ 1 , . . . , τ p -τ p-1 ), ∆t = diag(∆t 1 , . . . , ∆t c ), e(τ ) is given by the function e evaluated at time instants (t 0 + τ 1 , . . . , t 0 + τ p ), and λ > 0.

Typically, e is given by the difference between the output variable and a reference signal. In this way, the optimisation of the cost function leads to a trade-off between the minimisation of the control error (first term in ( 38)) and the control effort (second term in ( 38)). Moreover, if the control variable u has to satisfy only linear constraints, the optimisation of ( 38) can be carried out by means of typical QP procedures. Applying the same approach to the inverse kinematics, it is possible to minimise a generic quadratic cost function with respect to the control variable u or to the state variables ξ (reminding the linear relationship between them, given by ( 29)).

The Cartesian task is satisfied by imposing the linearised constraint (33). Moreover, (33) requires the task constraint to be satisfied at the future control times as well, and this permits to obtain a prediction of the states which takes into account the task constraint within the control horizon. Finally, the inequality constraints (36) ensure that the kinematic variables remain between their boundaries within the control horizon. In this way, a QP with linear constraints results, that intrinsically takes into account the Cartesian task and the limits of the kinematic variables in an efficient and elegant manner.

Typical quadratic cost function in redundancy resolution

Typical quadratic functions used in the redundancy resolution are:

• Minimum control effort In this case only the second term appears in the cost function (37), then e(t) = 0 and

η(u) = τp 0 u 2 (t)dt ≈ u T ( t) ∆t u ( t) . (39)
Generally, if the order of the system is greater than 1, it is convenient to include in the cost function the minimisation of all the kinematic variables. The resulting cost function will be therefore:

η(u) = τp 0 k d=1 q (k) (t) 2 dt ≈ k-1 d=1 (L d ξ 0 + F d u) T ∆τ (L d ξ 0 + F d u) + u T ( t) ∆t u ( t) (40) • Joint-range availability
The distance between the joint positions and a given reference value q (usually the middle value of the joint range) has to be minimised, thus the first term e in ( 37) and ( 38) becomes: 41) where W is an appropriate weighting matrix.

e(t) = W (q(t) -q) ≈ W (L 0 ξ 0 + F 0 u ( t) -q) (

SIMULATION RESULTS

As an illustrative example, consider a 4-degree-of-freedom planar arm with unitary length axis. As shown in Figure 1, the base of the robot corresponds to the origin of the Cartesian plane and the task consists of a circle of centre C(0, 2) and radius r = 1 m. The end-effector has to execute the circle, starting from and returning to the point (1, 2) in a total execution time equal to 3 s. The scalar Cartesian velocity of the end-effector is requested to follow a trapezoidal timing law. In this example, we choose to minimise the cost function (41). In this way, the control action will be calculated in order to attract the joint positions in the middle of their ranges. At the same time, the control actions ensure that the joint variables respect boundaries and constraints along the control horizon. Thus, at each cycle, the optimal solution corresponding to the first time instant is applied to the system, following the receding horizon principle. This simulation aims at demonstrating that, using the proposed method, the resulting kinematic profiles give better global results in terms of satisfaction of the boundaries and minimisation of the control action. Similar conclusions can be obtained by minimising other performance indices (e.g. ( 39)). The results will be omitted for the sake of brevity. The method has been tested at both first and second kinematic levels, that is, k equal to 1 and 2, with the tuning parameters shown in Table 1. The number and the distribution of the prediction and control instants have been chosen empirically. As a general approach, the first time instants are chosen at a distance equal to the sampling time of the system, while the following ones are more and more sparse. The proposed method has been compared against the local quadratic optimisations (3) and ( 4), setting the joint range limits as reported in Table 2. The resulting joint kinematic states are shown in Figure 2 for the velocity level resolution, and in Figures 3 and4 for the acceleration level.

Although the QP solver considers the bounds on position, velocity and acceleration of the joints, the local optimisation causes abrupt changes in the control variable, since it is not able to foresee the activation of the constraints (Figures 2 and4, magenta lines). This situation often leads to the impossibility of satisfying both the primary task and the joint bounds, and task scaling techniques would be needed to preserve at least the geometric path. The proposed method, on the contrary, considers the presence of the boundaries along the whole control horizon. In this way, the algorithm calculates a smoother control action that tends to avoid its saturation. As a consequence, if the task is feasible for the robot capabilities, it is more likely to find a solution that does not cause task deformation.

Of course the better performance of the proposed method implies a heavier computational burden. Nonetheless, the calculation of the control action is still a QP problem and this permits to obtain a much lighter optimisation problem if compared to typical NMPC schemes, as [START_REF] Zube | Cartesian nonlinear model predictive control of redundant manipulators considering obstacles[END_REF]. The optimisation has been performed in Matlab by using the QP solver described in [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF].

The optimisation takes about 0.45 ms on average (using a 2.5 GHz Intel Core i5-2520M processor). During the

Time [s] 0 0.5 1 1.5 2 2.5 3 q1 [rad/s] -4 -2 0 2 Time [s] 0 0.5 1 1.5 2 2.5 3 q2 [rad/s] -3 -2 -1 0 1 Time [s] 0 0.5 1 1.5 2 2.5 3 q3 [rad/s] -2 -1 0 1 2 Time [s] 0 0.5 1 1.5 2 2.5 3 q4 [rad/s] -2 -1 0 1 2 Fig. 2.
Comparison between the joint velocities given by the local-QP (dashed magenta) and the proposed (solid blue) methods in the constrained case (k = 1).

Time [s] 0 0.5 1 1.5 2 2.5 3 q1 [rad/s] -4 -2 0 2 Time [s] 0 0.5 1 1.5 2 2.5 3 q2 [rad/s] -4 -2 0 2 Time [s] 0 0.5 1 1.5 2 2.5 3 q3 [rad/s] -1 0 1 2 Time [s] 0 0.5 1 1.5 2 2.5 3 q4 [rad/s] -2 0 2 4
Fig. 3. Comparison between the joint velocities given by the local-QP (dashed magenta) and the proposed (solid blue) methods in the constrained case (k = 2). Note the different axis scales. first cycle the optimisation takes about 3.8 ms, since the solver has to initialise the problem. Existing algorithms using NMPC are, at least, one order of magnitude slower, even on a dedicated hardware and with a C/Fortran implementation.

CONCLUSIONS

This paper presented a new method for the online redundancy handling of robot manipulators. It has been shown that the use of a model-based predictive strategy permits to take into account the future evolution of the kinematic variables of the robot. In this way, the control actions are optimised with respect to a desired cost function (e.g. joint range availability) and to the presence of joint boundaries, along a certain time horizon. Considering the presence of the boundaries in the future permits to take action before their activation and this leads to smoother trends of the kinematic states. Moreover, if the desired task is feasible for the manipulator, the algorithm is more likely to find a solution without deforming the original trajectory. The irregular distribution of control and prediction instants permits to have longer horizons without increasing the computational burden. Finally, the reduction of the Cartesian task to a linear constraint allows the use of linear MPC techniques, ensuring a computational complexity lower than the classic NMPC strategies used in the literature to the same purpose.

Fig. 1 .

 1 Fig. 1. Desired task and initial configuration of the robot.

Fig. 4 .

 4 Fig. 4. Comparison between the joint accelerations given by the local-QP (dashed magenta) and the proposed (solid blue) methods in the constrained case (k = 2). Note the different axis scales. first cycle the optimisation takes about 3.8 ms, since the solver has to initialise the problem. Existing algorithms using NMPC are, at least, one order of magnitude slower, even on a dedicated hardware and with a C/Fortran implementation.

Table 1 .

 1 Tuning parameters used in simulation.

	T [ms]	1
	Np	40
	Nc	10
	τ [ms] T, 2T, ..., 5T, 10T, ..., 30T, 40T, ..., 80T, 105T, ..., 705T
	t [ms]	0, T, 2T, ..., 5T, 10T, ..., 25T
	λ	5 • 10 -3

Table 2 .

 2 Position, velocity and accelerations boundaries for the simulated robot.

	Joint	q min qmax qmin qmax qmin [rad] [rad] [rad/s] [rad/s] [rad/s 2 ] [rad/s 2 ] qmax
	1	0	π	-4	4	-100	100
	2 -π/2 π/2	-4	4	-100	100
	3 -π/2 π/2	-4	4	-100	100
	4 -π/2 π/2	-4	4	-100	100

ACKNOWLEDGEMENTS

The research leading to these results has received funding from the European Union H2020 program under grant agreement n. 637095 (FourByThree) and ECSEL-2016-1 under grant agreement n. 737453 (I-MECH).


Marco Faroni
email: m.faroni003@unibs.it

Manuel Beschi
email: manuel.beschi@itia.cnr.it

Manuel Berenguel

Antonio Visioli
email: antonio.visioli@unibs.it

Fast MPC with Staircase Parametrization of the Inputs: Continuous Input Blocking

Fast MPC with Staircase Parametrization of the Inputs: Continuous Input Blocking

Marco Faroni, Manuel Beschi, Manuel Berenguel, Antonio

I. INTRODUCTION

Although microprocessors are every day more and more powerful, the reduction of computational complexity of Model Predictive Control (MPC) algorithms is still an open topic for both the research and the industrial community. Quoting Y. Wang and S. Boyd, "a widely recognized shortcoming of model predictive control is that it can be used in applications with slow dynamics, where the sample time is measured in seconds or minutes" [START_REF] Wang | Fast model predictive control using online optimization[END_REF]. This is due to the high computational burden required to calculate the optimal control action. As a matter of fact, the optimal control inputs are computed by solving a Quadratic Program (QP). Then, in the spirit of receding horizon control, the first element of the input sequence is applied to the system, and the procedure is iterated at each time-step [START_REF] Camacho | Model Predictive Control[END_REF], [START_REF] Maciejowski | Predictive control with constraints[END_REF]. The computational complexity of such algorithm may be prohibitive for small sampling periods and this limits the use of MPC in the control of fast-dynamic systems.

Several methods have been proposed to speed up the resolution of the online optimization [START_REF] Rossiter | Efficient algorithms for trading off feasibility and performance in predictive control[END_REF]. Many of them focus on improving the efficiency of the QP solver, by exploiting the peculiar structure of the problem [START_REF] Ferreau | An online active set strategy to overcome the limitations of explicit MPC[END_REF], [START_REF] Bartlett | Quadratic programming algorithms for large-scale model predictive control[END_REF]. Others approximate the constrained problem in order to quickly find a suboptimal solution [START_REF] Kouvaritakis | Who needs QP for linear MPC anyway?[END_REF], [START_REF] Kouvaritakis | Efficient robust predictive control[END_REF]. In [START_REF] Wang | Fast model predictive control using online optimization[END_REF] it is shown that a (very) approximated solution of the QP still leads to satisfactory control performance, alleviating the required computational time. An alternative approach for complexity reduction is the well-known Explicit MPC [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF], in which the state space is explored offline, in order to obtain an explicit solution for the control action. The online computation of the input, then, is reduced to the evaluation of a piecewise linear function. Another interesting technique is represented by the Eventbased MPC, which permits to reduce the occurrences of the computation of the control action by including the predictive algorithm in an event-based framework [START_REF] Pawłowski | Event-based gpc for multivariable processes: A practical approach with sensor deadband[END_REF].

However, the most intuitive and widespread approach for reducing the computational complexity consists in the reduction of the number of degrees of freedom of the controller, that is, the number of optimization variables of the QP problem.

In the classic formulation of MPC, the degrees of freedom of the optimization are equal to the number of inputs multiplied by the lengths of the control horizon (in terms of number of sampling periods). A common practice to reduce the number of variables consists in choosing a control horizon shorter than the prediction one [START_REF] Garriga | Model predictive control tuning methods: A review[END_REF]. Another widespread technique for reducing the number of variables involved in the optimization is the parametrization of the input functions [START_REF] Cagienard | Move blocking strategies in receding horizon control[END_REF]- [START_REF] Qin | A survey of industrial model predictive control technology[END_REF]. This is the case, for instance, of the so-called Predictive Functional Control [START_REF] Richalet | Predictive functional control: principles and industrial applications[END_REF]. The parametrization consists in approximating the control sequence by means of a function with a minor number of degrees of freedom. This is often implemented by approximating every input with a piecewise constant function (with less steps than the control horizon). In other words, the input sequence is forced to be constant along several timesteps in the resolution of the QP problem. This technique is usually referred to as Input Blocking (IB).

However, such parametrization is strongly limited by the sampling period of the system. As a matter of fact, the steps of the staircase parametrization must occur at a multiple of the sampling period. Considering that common prediction horizons are typically in the order of ten sampling periods, such constraint limits the goodness of the approximation. This could be overcome by formulating the predictive equations in the continuous-time domain.

Continuous-time MPC was extensively studied in several works [START_REF] Demircioǧlu | Continuous-time generalized predictive control (CGPC)[END_REF], [START_REF] Wang | Continuous time model predictive control design using orthonormal functions[END_REF]. Generally, in its formulation, the inputs are parametrized by means of a sum of exponential function, allowing simpler QPs compared to the discrete case. Of course, the resolution of the online optimization problem needs a finite time to be performed. During the time required by the optimization, the open-loop trajectory is implemented and an intermittent feedback approach can be used, giving rise to the so-called Intermittent Predictive Control, which somehow lies between the continuous-time and the discrete-time approaches [START_REF] Gawthrop | Intermittent predictive control of an inverted pendulum[END_REF], [START_REF]Constrained intermittent model predictive control[END_REF].

In this paper, we exploit the advantages of the continuoustime approach to formulate a predictive model in which the control inputs are parametrized as staircase functions (similarly to IB). This permits to set the control and prediction timeinstants arbitrarily, disregarding the sampling period of the system. In this way, the steps of the staircase parametrization do not need to coincide with a multiple of the sampling period, as happens with classic IB techniques. The same reasoning can be applied to the prediction time-instants: in classic MPC the output prediction is evaluated in the multiples of the sampling period, while our formulation allows the arbitrary choice of the prediction instants. With respect to IB, this allows a more flexible choice in the parameters, ensuring better flexibility in the controller. As a matter of fact, the piecewise constant function better approximates the input sequence, while keeping the same structure of the IB parametrization (staircase function and same number of variables). The basic principle of such method was successfully implemented in [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF], in order to cope with the fast sampling period of robotic systems. The method is also formulated in terms of relative increment of the control action (i.e. Δ𝑢), by imposing that such increment must be piecewise constant along the control horizon.

Note that even though the predictive equations are derived in the continuous-time domain, the resulting controller can be directly implemented in discrete-time, since it is based on the assumption that the control input is a piecewise constant function.

We also provide a tuning method in which we choose the optimal distribution of the prediction and control instants along the horizon, by minimizing the deviation, in terms of integral of the absolute error with respect to the full-degree-of-freedom (FDOF) controller, that is, the standard MPC without any parametrization. By applying the same optimal tuning rule to the IB technique, we show that our method guarantees a better approximation of the FDOF controller. In other words, given the same number of variables for both the staircase parametrizations and aiming at minimizing the deviation with respect to the FDOF controller, our method ensures better results than the ones of classic IB.

The paper is organized as follows. Section II recalls the basic notions of standard MPC and Input Blocking. In Section III the proposed method is presented; then Section IV discusses the tuning and the choice of the prediction and control timeinstants. In Section V a second-order fast-dynamics system is analyzed as a case study, in order to demonstrate the effectiveness of the proposed strategy. Finally, conclusions are exposed in Section VI.

II. PRELIMINARIES

A. Generalities

Consider a discrete-time linear dynamic model, with 𝑛 states, 𝑚 inputs and 𝑞 outputs, such that:

𝑥(𝑘 + 1) = 𝐴 𝑑 𝑥(𝑘) + 𝐵 𝑑 𝑢(𝑘) 𝑦(𝑘) = 𝐶 𝑑 𝑥(𝑘) (1)
where 𝐴 𝑑 ∈ ℝ 𝑛×𝑛 , 𝐵 𝑑 ∈ ℝ 𝑛×𝑚 , and 𝐶 𝑑 ∈ ℝ 𝑞×𝑛 . Now, consider the system at sampling time 𝑘 > 0 and the corresponding state 𝑥(𝑘). The prediction of the output can be calculated by forward solving the difference equation (1) from the current state 𝑥(𝑘). In particular, consider:

• a prediction horizon 𝑁 𝑝 , which represents the number of future sampling times in which the output prediction has to be calculated; • a control horizon 𝑁 𝑐 , which represents the number of future sampling times for which the control actions are taken into account (𝑁 𝑐 ≤ 𝑁 𝑝 ). The output prediction in the future 𝑁 𝑝 sampling times, calculated at time 𝑖, is given by:

𝑦(𝑘 + 1|𝑘) = 𝐶 𝑑 𝐴 𝑑 𝑥(𝑘) + 𝐶 𝑑 𝐵 𝑑 𝑢(𝑘) 𝑦(𝑘 + 2|𝑘) = 𝐶 𝑑 𝐴 2 𝑑 𝑥(𝑘) + 𝐶 𝑑 𝐴 𝑑 𝐵 𝑑 𝑢(𝑘) + 𝐶 𝑑 𝐵 𝑑 𝑢(𝑘 + 1) . . . 𝑦(𝑘 + 𝑁 𝑝 |𝑘) = 𝐶 𝑑 𝐴 𝑁𝑝 𝑑 𝑥(𝑘) +⋅ ⋅ ⋅ + 𝐶 𝑑 𝐴 𝑁𝑝-𝑁𝑐 𝑑 𝐵 𝑑 𝑢(𝑘 + 𝑁 𝑐 -1)
(2) or, in matrix form:

𝑌 = 𝐹 𝑥(𝑘) + 𝐺𝑈 ( 3
)
where 𝐹 ∈ ℝ 𝑞𝑁𝑝×𝑛 and 𝐺 ∈ ℝ 𝑞𝑁𝑝×𝑚𝑁𝑐 are defined as

𝐹 := ⎡ ⎢ ⎢ ⎢ ⎣ 𝐶 𝑑 𝐴 𝑑 𝐶 𝑑 𝐴 2 𝑑 . . . 𝐶 𝑑 𝐴 𝑁𝑝 𝑑 ⎤ ⎥ ⎥ ⎥ ⎦ 𝐺 := ⎡ ⎢ ⎢ ⎢ ⎣ 𝐶 𝑑 𝐵 𝑑 0 0 ⋅ ⋅ ⋅ 0 𝐶 𝑑 𝐴 𝑑 𝐵 𝑑 𝐶 𝑑 𝐵 𝑑 0 ⋅ ⋅ ⋅ 0 . . . 𝐶 𝑑 𝐴 𝑁𝑝-1 𝑑 𝐵 𝑑 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 𝐶 𝑑 𝐴 𝑁𝑝-𝑁𝑐 𝑑 𝐵 𝑑 ⎤ ⎥ ⎥ ⎥ ⎦
Then, 𝑌 ∈ ℝ 𝑞𝑁𝑝 is the predicted output, and 𝑈 := [𝑢(𝑘), . . . , 𝑢(𝑘 +𝑁 𝑐 -1)] is the future output along the control horizon. Now, given a reference signal 𝑊 = [𝑤(𝑘 + 1), . . . , 𝑤(𝑘 + 𝑁 𝑝 )], it is possible to write an objective function 𝐽 that accounts for the predicted control error and control effort, that is:

𝐽 = (𝑊 -𝑌 ) 𝑇 𝑄(𝑊 -𝑌 ) + 𝑈 𝑇 𝑅𝑈. ( 4
)
where 𝑄 ∈ ℝ 𝑞𝑁𝑝×𝑞𝑁𝑝 and 𝑅 ∈ ℝ 𝑚𝑁𝑐×𝑚𝑁𝑐 are suitable weighting matrices. Using the predictive equation ( 3), it is possible to rewrite (4) as a function of the sole control inputs 𝑈 :

𝐽 = ( 𝑊 -𝐹 𝑥(𝑘)-𝐺𝑈 ) 𝑇 𝑄 ( 𝑊 -𝐹 𝑥(𝑘)-𝐺𝑈 ) +𝑈 𝑇 𝑅𝑈 (5)
Therefore, at each sampling time, the following QP problem has to be solved:

minimize 1 2 𝑈 𝑇 𝐻𝑈 + 𝑓 𝑇 𝑈 subject to 𝐸 𝑈 ≤ 𝑑 (6)
where 𝐻 = 2 𝐺 𝑇 𝑄𝐺 + 𝑅, 𝑓 𝑇 = 2 (𝐹 -𝑊 ) 𝑇 𝑄𝐺, and 𝐸 𝑈 ≤ 𝑑 is a generic set of linear inequalities, which typically represent the constraints on the amplitude and slew rate of inputs and outputs. The first element of the optimal control sequence 𝑈 is applied to the system, and the procedure is repeated, in the spirit of receding horizon control. Often, the MPC problem is formulated in terms of relative increment of the control input Δ𝑢, instead of 𝑢 [START_REF] Wang | Model predictive control system design and implementation using MATLAB®[END_REF]. Given the model ( 1), this can be obtained by means of a state augmentation, which adds an integrator to the original model. Thus, defining the new state as

𝑧 := [ Δ𝑥 𝑥 ] (7)
the new model results to be:

𝑧(𝑖 + 1) = Ã𝑑 𝑧(𝑖) + B𝑑 Δ𝑢(𝑖) (8) 𝑦(𝑖) = C𝑑 𝑧(𝑖) (9)
where:

Ã𝑑 = [ 𝐴 𝑑 𝑜 𝑛×𝑛 𝐼 𝑛 𝑜 𝑛×𝑛 ] , B𝑑 = [ 𝐵 𝑑 𝑜 𝑛×𝑚 ] , C𝑑 = [ 𝑜 𝑞×𝑛 𝐶 𝑑 ] . (10)
Using the system ( Ã𝑑 , B𝑑 , C𝑑 ) instead of (𝐴 𝑑 , 𝐵 𝑑 , 𝐶 𝑑 ) the derivation of the equations remains exactly the same. In this way, the resulting controller calculates the optimized control actions in terms of Δ𝑈 . Then, the control action 𝑢(𝑖) is updated by adding the first element of the sequence Δ𝑈 and the resulting value is applied to the process.

B. Input Blocking

A method to reduce the computational complexity of the QP problem (6) consists in reducing the number of degrees of freedom of the optimization. This can be done by choosing a control horizon shorter than the prediction one. Another widespread strategy is the so-called Input Blocking [START_REF] Cagienard | Move blocking strategies in receding horizon control[END_REF], which consists in forcing the control inputs to be constant along several sampling times. The method is equivalent to parametrizing the input by means of a staircase function, in which the steps occur at multiples of the sampling period. To this end, let Û := [û 1 , . . . , û𝑀 ], with 𝑀 < 𝑁 𝑐 , be the new optimization variable, such that:

𝑈 = 𝑇 Û (11)
where 𝑇 ∈ ℝ 𝑁𝑐×𝑀 is the so-called blocking matrix. The matrix 𝑇 is a Boolean matrix, which represents the mapping between the new and the old optimization variables. For instance, considering a SISO system and aiming at implementing an IB strategy in which 𝑁 𝑐 = 3, 𝑀 = 2, and 𝑢 2 = 𝑢 3 , the blocking matrix results to be:

𝑇 = ⎡ ⎣ 1 0 0 1 0 1 ⎤ ⎦
Therefore, the QP problem ( 6) can be reformulated with respect to the new optimization variable Û , becoming:

minimize 1 2 Û 𝑇 Ĥ Û + f 𝑇 Û subject to Ê Û ≤ 𝑑 ( 12
)
where Ĥ = 𝑇 𝑇 𝐻 𝑇 , f 𝑇 = 𝑓 𝑇 , Ê = 𝐸 𝑇 . Then, the original optimization variable can be obtained as 𝑈 = 𝑇 Û and the first element of the sequence is applied to the system. This strategy permits to reduce the computational complexity of the optimization problem, since the number of variables is lowered from 𝑚𝑁 𝑐 to 𝑚𝑀 . However, the solution to the problem ( 12) is suboptimal for the problem ( 6) and the goodness of such approximation depends on the choice of the matrix 𝑇 . Moreover, note that this parametrization forces the steps of the staircase function to happen at a multiple of the sampling period and this represents a limit in the approximation, considering that common prediction horizons are typically in the order of 𝑁 𝑝 = 10. For this reason we propose to formulate the predictive model in the continuoustime domain, as it will be explained in the next section.

III. PROPOSED METHOD

As mentioned above, the IB technique is basically a stepwise parametrization of the input, in which the steps of the staircase function necessarily happen at a multiple of the sampling period. This is due to the formulation of the predictive equations in the discrete-time domain. In order to overcome this limitation, we formulate the predictive equations by using a continuous-time model, but imposing the input to be piecewise constant, as in the discrete case. In this way the staircase parametrization will result somehow more flexible, since the choice of the time-instants in which each step occurs will not be constrained by the sampling period. Note that the resulting controller will be directly implemented in discretetime, since it is based on the assumption that the input is piecewise constant.

A. Formulation of the predictive model

Consider the continuous-time state-space formulation of a given process with 𝑚 inputs, 𝑞 outputs and an 𝑛-dimensional state vector:

ẋ(𝑡) = 𝐴 𝑥(𝑡) + 𝐵 𝑢(𝑡) 𝑦(𝑡) = 𝐶 𝑥(𝑡) ( 13
)
where 𝐴 ∈ ℝ 𝑛×𝑛 , 𝐵 ∈ ℝ 𝑛×𝑚 , and 𝐶 ∈ ℝ 𝑞×𝑛 . Now, consider a time instant 𝑡 0 and the corresponding state 𝑥 0 = 𝑥(𝑡 0 ). The predicted output at a generic time 𝑡 = 𝑡 0 + 𝜏 , 𝜏 > 0, calculated at time 𝑡 0 is given by: Assume that 𝐵 is constant and 𝑢 is a staircase function with a number of steps equal to 𝑐. Let t𝑖 be the time instants, with respect to 𝑡 0 , in which the steps of the function are applied and let 𝑢(𝑡 0 + t𝑖 ) be the input value in the time interval between t𝑖-1 and t𝑖 . These intervals are defined, then, as Δ t := [ t1 -𝑡 0 , t2 -t1 , , . . . , t𝑐 -t𝑐-1 ]. A sketch of the input function 𝑢(𝑡) is drawn in Figure 1. Note that the lengths of the intervals do not need to be the same. Thus, ( 14) becomes:

𝑦(𝑡 0 + 𝜏 |𝑡 0 ) = 𝐶e 𝐴𝜏 𝑥 0 + 𝐶e 𝐴𝜏 ∫ 𝜏 0 e -𝐴𝛾 𝐵𝑢(𝛾) 𝑑𝛾 (14)
𝑦(𝑡 0 + 𝜏 |𝑡 0 ) = 𝐶e 𝐴𝜏 𝑥 0 + 𝐶 𝑐 ∑ 𝑖=1 Γ 𝑖 𝑢(𝑡 0 + t𝑖 ) (15)
where, assuming t0 = 0 for notation simplicity,

Γ 𝑖 := ⎧  ⎨  ⎩ e 𝐴(𝜏 -t𝑖) ∫ Δ t𝑖 0 e 𝐴𝛾 𝑑𝛾 𝐵, if 𝜏 ≥ t𝑖 ∫ 𝜏 -t𝑖-1 0 e 𝐴𝛾 𝑑𝛾 𝐵, if t𝑖-1 < 𝜏 < t𝑖 0, if 𝜏 < t𝑖 ∧ 𝜏 < t𝑖-1 (16)
Having defined the control time instants t𝑖 , a constant matrix Φ 𝜏 can be defined for a given prediction time 𝜏 , such that:

Φ 𝜏 := [ 𝐶 Γ 1 , 𝐶 Γ 2 , . . . , 𝐶 Γ 𝑐 ] (17)
Substituting [START_REF] Wang | Continuous time model predictive control design using orthonormal functions[END_REF] in [START_REF] Richalet | Predictive functional control: principles and industrial applications[END_REF] gives:

𝑦(𝑡 0 + 𝜏 |𝑡 0 ) = 𝐶e 𝐴𝜏 𝑥 0 + Φ 𝜏 𝑈 ( 18
)
where 𝑈 := [𝑢(𝑡 0 + t1 ), . . . , 𝑢(𝑡 0 + t𝑐 )] 𝑇 . Extending this reasoning to a generic number of prediction times (𝜏 1 , . . . , 𝜏 𝑝 ), 𝑝 > 0, we obtain the following system of equations:

⎧   ⎨   ⎩ 𝑦(𝑡 0 + 𝜏 1 |𝑡 0 ) = 𝐶e 𝐴𝜏1 𝑥 0 + Φ 𝜏1 𝑈 . . . 𝑦(𝑡 0 + 𝜏 𝑝 |𝑡 0 ) = 𝐶e 𝐴𝜏𝑝 𝑥 0 + Φ 𝜏𝑝 𝑈 ( 19
)
As for the control time instants, the prediction instants (𝜏 1 , . . . , 𝜏 𝑝 ) are not necessarily equally spaced. System [START_REF]Constrained intermittent model predictive control[END_REF] can be rewritten in matrix form as:

𝑌 = 𝐹 𝑥 0 + 𝐺𝑈 ( 20
)
where 𝑌 ∈ ℝ 𝑞𝑝×1 is the predicted output vector at times (𝑡 0 + 𝜏 1 , . . . , 𝑡 0 +𝜏 𝑝 ), 𝑈 ∈ ℝ 𝑛𝑐×1 is the vector of the control inputs, and:

𝐹 = ⎡ ⎢ ⎣ 𝐶e 𝐴𝜏1 . . . 𝐶e 𝐴𝜏𝑝 ⎤ ⎥ ⎦ ∈ ℝ 𝑞𝑝×𝑛𝑝 , 𝐺 = ⎡ ⎢ ⎣ Φ 𝜏1 . . . Φ 𝜏𝑝 ⎤ ⎥ ⎦ ∈ ℝ 𝑞𝑝×𝑚𝑐 .
In classic discrete MPC, prediction and control instants are equally distributed along the respective horizons and correspond to the sampling times. However, although the derivation of the predictive equations in continuous time might seem less intuitive, its great advantage consists in the possibility of using a free distribution of prediction and control instants along the horizons, disregarding the real sampling period of the system. As a matter of fact, it is not necessary for the time instants 𝜏 𝑖 and t𝑖 to be neither equally distributed nor coincident with a sampling time-instant. Anyhow, note that the controller can be implemented using a receding horizon approach with a given sampling period 𝑇 , under the only condition t1 ≥ 𝑇 (as better explained in Section IV).

B. Handling of constraints

Constraints expressed as linear combinations of the control actions can be included in the MPC framework, just like in the classic MPC formulation, resulting in a constrained LQ problem.

C. Optimal control action

As mentioned in Section II, in classic linear MPC, the control sequence is obtained by minimizing the cost function (4). In the continuous-time domain this can be rewritten as:

𝐽 = ∫ 𝑡0+𝜏𝑝 𝑡0 ( 𝑦(𝑡) -𝑤(𝑡) ) 2 𝑑𝑡 + 𝜆 ∫ t𝑐 𝑡0 𝑢 2 (𝑡)𝑑𝑡 (21)
This cost function can be approximated as:

𝐽 ≈ (𝑌 -𝑊 ) 𝑇 Δ𝜏 (𝑌 -𝑊 ) + 𝜆𝑈 𝑇 Δ𝑡 𝑈 (22)
where 𝑊 ∈ ℝ 𝑞𝑝×1 is the reference signal at times (𝑡 0 + 𝜏 1 , . . . , 𝑡 0 + 𝜏 𝑝 ), Δ𝜏 = diag(𝜏 1 , 𝜏 2 -𝜏 1 , . . . , 𝜏 𝑝 -𝜏 𝑝-1 ), Δ𝑡 = diag(Δ𝑡 1 , . . . , Δ𝑡 𝑐 ), and 𝜆 > 0.

Note that ( 22) is derived from ( 21) by substituting the integrals with a discrete approximation, with step size equal to the intervals (𝜏 1 , 𝜏 2 -𝜏 1 , . . . , 𝜏 𝑝 -𝜏 𝑝-1 ). Therefore, using [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF], the following QP results (in terms of the sole control action 𝑈 ∈ ℝ 𝑚𝑐×1 ):

minimize 1 2 𝑈 𝑇 𝐻𝑈 + 𝑓 𝑇 𝑈 subject to 𝐸 𝑈 ≤ 𝑑 ( 23
)
where 𝐻 = 2 𝐺 𝑇 Δ𝜏 𝐺 + 𝜆Δ𝑡, 𝑓 𝑇 = 2 (𝐹 -𝑊 ) 𝑇 Δ𝜏 𝐺, and 𝐸 𝑈 ≤ 𝑑 is a generic set of linear inequalities, as in [START_REF] Bartlett | Quadratic programming algorithms for large-scale model predictive control[END_REF].

D. Formulation in terms of Δ𝑈

As already mentioned, classic MPC controllers are usually formulated in terms of increments of the control action (i.e., Δ𝑢 instead of 𝑢). The first optimal Δ𝑢 is then added to the previous control action and this value is applied to the system. Given a discrete-time state-space system, this formulation can be achieved by means of a state augmentation, as mentioned in Section II. In our case, an analogous state augmentation of the continuous-time model [START_REF] Tøndel | Complexity reduction in explicit model predictive control[END_REF] would lead to a formulation of the predictive equations in which the derivative of the inputs is a staircase function. This means that the control input would result to be a piecewise linear function. However, we want to formulate the predictive equations by imposing the increment of the control action (Δ𝑢) to be piecewise constant along the correspondent time-intervals Δ t𝑖 (as in discrete MPC). To this aim, the formulation retraces the derivation of the predictive equations from [START_REF] Qin | A survey of industrial model predictive control technology[END_REF]. By assuming 𝑢 = 𝑢 0 + Δ𝑢(𝑡), where 𝑢 0 = 𝑢(𝑡 0 ) and Δ𝑢 is a staircase function with 𝑐 steps, it results:

𝑦(𝑡 0 + 𝜏 |𝑡 0 ) = 𝐶e 𝐴𝜏 𝑥 0 + Ω 𝜏 𝑢 0 + 𝐶 𝑐 ∑ 𝑖=1 Γ𝑖 Δ𝑢(𝑡 0 + t𝑖 ) (24)
where Δ𝑢(𝑡 0 + t𝑖 ) = 𝑢(𝑡 0 + t𝑖 ) -𝑢(𝑡 0 + t𝑖-1 ) and

Ω 𝜏 := 𝐶 ∫ 𝜏 -𝑡0 0 e 𝐴𝛾 𝑑𝛾 𝐵 Γ𝑖 := { ∫ 𝜏 -t𝑖-1 0 e 𝐴𝛾 𝑑𝛾 𝐵, if 𝜏 > t𝑖-1 0, if 𝜏 ≤ t𝑖-1
Therefore:

Φ𝜏 := [ 𝐶 Γ1 , 𝐶 Γ2 , . . . , 𝐶 Γ𝑐 ] ( 25
)
and finally:

𝑌 = F [ 𝑥 0 𝑢 0 ] + GΔ𝑈 (26)
where Δ𝑈 = [Δ𝑢(𝑡 0 + t1 ), . . . , Δ𝑢(𝑡 0 + t𝑐 ) ] 𝑇 and

F = ⎡ ⎢ ⎣ 𝐶e 𝐴𝜏1 , Ω 𝜏1 . . . . . . 𝐶e 𝐴𝜏𝑝 , Ω 𝜏𝑝 ⎤ ⎥ ⎦ ∈ ℝ 𝑞𝑝×𝑛𝑝 , G = ⎡ ⎢ ⎣ Φ𝜏1 . . . Φ𝜏𝑝 ⎤ ⎥ ⎦ ∈ ℝ 𝑞𝑝×𝑚𝑐 .
Once the new predictive equations have been derived, the optimal control problem can be set up as in Section III-C, but in terms of Δ𝑈 . Then, at each time-step, the current value of the control action 𝑢 0 will be updated by summing the first element of the optimal sequence Δ𝑈 and will be applied to the system.

IV. TUNING OF THE PREDICTION AND CONTROL HORIZONS

The length of the prediction horizon and the weight on the control action 𝜆 can be chosen by following standard MPC tuning rules (e.g. [START_REF] Garriga | Model predictive control tuning methods: A review[END_REF]). Therefore, the main issue in the practical application of the proposed method is the choice of an optimal distribution of the prediction and control timeinstants along the respective horizons. As a matter of fact, in discrete MPC, such instants are fixed and correspond to the future sampling times. However, the proposed method requires to explicit the time-instants 𝜏 𝑖 and t𝑗 , necessary for the approximation of the integral cost function [START_REF] Wang | Model predictive control system design and implementation using MATLAB®[END_REF] and the parametrization of the control inputs, respectively. This choice must be made during the design of the controller, then the distribution of time-instants will always be the same.

The tuning strategy is based on the assumption that the proposed approximated control law should mimic at best the FDOF controller (i.e. the discrete-time MPC without any approximation). Thus, given a prediction horizon 𝑇 𝑝 in terms of seconds, the number of prediction time-instants 𝑝 along the predictive horizon, and the number of steps 𝑐 of the approximated input function; the distribution of 𝜏 𝑖 and t𝑗 should minimize the integral of the value of the difference between the control laws given by the FDOF controller and the proposed one, that is: minimize

∫ 2 𝑇𝑝 0 ( 𝑢 𝐹 𝐷𝑂𝐹 (𝑡) -𝑢(𝑡) ) 𝑑𝑡 subject to ∑ 𝑝 𝑖=1 𝜏 𝑖 ≤ 𝑇 𝑝 ; 𝜏 𝑖 ≥ 𝜏 𝑖-1 ∑ 𝑐 𝑗=1 t𝑗 ≤ 𝑇 𝑝 ; t𝑗 ≥ t𝑗-1 t1 ≥ 𝑇 𝑠 ( 27
)
where 𝜏 𝑖 and t𝑗 are the optimization variables, while 𝑢 𝐹 𝐷𝑂𝐹 and 𝑢 are the control inputs given by the FDOF and the proposed controller, respectively. 𝑇 𝑠 represents the sampling period at which the receding horizon is implemented. The upper integration extreme was chosen equal to 2 𝑇 𝑝 , in order to account for all the transient in the minimization (since 𝑇 𝑝 is generally shorter than the rise time of the system). Moreover, note that the last condition t1 ≥ 𝑇 𝑠 is necessary for the proposed controller to be feasible for the use of receding horizon with a given sampling period 𝑇 𝑠 . As a matter of fact, in order for this to happen, the successive sampling time must occur within the first step of the input parametrization.

Of course the goodness of the approximation depends on 𝑐 and 𝑝: the closer they are to the FDOF's values, the smaller is the deviation between the two inputs. Figure 2 shows the Pareto curve of the objective function of Problem (27) over the number of control steps 𝑐 for the case study (28), given a fixed number of prediction instants 𝑝 = 10. As expected, the cost value decreases for 𝑐 getting closer to the FDOF's value.

It is important to point out that a limitation of this tuning approach is that the choice of the prediction and control timeinstants is optimal only with respect to the reference signal used to obtain the control input 𝑢 𝐹 𝐷𝑂𝐹 . In this work only step responses were considered. The same limitation applies to the presence of disturbances or unmodeled dynamics, since they are not considered in the generation of 𝑢 𝐹 𝐷𝑂𝐹 .

Note that the optimization problem (27) can be addressed by means of gradient-based algorithms (e.g. fmincon in Matlab). As a matter of fact, the optimization variables can span continuously the prediction horizon, without any constraint given by the sampling period. On the contrary, the analogous discrete optimization problem would constrain the time-instants to be a multiple of the sampling period. This happens, for instance, when the same tuning rule is applied to IB: the input is forced to be constant along a whole number of sampling periods, leading to a discrete optimization problem. The possibility of using gradient-based algorithms is an advantage in case a high number of optimization variables is involved, since the discrete problem should be solved by means of combinatorial or stochastic methods (e.g. genetic algorithms) and this could be much more time-consuming.

V. NUMERICAL RESULTS

As an illustrative example, consider the fast-dynamics second-order system

𝐺(𝑠) = 2500 (𝑠 + 50) 2 (28)
and a sampling time 𝑇 𝑠 = 0.01 s. In order to control the system by means of a model predictive algorithm, consider a prediction horizon 𝑁 𝑝 = 10, accordingly with common tuning rules used in MPC practice [START_REF] Garriga | Model predictive control tuning methods: A review[END_REF]. In terms of seconds, this corresponds to a prediction horizon equal to 𝑇 𝑝 = 0.1 s. To prove the effectiveness of the proposed method, four MPC strategies are analyzed:

• FDOF controller, that is, classic MPC with 𝑁 𝑐 = 𝑁 𝑝 (this controller represents the reference case for the following approximated strategies); • classic MPC with 𝑁 𝑐 < 𝑁 𝑝 ;

• Input Blocking (IB); • the proposed method (in the remainder of the section it will also be referred to as Continuous Input Blocking (CIB)). In particular, the aim of this section is the comparison of the last three methods (𝑁 𝑐 < 𝑁 𝑝 , IB, CIB), in order to demonstrate which one ensures the best approximation of the FDOF case, being equal the reduction of computational complexity (i.e., with the same number of optimization variables). Therefore, given the same number of control instants along the horizon, the three approximated controllers are set as follows:

• CIB: the time-instants 𝜏 𝑖 and t𝑗 are chosen by following the tuning rule (27). • IB: similarly to the previous case, the Blocking Matrix 𝑇 is selected by minimizing the same objective function as in (27). • 𝑁 𝑐 < 𝑁 𝑝 : since it is the classic MPC with 𝑁 𝑐 < 𝑁 𝑝 , the control instants correspond to the first 𝑁 𝑐 sampling times. Figure 3 shows the different choices in the distribution of the prediction and control time-instants between the IB and the CIB strategies following the tuning rule (27) for a number of control moves equal to 3 (i.e. 𝑐 = 𝑀 = 3) and with 𝜆 = 5. Of course the prediction instants in the IB technique are the multiples of the sampling period, while for the CIB method they can be chosen arbitrarily. Similarly, the control instants in the CIB method can be chosen freely, whereas in IB the control inputs must be kept constant along a whole number of sampling periods (in this case it results 𝑢 2 = 𝑢 3 and 𝑢 4 = 𝑢 5 = 𝑢 6 ).

All the controllers are implemented in terms of relative increments of the control input (i.e. in terms of Δ𝑢) and implementing the receding horizon approach. The simulations are performed for a number of control moves (𝑐, 𝑀, 𝑁 𝑐 in the corresponding cases) equal to 2 and 3. For each controller, the goodness of the approximation of the FDOF case is evaluated by means of the following indices:

𝜂 1 := IAE(𝑦 -𝑦 𝐹 𝐷𝑂𝐹 ) IAE(𝑤 -𝑦 𝐹 𝐷𝑂𝐹 ) ⋅ 100 (29)
𝜂 2 := √ ISE(𝑦 -𝑦 𝐹 𝐷𝑂𝐹 ) ISE(𝑤 -𝑦 𝐹 𝐷𝑂𝐹 ) ⋅ 100 (30)
𝜂 3 := IAE(𝑢 -𝑢 𝐹 𝐷𝑂𝐹 ) IAE(𝑢 ∞ -𝑢 𝐹 𝐷𝑂𝐹 ) ⋅ 100 (31) 𝜂 4 := √ ISE(𝑢 -𝑢 𝐹 𝐷𝑂𝐹 ) ISE(𝑢 ∞ -𝑢 𝐹 𝐷𝑂𝐹 ) ⋅ 100 ( 32
)
where 𝑦 is the step-response of the evaluated method, 𝑦 𝐹 𝐷𝑂𝐹 is the step-response in the FDOF case, 𝑤 is the unitary step, 𝑢 is the control input of the evaluated method, 𝑢 𝐹 𝐷𝑂𝐹 is the control input given by the FDOF controller, 𝑢 ∞ is the steady-state value of the control input, and IAE(⋅) and ISE(⋅) calculate the integral of the absolute value and the integral of the square value of a signal, that is:

IAE(𝑧) := ∫ ∞ 0 |𝑧(𝑡)|𝑑𝑡 ISE(𝑧) := ∫ ∞ 0 ( 𝑧(𝑡) ) 2 𝑑𝑡
The indices 𝜂 1 and 𝜂 2 measure how well the step-response of the method under evaluation approximates the FDOF stepresponse. Similarly, 𝜂 3 and 𝜂 4 evaluate how well the control action mimics the FDOF control law. Such indices can be considered as a percentage decay-measure with respect to the ideal case: the smaller is the index value, the better is the approximation.

The results of the simulations are reported in Table I and II. The weight on the control action was chosen as 𝜆 = 5 and . Most important, compared to the other strategies, the proposed method gives better values for all the indices and for all the considered cases. This is also visible in Figures 4 and6, which show the step response of the four controllers, with 𝜆 = 5 and 0.2 respectively. As expected, the proposed controller (solid black line) follows the FDOF one (dashed grey line) more accurately. Similarly, Figures 5 and7 show the control inputs given by the four controllers, with 𝜆 = 5 and 0.2, respectively. As in the previous case, the proposed method (solid black line) is evidently closer to the reference control signal (dashed grey line). The use of a number of control instants greater than 3 would lead to similar considerations. However, the use of 𝑐 > 3 for this case study would not give significant improvements, whereas it would unnecessarily increase the computational burden of the algorithm. As a matter of fact, the Pareto curve depicted in Figure 2 shows that a greater number of control instants does not lead to a remarkable reduction in the minimum of the objective function in (27). Step-response comparison between the methods with two control instants (𝜆 = 0.2). FDOF (dashed grey), CIB (solid black), IB (dotted blue), and 𝑁𝑐 < 𝑁𝑝 (dashed-dotted red).

VI. CONCLUSIONS

This paper presented a new method to reduce the computational complexity of MPC with online optimization. The core of the proposed method is the formulation of the predictive equations and the parametrization of the inputs in the continuous-time domain. The control sequence is forced to be piecewise constant along the control horizon, but with less steps than in the FDOF case. The spirit of such parametrization is similar to the widespread IB technique. However, the continuous-time formulation permits to set the prediction and control time-instants disregarding the sampling period of the system. This ensures a better flexibility of the parametrization with respect to IB. Note that the staircase parametrization of the input permits to consider the resulting controller directly in the discrete-time domain.

The optimal tuning of the proposed method has been addressed by minimizing the deviation between the approximated controller and the ideal one (i.e. discrete MPC without any parametrization). It should be pointed out that such minimization can be performed by means of gradient-based algorithms, since the optimization variables can continuously span the time-interval. This is an advantage in case a high number of variables is involved in the optimization. As a matter of fact, the discrete counterpart of such minimization should be carried out by means of combinatorial or stochastic algorithms, which could be much more time-consuming and difficult to tune.

The analysis of a case study has demonstrated the effectiveness of the proposed method against two common strategies for complexity reduction (IB and control horizon shorter than the predictive one). The numerical results show that the proposed method guarantees a better approximation of the ideal controller, in terms of both output response and control law.

Future work will include the analysis of the effects of modeling uncertainties and disturbances on the optimal tuning methodology.

Fig. 1 .

 1 Fig. 1. Sketch of the parametrization of the input as a staircase function.

Fig. 2 .

 2 Fig. 2. Pareto curve of the cost function (27) over 𝑐, with 𝑝 = 10, for the system (28). The y-axis is in logarithmic scale.

Fig. 3 .

 3 Fig. 3. Different choices in the distribution of the prediction and control time-instants between CIB and IB (𝜆 = 5).

Fig. 6 .

 6 Fig.6. Step-response comparison between the methods with two control instants (𝜆 = 0.2). FDOF (dashed grey), CIB (solid black), IB (dotted blue), and 𝑁𝑐 < 𝑁𝑝 (dashed-dotted red).

TABLE I NUMERICAL

 I COMPARISON OF THE APPROXIMATED METHODS (CIB = PROPOSED METHOD, IB = INPUT BLOCKING, 𝑁𝑐 < 𝑁𝑝 = CLASSIC MPC WITH SHORTER CONTROL HORIZON) WITH 𝜆 = 5 .

	N. control Method instants	𝜂 1	𝜂 2	𝜂 3	𝜂 4
		CIB	0.91	0.70	2.89	1.74
	2	IB	1.52	1.30	5.17	3.71
		𝑁𝑐 < 𝑁𝑝 6.73	5.99	23.4	1 7 .6
		CIB	5 ⋅ 10 -5 4 ⋅ 10 -5 1.5 ⋅ 10 -4 8 ⋅ 10 -5
	3	IB	0.55	0.47	1.74	1.26
		𝑁𝑐 < 𝑁𝑝 2.57	2.18	9.14	6.51
				TABLE II
	NUMERICAL COMPARISON OF THE APPROXIMATED METHODS (CIB =
	PROPOSED METHOD, IB = INPUT BLOCKING, 𝑁𝑐 < 𝑁𝑝 = CLASSIC MPC
		WITH SHORTER CONTROL HORIZON) WITH 𝜆 = 0.2 .
	N. control Method instants	𝜂 1	𝜂 2	𝜂 3	𝜂 4
		CIB	8 ⋅ 10 -5 5 ⋅ 10 -5 2.4 ⋅ 10 -4 1.4 ⋅ 10 -4
	2	IB	3.66	2.95	13.5	1 2 .8
		𝑁𝑐 < 𝑁𝑝 6.24	5.10	14.6	1 3 .0
		CIB	6 ⋅ 10 -5 4 ⋅ 10 -5 2.0 ⋅ 10 -4 1.3 ⋅ 10 -4
	3	IB	1.19	0.87	4.54	4.24
		𝑁𝑐 < 𝑁𝑝 1.36	1.06	4.83	4.54

𝜆 = 0.2 respectively. As expected, the classic strategy with 𝑁 𝑐 < 𝑁 𝑝 generally leads to worse results with respect to the other two methods (being equal the number of optimization variables)

  Fig.5. Control-input comparison between the methods with two control instants (𝜆 = 5). FDOF (dashed grey), CIB (solid black), IB (dotted blue), and 𝑁𝑐 < 𝑁𝑝 (dashed-dotted red).

		1.1
		1
		0.9
		0.8
		0.7			0.75
	y	0.6 0.5			0.7
		0.4			0.65
		0.3
		0.2			0.6
		0.1			0.05 0.55		0.055		0.06
		0
		0	0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
					Time [s]
	Fig. 4. Step-response comparison between the methods with two control
	instants (𝜆 = 5). FDOF (dashed grey), CIB (solid black), IB (dotted blue),
	and 𝑁𝑐 < 𝑁𝑝 (dashed-dotted red).
		1.1
		1
		0.9
		0.8			1.1
	u	0.7
					1
		0.6
		0.5			0.9
		0.4			0.8
					0.02	0.03	0.04	0.05	0.06
		0.3
		0	0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
					Time [s]
		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14
					Time [s]

  Fig.7. Control-input comparison between the methods with two control instants (𝜆 = 0.2). FDOF (dashed grey), CIB (solid black), IB (dotted blue), and 𝑁𝑐 < 𝑁𝑝 (dashed-dotted red).

		1.9
		1.8
		1.7
		1.6
		1.5
	u	1.4
		1.3
		1.2
		1.1
		1
		0.9
		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14
					Time [s]

ACKNOWLEDGMENT

This work has been partially funded by the National R+D+i Plan Project DPI2014-56364-C2-1-R of the Spanish Ministry of Economy and Competitiveness and ERDF funds.