João Cavalcanti Santos

Member, IEEE Marc Gouttefarde
email: gouttefarde@lirmm.fr

A Nonlinear Model Predictive Control for the Position Tracking of Cable-Driven Parallel Robots

Keywords:

The present paper proposes a Nonlinear Model Predictive Control (NMPC) strategy for the position tracking of Cable-Driven Parallel Robots (CDPRs). The NMPC formulation handles explicitly the cable tensions and their limits. Accordingly, the cable tension distribution is performed as an integral part of the NMPC feedback control strategy, which notably allows the CDPR to operate on the wrench-feasible workspace boundaries without failure. In order to integrate the cable tension minimization within the NMPC formulation, the concept of Wrench Equivalent Optimality (WEO) is introduced. The WEO is a non-negative measure able to evaluate if the wrench generated by a given cable tension vector can be generated by an alternative tension vector with smaller 2-norm. The redundancy resolution performed by means of the minimization of the WEO enables the stability of the closed-loop system to be proved. More precisely, sufficient conditions for the uniform asymptotic stability are deduced using results from the analysis of NMPC schemes without terminal constraints and costs. Furthermore, the proposed NMPC strategy is validated experimentally on a fully-constrained 6 degree-of-freedom CDPR.

I. INTRODUCTION

P OSITION tracking control of Cable-Driven Parallel Robots (CDPRs) addresses the problem of generating a mobile platform trajectory as close as possible to a desired one. The majority of previous works dealing with this problem implemented model-based strategies including a feedforward input based on the system dynamics [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF]. Acceptable results can notably be obtained with a PID feedback correction combined to computed torque control [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] or with a similar approach with position-controlled winches [START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF]. Nevertheless, advanced control techniques may lead to improved results. Sliding Mode Control (SMC) [START_REF] Slotine | Sliding controller design for non-linear systems[END_REF] received particular attention in this context [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF]- [START_REF] Oh | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF]. SMC is well-known to present easy implementation, good robustness against modeling uncertainties and, generally, finite time convergence. Babaghasabha et al. [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF] applied this control strategy to a planar CDPR. Terminal sliding mode [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF] and super twisting controller [START_REF] Schenk | Application of a differentiator-based adaptive super-twisting controller for a redundant cable-driven parallel robot[END_REF] are improved versions of SMC that were used to control CDPRs. Some other examples of advanced control techniques used with the same goal are H ∞ [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF], [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF], Lyapunov-based controllers [START_REF] Oh | Generation of feasible set points and control of a cable robot[END_REF], [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF] , and learning-based tracking controllers [START_REF] Li | Tracking control of fully-constrained cabledriven parallel robots using adaptive dynamic programming[END_REF].

Although minimum and maximum cable tension limits play a crucial role in the operation of a CDPR, the control strategies proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF] are not able to handle these constraints as an integral part of the feedback motion controller. These constraints are typically handled in a cable tension distribution algorithm, e.g. [START_REF] Ueland | Optimal Force Allocation for Overconstrained Cable-Driven Parallel Robots: Continuously Differentiable Solutions With Assessment of Computational Efficiency[END_REF]- [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF], which takes place once the desired wrench is already computed by the controller. As a consequence, an unfeasible desired wrench taken as the input to a tension distribution algorithm generally leads to a control breakdown. This kind of feasibility issue typically arises close to the boundaries of the CDPR wrench-feasible workspace (WFW) [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF] or in presence of significant modeling uncertainties or disturbances.

To overcome this problem, a Model Predictive Control (MPC) strategy was introduced in our preliminary work [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF]. As highlighted in [START_REF] Maciejowski | Predictive control: with constraints[END_REF], MPC is one of the few control strategies able to handle system constraints explicitly. Thus, the Linear MPC (LMPC) introduced in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] integrates the cable tension distribution calculation within its formulation. Experimental results showed that this MPC is not prone to the aforementioned feasibility issues since it can operate on the WFW boundary without failure.

In spite of the significant improvements obtained in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], the stability analysis of the corresponding closed-loop system is hindered by the linear approximation of the actual CDPR nonlinear dynamics. Furthermore, the obtained tracking error leaves room for improvement. In this context, the Nonlinear Model Predictive Control (NMPC) strategy introduced in the present paper seeks to ally superior tracking precision and a closed-loop stability proof, in addition to the aforementioned ability to work on the WFW boundary without failure.

Apart from [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], studies addressing MPC schemes for CDPRs are very few [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF], [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF]. Katliar et al. proposed in [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF] an NMPC to generate a desired set of accelerations and velocities for a motion simulator. Its performance was investigated through numerical simulations while the closedloop system stability and the ability to work close to the WFW boundaries were not addressed. Moreover, the present paper deals with trajectory tracking and prioritizes positioning accuracy. Numerical simulations and experimental results are presented in [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF]. However, the proposed LMPC is meant for vibration attenuation of a 2-DoF planar CDPR where only some of the cables are controlled in tension and the closedloop system stability is not addressed. The NMPC proposed in the present paper is compatible with any CDPR actuated by at least as many cables as number of DoFs. Experimental validation is done on a 6-DoF CDPR whose eight cable tensions are controlled by the NMPC scheme in real time. Additionally, the focus of the present work is position tracking control, instead of vibration attenuation.

Regarding the stability of systems controlled with MPC schemes, [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] advocates that terminal costs and constraints (shortly referred to as terminal conditions) represent a powerful tool since the stability of the closed-loop system is deduced as a direct consequence of Lyapunov stability theory [START_REF] Rawlings | Model predictive control: Theory, Computation, and Design[END_REF]Section 2.4]. Accordingly, numerous applications of MPC in robotics used terminal conditions [START_REF] Chipalkatty | Less Is More: Mixed-Initiative Model-Predictive Control With Human Inputs[END_REF]- [START_REF] Heshmati-Alamdari | A self-triggered visual servoing model predictive control scheme for under-actuated underwater robotic vehicles[END_REF]. Nevertheless, the use of terminal conditions in an MPC scheme may lead to important drawbacks such as terminal constraint construction issues, increased computational burden, and operating range limitations [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]. More specifically, terminal constraints may render an MPC scheme more prone to unfeasible optimization problems in particular when operating close to system constraints, which goes against the ability to work close to the WFW boundary without failure such as demonstrated in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF].

Advantages and drawbacks of terminal conditions are further discussed in [START_REF] Grüne | NMPC without terminal constraints[END_REF], [START_REF] Mayne | An apologia for stabilising terminal conditions in model predictive control[END_REF] and [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 7.4].

Alternatively, MPC schemes without terminal conditions can be used. The formulation and real-time solution of the corresponding optimization problems without additional constraints are simplified. However, more advanced techniques are necessary in order to guarantee the stability of the closed-loop system since it cannot be deduced as a direct consequence of Lyapunov stability theory. Hence, several studies used MPC schemes without terminal conditions in robotics [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF], [START_REF] Khadem | Autonomous Steering of Concentric Tube Robots via Nonlinear Model Predictive Control[END_REF]- [START_REF] Wen | The study of model predictive control algorithm based on the force/position control scheme of the 5-DOF redundant actuation parallel robot[END_REF], but very few performed a stability analysis, e.g., Worthmann et al. [START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF] use a non-quadratic stage cost for the NMPC without terminal conditions of non-holonomic mobile robots and attain a stable closed-loop system. Position tracking control of CDPRs presents paramount differences with the steering of non-holonomic mobile robots so that the results obtained in [START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF] are not applicable in the present paper. Nevertheless, similarly to [START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF], the proposed NMPC scheme stability analysis will be based on [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 6].

In order to apply the theory of Grüne and Pannek [START_REF] Grüne | Nonlinear model predictive control[END_REF], the NMPC scheme proposed in this paper applies an original formulation of the cost functional. Many MPC applications in robotics use the penalization of (i) the tracking errors, (ii) the amplitude of the control inputs and/or (iii) the variation of the control inputs. Typically, the corresponding cost function consists of one of these three variables or of a weighted sum of two of them [START_REF] Incremona | MPC for Robot Manipulators With Integral Sliding Modes Generation[END_REF], [START_REF] Khadem | Autonomous Steering of Concentric Tube Robots via Nonlinear Model Predictive Control[END_REF], [START_REF] Jurado Realpe | Model Predictive Controller for a Planar Tensegrity Mechanism with Decoupled Position and Stiffness Control[END_REF], [START_REF] Kayacan | Robust Trajectory Tracking Error Model-Based Predictive Control for Unmanned Ground Vehicles[END_REF]- [START_REF] Wen | The study of model predictive control algorithm based on the force/position control scheme of the 5-DOF redundant actuation parallel robot[END_REF]. Indeed, this was the case for the above mentioned MPC schemes used to control CDPRs [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF], [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF], in which the cable tensions were taken as the control inputs. However, for the position tracking control of CDPRs, this formulation is not consistent with the methodology used for the stability analysis of NMPC schemes presented in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 6] and may lead to deteriorated tracking errors. Indeed, an NMPC scheme considering (ii) and/or (iii) is generally not able to track a feasible desired trajectory with null tracking errors since the optimization of a cost functional consisting of a weighted sum of (i)-(iii) leads to a compromise between these three aspects, which generally have different individual minima. In other words, the minimization of a weighted sum of (i)-(iii) seeks to minimize (ii) and/or (iii) to the detriment of (i), the tracking error. The control scheme proposed in this work provides a solution to this issue.

The main contribution of the present paper is the formulation of an NMPC scheme for the position tracking control of CDPRs using the concept of Wrench Equivalent Optimality (WEO). For a given platform pose and a set of cable tensions, the WEO is a non-negative measure that evaluates whether the wrench generated by these cable tensions can be generated by an alternative cable tension vector with smaller 2-norm. The penalization of the WEO is consistent with the minimization of the tracking errors in the sense that both tracking errors and WEO can be null. By means of this original formulation, the cable tension distribution is implicitly performed in real time and the stability of the closed-loop system is analyzed using methods for the analysis of NMPC schemes without terminal conditions. In addition, since the cable tensions are explicitly handled in the NMPC controller, the proposed strategy allows the CDPR to operate on the boundary of the WFW without failure. Finally, compared to the state-of-the-art LMPC scheme [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], the proposed strategy leads to improved tracking errors.

The tracking error improvement is validated experimentally on a 6-DoF CDPR whose mobile platform follows a typical pick-and-place trajectory. In these experiments, in order to show that the proposed scheme can operate on the WFW boundary without failing, the same pick-and-place path is performed with maximum allowed tensions drastically reduced. The desired trajectory escapes from the WFW defined with these new cable tension limits. The proposed NMPC is able to track a trajectory as close as possible to the desired one without violating the cable tension limits, some tensions attaining their maximum limit. In order to highlight the relevance of this feature, the same trajectory is performed with the control scheme proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]. In this case, as soon as the desired trajectory reaches an unfeasible pose, the CDPR operation is terminated because the controller requires the tension distribution module to generate an unfeasible wrench. Note that this infeasibility issue would also be obtained with advanced control schemes such as robust [START_REF] Korayem | Dynamic load carrying capacity of flexible cable suspended robot: Robust feedback linearization control approach[END_REF], [START_REF] Khosravi | Robust PID control of fullyconstrained cable driven parallel robots[END_REF] and adaptive [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF]- [START_REF] Schenk | Application of a differentiator-based adaptive super-twisting controller for a redundant cable-driven parallel robot[END_REF] strategies.

This paper is organized as follows. Section II presents the CDPR continuous and discrete-time dynamic models. Section III introduces the proposed NMPC scheme. In Section IV, the stability of the corresponding closed-loop system is analyzed. In Section V, the proposed NMPC scheme is compared to an LMPC in numerical simulations. Experimental results are presented in Section VI and conclusions are drawn in Section VII. Appendices A, B and C discuss some details on the dynamic model and on the stability analysis. Appendix D briefly describes the used numerical optimization methods.

II. PRELIMINARIES AND DYNAMIC MODELS

This section presents the continuous dynamic model of a CDPR consisting of an n-DoF mobile platform driven by m cables, where n m. An example of a 6-DoF CDPR driven by 8 cables is depicted in Figure 2. As for most of the typical applications of MPC in robotics, the control scheme proposed in Section III is based on a discrete-time system. Accordingly, a discrete-time CDPR dynamic model based on its continuous counterpart is also presented.

The platform pose is given by the vector-valued function

x x x : R → R n , so that x x x(t) = p p p(t) T ψ ψ ψ(t) T T is the pose vector containing both the platform position p p p(t) and orientation ψ ψ ψ(t) at instant t. These vectors are written with respect to the fixed reference frame O r , as depicted Figure 2. Typically, ψ ψ ψ(t) consists of Euler angles. The dependence on time is dropped leading to the shorthand notation x x x whenever there is no risk of confusion. The same shorthand notation is used for other vector-valued and matrix-valued functions in the remainder of this paper. The common dot notation is used to refer to the time derivatives, so that dx x x/dt = ẋ x x and d 2 x x x/dt 2 = ẍ x x. The vector-valued function τ τ τ c : R → R m represents the vector of cable tensions in function of time

τ τ τ c (t) = τ 1 (t) . . . τ m (t) T
. Each cable force is applied on the corresponding attachment point of the mobile platform (see Figure 1). The mobile platform equations of motion can be obtained with the Newton-Euler formalism, which, neglecting cable distributed mass, leads to the following dynamic model

M(x x x) ẍ x x + C(x x x, ẋ x x) ẋ x x = g(x x x) + W(x x x) τ τ τ c , (1)
where M(x x x) is the mass matrix, C(x x x, ẋ x x) is the vector of Coriolis and centripetal forces, g(x x x) is the vector of gravitational forces and W(x x x) is the wrench matrix. Details on the continuous dynamic model ( 1) are discussed in Appendix A.

Based on the nonlinear continuous-time system (1), it is necessary to define the discrete-time dynamic model that will be used in Section III. This discrete-time model is written in terms of sequences of vectors, which are denoted with bold (non-italic) lowercase letters followed by (•) (such as, for instance, s(•)). The set of sequences of vectors with dimension n s and infinite number of elements is denoted as S n s . More precisely,

S n s = {s(•) : N 0 → R n s }, (2)
with N 0 = N ∪ {0} the set of non-negative integers and N the set of strictly positive integers. Accordingly, the k th vector of a sequence s(•) ∈ S n s is denoted as

s(k) = s k ∈ R n s .
Consider the sequences of vectors x(•), ẋ(•), ẍ(•) ∈ S n and y(•) ∈ S 2n defined as

x(k) = x k = x x x(t 0 + k ∆t) ẋ(k) = ẋk = ẋ x x(t 0 + k ∆t) ẍ(k) = ẍk = ẍ x x(t 0 + k ∆t) y(k) = y k = x T k ẋT k T (3)
for all k ∈ N 0 , an initial time t 0 and a sampling period ∆t. Without loss of generality, since the continuous system (1) is time-invariant (according to the definition in [51, Chapter 1]) the initial time is considered to be t 0 = 0. The subsequent time instants are denoted as

t 1 = ∆t, t 2 = 2 ∆t, . . . , t k = k ∆t, for k ∈ N 0 .
Similarly, considering a digital control approach, the continuous-time representation τ τ τ c of the actual cable tensions is considered piece-wise constant and given by its discrete counterpart τ τ τ(•) ∈ S m , i.e.:

τ τ τ(k) = τ τ τ k = τ τ τ c (t), for k ∆t < t (k + 1) ∆t and k ∈ N 0 . (4)
Typically, every cable tension vector τ τ τ k , ∀ k ∈ N 0 should satisfy τ τ τ min τ τ τ k τ τ τ max , for constant τ τ τ min , τ τ τ max ∈ R m + and

R m + = {v ∈ R m | v > 0}.
The set of admissible cable tensions U ⊂ R m + is thus defined as

U = {τ τ τ ∈ R m | 0 < τ τ τ min τ τ τ τ τ τ max }. (5)
Consider a sequence of feasible cable tensions τ τ τ(•) ∈ S m represented by an infinite number of vectors {τ τ τ 0 , τ τ τ 1 , . . . }. For a given time instant k ∈ N 0 , special attention will be devoted to the finite sequence of cable tensions {τ τ τ k , τ τ τ k+1 , . . . , τ τ τ k+h p -1 }, with a positive integer h p . Accordingly, the following set of truncated sequences with a finite number of feasible cable tension vectors is introduced:

U h p k = τ τ τ(•) : N k,k+h p -1 → U , (6)
where,

N n 1 ,n 2 = {i ∈ N 0 | n 1 i n 2 }, for given n 1 , n 2 ∈ N 0 .
In order to define a discrete-time dynamic model based on (1), variables x k+1 and ẋk+1 should be computed based on known x k , ẋk and τ τ τ k . The following discrete-time dynamic system is considered:

y k+1 = A I ∆t I 0 I y k + B(y k ) 0 ∆t M(x k ) -1 W(x k ) τ τ τ k + 0 ∆t M(x y k ) -1 g(x k ) -C(x k , ẋk )ẋ k v(y k ) (7)
or, using a more compact notation,

y k+1 = φ φ φ y (y k , τ τ τ k ) = A y k + B(y k ) τ τ τ k + v(y k ) (8)
with φ φ φ y (y k , τ τ τ k ) denoting the transition mapping that represents the discrete-time dynamic model. Details on the deduction of this model are given in Appendix A. The transition mapping in (8) can be used in the prediction step of an NMPC scheme. For a state y k ∈ R 2n , the sequence of h p ∈ N future states is predicted for a known sequence of future cable tensions τ τ τ(•) ∈ U h p k by using the transition mapping φ φ φ y recursively. More precisely,

y τ τ τ(•) (0, y k ) = y k y k+1 = y τ τ τ(•) (1, y k ) = φ φ φ y (y k , τ τ τ k ) y k+2 = y τ τ τ(•) (2, y k ) = φ φ φ y (y k+1 , τ τ τ k+1 )
. . .

y k+h p = y τ τ τ(•) (h p , y k ) = φ φ φ y y k+h p -1 , τ τ τ k+h p -1 (9)
where the state y τ τ τ(•) ( j, y k ) represents the j th term of the sequence of predicted future states y τ τ τ( Although the domain of φ φ φ y may be defined as R 2n × R m , a proper CDPR operation should be constrained to a limited set of cable tensions. For this reason, the set of feasible cable tensions U is used instead of R m . Similarly, a constrained set of platform poses and velocities should be defined. The set Y ⊂ R 2n is thus defined as the set of all admissible states y ∈ R 2n within which the controller is able to operate. The definition of Y should use tools related to the computation of the workspace of CDPRs, e.g. [START_REF] Gagliardini | Determination of a Dynamic Feasible Workspace for Cable-Driven Parallel Robots[END_REF]- [START_REF] Pott | Determination of the wrench-closure translational workspace in closed-form for cable-driven parallel robots[END_REF]. More details on this matter are discussed in Section VI-A2. The set X ⊂ R n is defined as the projection of the set Y on the space of platform

poses, i.e. X = {x ∈ R n | ∃ ẋ ∈ R n : x T ẋT T ∈ Y}.
In summary, the domain of φ φ φ y is defined as Y × U.

III. NMPC FOR CDPR POSITION TRACKING

The block diagram in Figure 3 outlines an overall position tracking control scheme integrating the NMPC strategy proposed in this paper. Based on the desired and estimated platform motions (x d (•), ẋd (•) ∈ S n and x(•), ẋ(•) ∈ S n , respectively), the NMPC scheme (block (a)) defines a set of desired cable tensions τ τ τ d (•) ∈ S m . The cable tension control in block (d) is responsible for generating these cable tensions by means of desired motor velocities qd . As described in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] and [START_REF] Santos | A Simple and Efficient Non-Model Based Cable Tension Control[END_REF], qd is computed applying a feedforward term q f f related to the platform velocity in addition to a PI feedback correction based on the error between the actual cable tensions and the desired ones. In addition, the authors showed in [START_REF] Santos | A Simple and Efficient Non-Model Based Cable Tension Control[END_REF] that this simple non model-based control scheme leads to small tension tracking errors even in the presence of friction. The actual platform pose and velocity are estimated based on the measured motor positions and velocities. Blocks (b) and (c) correspond to the kinematic modeling of the CDPR and can be implemented using well-known methods, e.g. [START_REF] Pott | On the forward kinematics of cable-driven parallel robots[END_REF]- [START_REF] Santos | A Real-Time Capable Forward Kinematics Algorithm for Cable-Driven Parallel Robots Considering Pulley Kinematics[END_REF]. Note that the NMPC in block (a) is independent of blocks (b), (c) and (d) and thus independent of cable properties (e.g. elasticity) and friction at the winches and routing pulleys.

The predicted platform velocity ẋ f f is obtained within the NMPC block applying [START_REF] Alikhani | Sliding Mode Control of a Cable-driven Robot via Double-Integrator Sliding Surface[END_REF] with τ τ τ d (k).

The NMPC scheme introduced in this section is consistent with the standard formulation of NMPC schemes without terminal conditions [START_REF] Grüne | Nonlinear model predictive control[END_REF]. In accordance with the definition of receding horizon MPC schemes [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], the proposed NMPC strategy solves an Optimal Control Problem (OCP) at each controller cycle such that a cost functional is minimized while satisfying system constraints. Since no terminal condition is used, the cost functional, denoted by J h p : N 0 × Y × U h p k → R for a given k ∈ N 0 , consists of a sum of individual costs associated to each sampling time within the prediction horizon h p ∈ N. The function that computes these individual costs is called stage cost and is denoted by : N 0 × Y × U → R. Therefore, for an instant k ∈ N 0 , actual state y k ∈ Y and sequence of cable tensions τ τ τ(•) ∈ U h p k , the cost functional is given by

J h p (k, y k , τ τ τ(•)) = h p -1 ∑ j=0 k + j, y τ τ τ(•) ( j, y k ), τ τ τ k+ j . ( 10
)
Algorithm 1 summarizes the common implementation of NMPC schemes without terminal conditions applying the notations used in this paper. Algorithm 1-(b) defines the control policy τ τ τ f b : N 0 ×Y → U, which computes the vector of desired cable tensions τ τ τ f b (k, y k ) at instant k ∈ N 0 and for state y k ∈ Y. This vector is computed based on the minimization of the cost functional J h p . The computation of J h p applies the discretetime system [START_REF] Alikhani | Sliding Mode Control of a Cable-driven Robot via Double-Integrator Sliding Surface[END_REF] in order to predict the h p future states, as in [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF]. The cost functional J h p is computed considering a given stage cost (which will be defined later in this section). The control policy τ τ τ f b is used in Algorithm 1-(a) such that the minimization of J h p is solved at each controller cycle with updated estimations of the states. As main output, the NMPC scheme defines a set of desired cable tensions τ τ τ d (k) to be applied in the time interval t ∈ [k ∆t, (k + 1) ∆t). These outputs are used as setpoints for the cable tension control in block (d) shown in Figure 3. The numerical procedure used for the minimization of

J h p is discussed in Appendix D.
The definition of a pertinent stage cost is a crucial step in the design of a stable and effective NMPC scheme. In order to meet the theoretical requirements necessary to guarantee stability and obtain appropriate performances, one should seek some particular properties of the stage cost. Postponing the formal developments on this matter to Section IV, some key properties are sketched in the next paragraphs to explain the rationale used in the design of the proposed OCP.

Typically, NMPC schemes penalize states and control inputs with respect to the desired behavior of the system. The following definition rigorously formulates such desired behavior.

y τ τ τ nd (•) (k, y 0 ) = y d (k) ∈ Y and τ τ τ nd (k) ∈ U (11)
for all k ∈ N. Moreover, each vector of τ τ τ nd (•) is considered to have minimal 2-norm such that, for every k ∈ N,

τ τ τ ∈ U | W (x d (k)) τ τ τ nd (k) = W (x d (k)) τ τ τ and τ τ τ < τ τ τ nd (k) . ( 12
)
Algorithm 1 NMPC Algorithm (a) Overall NMPC Implementation

1: Set k ← 0; 2: loop 3:
Estimate the actual states y k using a forward kinematic model;

4: Set desired cable tensions τ τ τ d (k) ← τ τ τ f b (k, y k ) valid for t ∈ [k ∆t, (k + 1) ∆t); 5:
Set k ← k + 1; 6: end loop (b) Control policy -computation of the desired cable tensions

Inputs: x T k ẋT k T = y k ∈ Y and k ∈ N 0 ; Output: τ τ τ f b ∈ U; 1: function τ τ τ f b (k, y k ) 2: Find τ τ τ * (•) ∈ U h p k that minimizes J h p k, y k , τ τ τ(•) 3:
and satisfies constraints (16b);

4: Set τ τ τ f b ← τ τ τ * k = τ τ τ * (k); 5:
return τ τ τ f b . 6: end function A desired trajectory y d (•) ∈ S 2n that does not satisfy the aforementioned conditions is called unfeasible.

As indicated in (12) each vector of τ τ τ nd (•) presents minimal 2-norm. Indeed, the minimization of the 2-norm of the vector of cable tensions is widely used in the state-of-the-art control schemes, e.g. [START_REF] Oh | Cable suspended planar robots with redundant cables: Controllers with positive tensions[END_REF]- [START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF]. This choice leads to continuous cable tension while minimizing the energy consumption [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n + 2 Cables[END_REF].

The cost functional of a typical NMPC scheme should be null for state and cable tensions equal to the desired trajectory and nominal desired cable tensions, respectively, and strictly positive if the state is not equal to the desired one. More precisely, the stage cost function should satisfy (see [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 3.3])

(k, y, τ τ τ) = 0 if y = y d (k) and τ τ τ = τ τ τ nd (k), (k, y, τ τ τ) > 0 if y = y d (k). (13)
Aiming at a stage cost satisfying ( 13), (k, y, τ τ τ) could be defined as a sum of weighted norms such as

y -y d (k) K y + τ τ τ -τ τ τ nd (k) K τ τ τ , with v K = √ v T K v,
for any vector v and positive definite matrix K. Nevertheless, in order to obtain τ τ τ nd (•), one would need to solve the tension distribution problem along the whole trajectory beforehand, using a tension distribution algorithm. The nominal wrenches to be applied on the platform should then be computed based on the nominal inverse dynamics. Therefore, these wrenches could not be updated in the controller real-time operation. Due to the incidence of disturbances and modeling uncertainties, the wrench computed with the nominal inverse dynamics can be inefficient for the tracking control. Furthermore, the desired trajectory should be feasible in the sense that there exists a sequence of cable tension distributions in U able to generate this trajectory.

Conversely, the proposed NMPC scheme is able to define in real time the optimal cable tension distribution without fixing a desired wrench. For this purpose, the following stage cost is proposed:

(k, y, τ τ τ) = y -y d (k) 2 K y + τ τ τ -τ τ τ wo (x y , τ τ τ) 2 K τ τ τ , (14)
where

K y = diag(k y ) and K τ τ τ = diag(k τ τ τ ), with constant weighting vectors k y ∈ R 2n + and k τ τ τ ∈ R m + . The term τ τ τ -τ τ τ wo (x y , τ τ τ) 2
K τ τ τ represents the Wrench Equivalent Optimality (WEO) with respect to x y and τ τ τ. The WEO is a nonnegative measure that evaluates whether the wrench generated by τ τ τ can be generated by an alternative set of cable tensions with smaller 2-norm at the pose x y . More precisely, for given x ∈ X and τ τ τ ∈ U, the vector-valued function τ τ τ wo : X × U → U denotes the wrench equivalent optimal (WE-optimal) cable tensions, defined as

τ τ τ wo (x, τ τ τ) = arg min τ τ τ τ τ τ 2 (15a) s. t. W(x) τ τ τ = W(x) τ τ τ (15b) τ τ τ min τ τ τ τ τ τ max (15c)
In words, τ τ τ wo (x, τ τ τ) is the vector of cable tensions with minimal 2-norm that is able to generate the same wrench as τ τ τ at the pose x. Therefore, for given k ∈ N, y ∈ Y and τ τ τ ∈ U, the term τ τ ττ τ τ wo (x y , τ τ τ) 2 K τ τ τ = 0 iff τ τ τ has minimal 2-norm (τ τ τ = τ τ τ wo (x y , τ τ τ)). Otherwise, τ τ ττ τ τ wo (x y , τ τ τ) 2 K τ τ τ > 0, indicating that there exists a τ τ τ wo (x y , τ τ τ) ∈ U generating the same wrench (thanks to (15b)) with τ τ τ wo (x y , τ τ τ) < τ τ τ .

Accordingly, for given time instant k ∈ N 0 and state y k ∈ Y, the OCP solved in Algorithm 1-(b), which corresponds to the minimization of J h p (k, y k , τ τ τ(•)) with respect to τ τ τ(•), can be written as min

τ τ τ(•)∈U hp k h p -1 ∑ j=0 k + j, y τ τ τ(•) ( j, y k ) , τ τ τ k+ j (16a) s. t. y τ τ τ(•) j, y k ∈ Y ∀ j ∈ N 1,h p (16b)
with the stage cost defined in [START_REF] Oh | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF]. Note that N 1,h p is considered in (16b) while the summation in (16a) is taken from j = 0 to h p -1 since ( 9) is used in order to predict future states up to

y k+h p = y τ τ τ(•) (h p , y k ) = φ φ φ y y k+h p -1 , τ τ τ k+h p -1 and (16b)
ensures that all these predicted future states are feasible.

According to the definition of the proposed stage cost in [START_REF] Oh | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF], it is worth highlighting that ( 15) is an optimization problem nested within the main problem [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF]. More precisely, the proposed OCP can be classified as a bilevel optimization problem (see [START_REF] Dempe | Foundations of bilevel programming[END_REF], [START_REF] Sun | Fast UAV trajectory optimization using bilevel optimization with analytical gradients[END_REF]), in which the lower level minimization [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF] is embedded within the upper level [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF]. In practice, the solution of the upper level problem ( 16) cannot be decoupled from the solution of the lower level [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF] and each iteration on the numerical solution of ( 16) should consider updated values of the WE-optimal cable tensions computed according to [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF]. For further details on the numerical solution of ( 16), refer to Appendix D.

Note also that, as a consequence of ( 12) and ( 15),

τ τ τ wo (x d (k), τ τ τ nd (k))-τ τ τ nd (k) K τ τ τ + y d (k) -y d (k) K y = k, y d (k), τ τ τ nd (k) = 0 ∀ k ∈ N (17)
and conditions ( 13) are satisfied with the stage cost [START_REF] Oh | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF]. Thereby, the minimal 2-norm tension distribution is performed when solving the OCP ( 16) , i.e., the cable tension distribution is implicitly performed in real time within the OCP.

In order to exemplify the advantage of this approach, consider k ∈ N, y k ∈ Y and a sequence of cable tensions τ τ τ(•) ∈ U h p k that generates an optimal trajectory y * τ τ τ(•) (•, y k ) = {y * k+1 , . . . , y * k+h p } in the sense that the tracking error

h p ∑ j=1 y * k+ j -y d (k + j) K y (18)
is minimal. A sequence of WE-optimal cable tensions τ τ τ * (•) defined such that

τ τ τ * k+ j = τ τ τ wo (x y * k+ j , τ τ τ k+ j ), ∀ j ∈ N 0,h p -1 (19)
presents minimal 2-norm and generates the very same optimal trajectory

y * τ τ τ(•) (•, y k ). In accordance with Definition 1, if y k = y d (k) and y d (•) is feasible, τ τ τ * k+ j = τ τ τ nd (k + j) for all j ∈ N 0,h p -1 , J(k, y k , τ τ τ * (•)) = 0 and
the resulting tracking error is null. The use of the WEO in the stage cost ( 14) is thus consistent with the minimization of the tracking errors in the sense that WE-optimal cable tensions are able to generate the desired trajectory with null tracking errors. Furthermore, in case of disturbances, modeling uncertainties and unfeasible y d (•), the proposed OCP formulation is still able to find an optimal trajectory y * τ τ τ * (•) (•, y k ) with minimal tracking error and WE-optimal cable tensions τ τ τ * (•).

IV. STABILITY ANALYSIS

This section analyzes the stability of the closed-loop system

y k+1 = φ φ φ y y k , τ τ τ f b (k, y k ) (20)
obtained with the discrete-time system (8) and the NMPC control policy τ τ τ f b introduced in Algorithm 1-(b) and [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF].

To this end, some comparison functions are used. In accordance with commonly used notations [START_REF] Grüne | Nonlinear model predictive control[END_REF], [START_REF] Khalil | Nonlinear control[END_REF], the following classes of functions are considered

K := {α : R + 0 → R + 0 |α continuous, strictly increasing and α(0) = 0}; L := {δ : R + 0 → R + 0 | δ continuous, strictly decreasing and lim t→∞ δ (t) = 0}; K L := {β : R + 0 × R + 0 → R + 0 | β (r, •) ∈ L and β (•,t) ∈ K } K ∞ := {α ∈ K | α unbounded}, where R + 0 = {r ∈ R | r 0}.
Furthermore, for a given sequence s 0 (•) ∈ S n s and instant k ∈ N 0 , special attention will be devoted to a truncated sequence s k (•) = {s 0 (k), s 0 (k +1), . . . }. Accordingly, for given n s ∈ N and n i ∈ N 0 , the following set of sequences of vectors is introduced:

S n s n i = {s(•) : N n i ,∞ → R n s }. ( 21
)
The stability analysis presented in this section deduces sufficient conditions under which the closed-loop system [START_REF] Li | Tracking control of fully-constrained cabledriven parallel robots using adaptive dynamic programming[END_REF] is uniformly asymptotically stable according to the following definition.

y j+1 = φ φ φ y y j , τ τ τ f b ( j, y j ) , ∀ j ∈ N ( 22
)
is called uniformly asymptotically stable to

y d (•) on Y if, for each y k ∈ Y and k ∈ N, there exists a function β ∈ K L such that y f b (•) ∈ S 2n
k defined according to

y f b k = y k y f b i+1 = φ φ φ y y f b i , τ τ τ f b (i, y f b i ) , ∀ i ∈ N k,∞ (23)
satisfies the following relation

y f b k+ j -y d (k + j) β ( y k -y d (k) , j), ∀ j ∈ N. (24)
In words, [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n + 2 Cables[END_REF] constructs the sequence of states y f b (•) generated by the closed-loop ( 22) "departing" from y k at sampling time k. The closed-loop [START_REF] Oh | Cable suspended planar robots with redundant cables: Controllers with positive tensions[END_REF] takes as feedback control policy τ τ τ f b defined in Algorithm 1. The definition of uniform asymptotic stability presented in Definition 2 is based on an upper bound of the error y f b k+ j -y d (k + j) . Note that Definition 2 implicitly assumes τ τ τ k = τ τ τ d (k), ∀ k ∈ N, i.e. the cable tension control loop (d) in Figure 3 is able to keep a null error.

The upper bound used to limit the tracking error is established in terms of the function β ∈ K L . An illustrative example of a function belonging to this class is presented in Figure 4. In accordance with (IV), for a positive real constant r 0 , the illustrative function β ex (r 0 ,t) = r 3 0 /t tends to zero for t → ∞ and is strictly decreasing. Since the second argument of β in [START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF] represents the time, the tracking error asymptotically converges to zero. In contrast, β is strictly increasing with respect to the first argument. Therefore, for increased initial errors y k -y d (k) , the function of time

β ( y k -y d (k) , •) ∈ L is increased as well.
Definition 2 is used in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Definition 2.16] for the analysis of NMPC schemes and is closely related to [51, Definition 4.2], which is more often used in general nonlinear control. As a matter of fact, the analysis presented in this section is based on [START_REF] Grüne | Nonlinear model predictive control[END_REF] and the notations used are similar to those proposed in this reference.

Since the OCP [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF] does not use stabilizing terminal conditions, tools designed to analyze NMPC schemes without terminal conditions should be used in the stability analysis of the closed-loop [START_REF] Li | Tracking control of fully-constrained cabledriven parallel robots using adaptive dynamic programming[END_REF]. The minimal time-varying stage cost denoted by * : N × Y → R is fundamental in this context. This function is defined as Summarizing some of the results presented in [43, Chapter 6], the following theorem addressing the stability of NMPC schemes without terminal conditions can be stated.

* (k, y) = min τ τ τ∈U (k, y, τ τ τ). ( 25
)
Theorem 1. Consider the NMPC Algorithm 1 with prediction horizon h p ∈ N and minimal time-varying stage cost satisfying

α 1 ( y -y d (k) ) * (k, y) α 2 ( y -y d (k) ) ∀ k ∈ N and y ∈ Y, (26)
for suitable α 1 , α 2 ∈ K ∞ . Suppose that, for all y k ∈ Y, there exist a feasible τ τ τ e (•) ∈ U m k , real C < ∞ and σ ∈ (0, 1)

satisfying k + j, y τ τ τ e (•) ( j, y k ), τ τ τ e k+ j C σ j * (k, y k ) , (27)
for all k, j ∈ N. Then, the nominal closed-loop (20) is uniformly asymptotically stable on Y provided that h p is sufficiently large.

A draft of the proof of Theorem 1 is presented in Appendix B. This theorem describes sufficient conditions to obtain a stable closed-loop system with the proposed NMPC for h p sufficiently large.

Conditions (26) will be quite straightforwardly deduced at the end of this section based on [START_REF] Zeinali | Design and Application of Chattering-Free Sliding Mode Controller to Cable-Driven Parallel Robot Manipulator: Theory and Experiment[END_REF]. Nevertheless, inequality [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF] represents an important restriction to the properties of the stage cost. This inequality imposes, for each discrete instant k ∈ N and state y k ∈ Y, the existence of a feasible sequence τ τ τ e (•) that leads to a stage cost k + j, y τ τ τ e (•) ( j, y k ), τ τ τ e j converging exponentially to zero in time. The exponential convergence is a result of the term C σ j with σ ∈ (0, 1). Note that the upper bound function C σ j * (k, y k ) is proportional to the initial minimal stage cost * (k, y k ).

It is important to highlight that Theorem 1 is based on the existence of any feasible sequence τ τ τ e (•) that satisfies [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF]. The control policy τ τ τ f b is not considered explicitly in this theorem. In contrast to many control architectures, constrained MPC schemes generally do not have an analytical expression of the control policy. For this reason, the closed-loop system stability cannot be studied explicitly, but, instead, by means of suitable properties of the stage cost. In this section, these suitable property are given by ( 26)- [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF].

Since one of the terms in the stage cost ( 14) is the error yy d 2 K y , this error needs to be bounded. To this end, the following assumption establishes important properties of the controlled system in order to attain [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF].

For each y k ∈ Y and k ∈ N, there exist

τ τ τ(•) ∈ S m k , real C 1 < ∞ and σ 1 ∈ (0, 1) such that e y k + j, y τ τ τ(•) ( j, y k ) C 1 σ j 1 e y (k, y k ) (29a) τ τ τ k+ j ∈ U (29b)
for all j ∈ N 0 .

If Assumption 1 is true, then, for each state y k ∈ Y and instant k ∈ N, there exists a feasible sequence of cable tensions τ τ τ(•) that generates a trajectory in which the error y τ τ τ(•) ( j, y k ) -y d (k + j) K y exponentially converges to zero. Moreover, the upper bound function in (29a) is proportional to the initial errors e y (k, y k ). Note that, for given k and y k , the error e y (k, y k ) is constant with respect to the discrete time j. It is easy to show that this assumption is satisfied with a typical control scheme, as detailed in Appendix C. It is interesting to reiterate that the present stability analysis relies on the existence of any control sequence satisfying some particular properties. The existence of a control policy satisfying (29a) will be used in order to deduce sufficient conditions to satisfy [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF], which plays an important role in the stability analysis of [START_REF] Li | Tracking control of fully-constrained cabledriven parallel robots using adaptive dynamic programming[END_REF].

As important as (29a) imposing the exponential convergence of error e y (k, •), (29b) requires that the sequence τ τ τ(•) is feasible. Considering for instance one of the strategies proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF], this means that the wrench defined in the feedback motion control policy is feasible for any pose and velocity within the state constraint set Y. Clearly, the definition of the robot workspace and the set Y ⊂ R 2n play crucial roles in Assumption 1. The study of the feasibility of a given wrench considering different poses and velocities is addressed in several works (e.g. [START_REF] Gagliardini | Determination of a Dynamic Feasible Workspace for Cable-Driven Parallel Robots[END_REF]- [START_REF] Pott | Determination of the wrench-closure translational workspace in closed-form for cable-driven parallel robots[END_REF]) and is out of the scope of this paper. Let us also point out that Assumption (29b) is often implicitly taken. Studies such as [START_REF]Dynamic modeling and control of parallel robots with elastic cables: Singular perturbation approach[END_REF], [START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF], [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF], [START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF], [START_REF] Shang | Dual-space adaptive synchronization control of redundantly-actuated cable-driven parallel robots[END_REF] analyze the corresponding closed-loop systems considering that the feedback loop does not lead to unfeasible cable tensions and, therefore, satisfies (29b).

Since the term e y (k, y) = yy d (k) 2 K y in the expression (14) of the stage cost (k, y, τ τ τ) is independent of τ τ τ, the minimal time-varying stage cost is given by * (k, y) = e y (k, y) + min

τ τ τ∈U ( τ τ τ -τ τ τ wo (x y , τ τ τ) 2 K τ τ τ ). (30)
Lemma 2 will show that the minimum min τ τ τ∈U ( τ τ ττ τ τ wo (x y , τ τ τ) K τ τ τ ) is known and (30) can thus be simplified. Before stating and proving Lemma 2, a preparatory lemma is necessary.

Lemma 1. Consider τ τ τ a , τ τ τ b ∈ U and x ∈ X. If W(x) τ τ τ a = W(x) τ τ τ b , then τ τ τ wo (x, τ τ τ a ) = τ τ τ wo (x, τ τ τ b ). (31)
Proof. First, the existence of a τ τ τ wo (x, τ τ τ) for all x ∈ R n , τ τ τ ∈ U should be proved. The objective function τ τ τ 2 in (15a) is quadratic and strictly convex. The substitution of τ τ τ = τ τ τ ∈ U in the constraints (15b)-(15c) shows that τ τ τ itself is an element of the set defined by these constraints. Hence, this set is convex and non-empty. Therefore, ( 15) is a feasible strictly convex inequality constrained Quadratic Programming (QP) problem that possesses a global minimum. Denoting f ∈ R n the wrench such that f = W(x) τ τ τ a , Definition [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF] indicates that both τ τ τ wo (x, τ τ τ a ) and τ τ τ wo (x, τ τ τ b ) are obtained with

τ τ τ wo (x, τ τ τ a ) = τ τ τ wo (x, τ τ τ b ) = arg min τ τ τ τ τ τ 2 2 (32a) s. t. W(x) τ τ τ = f (32b) τ τ τ min τ τ τ τ τ τ max (32c)
and, therefore, (31) is true.

Lemma 1 shows that if two vectors of cable tensions generate the same wrench in a given pose, they lead to the same WE-optimal tension distribution. The following corollary is based on this result.

Corollary 1. For every τ τ τ a , τ τ τ b ∈ U and x ∈ X, if τ τ τ b = τ τ τ wo (x, τ τ τ a ), (33)
then τ τ τ b = τ τ τ wo (x, τ τ τ b ). (34)
Proof. Taking τ τ τ a , τ τ τ b and x from the corollary statement such that τ τ τ b = τ τ τ wo (x, τ τ τ a ), constraints (15b)-(15c) in the definition of τ τ τ wo implies that W(x) τ τ τ a = W(x) τ τ τ b and τ τ τ b ∈ U. Therefore, applying Lemma 1, τ τ τ wo (x, τ τ τ a ) = τ τ τ wo (x, τ τ τ b ), and thanks to [START_REF] Rawlings | Model predictive control: Theory, Computation, and Design[END_REF], this leads to τ τ τ b = τ τ τ wo (x, τ τ τ b ).

Lemma 2. For every y ∈ Y, τ τ τ ∈ U and k ∈ N,

min τ τ τ∈U ( τ τ τ -τ τ τ wo (x y , τ τ τ) 2 K τ τ τ ) = 0, ( 35
)
and the minimal time-varying stage cost is given by

* (k, y) = y -y d (k) 2 K y = e y (k, y). ( 36
)
Proof. For every y ∈ Y and τ τ τ ∈ U, one may set τ τ τ out = τ τ τ wo (x y , τ τ τ), and thanks to the constraints (15c), τ τ τ out ∈ U. Moreover, we have τ τ τ outτ τ τ wo (x y , τ τ τ) 2 K τ τ τ = 0. Since the expression τ τ ττ τ τ wo (x, τ τ τ) 2 K τ τ τ = 0 is non-negative, this proves [START_REF] Sun | Disturbance Rejection MPC for Tracking of Wheeled Mobile Robot[END_REF]. Substituting [START_REF] Sun | Disturbance Rejection MPC for Tracking of Wheeled Mobile Robot[END_REF] in [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF], equation ( 36) is obtained.

As discussed in Section III, τ τ τ wo (x, τ τ τ) is the vector of cable tensions with minimal 2-norm that generates the same wrench as τ τ τ in pose x. Loosely speaking, the influence of τ τ τ wo (x, τ τ τ) and τ τ τ on the system dynamics should be equivalent. This equivalence is described in a more precise manner in the following assumption. (37a)

W(x y ) τ τ τ b = W(x y ) τ τ τ a (37b)
then, the following relation also holds:

y + = φ φ φ y (y, τ τ τ b ) = φ φ φ y (y, τ τ τ a ). (38)
It is easy to demonstrate that Assumption 2 is satisfied with the transition mapping [START_REF] Williams | Planar translational cabledirect-driven robots[END_REF]. It can be noted that this assumption is implicitly taken if the stability of a digital control system is analyzed considering the robot controller as a continuous-time system, as in [START_REF]Dynamic modeling and control of parallel robots with elastic cables: Singular perturbation approach[END_REF], [START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF], [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF], [START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF], [START_REF] Shang | Dual-space adaptive synchronization control of redundantly-actuated cable-driven parallel robots[END_REF].

Based on Assumption 2, for a given initial state vector, two sequences of cable tensions that generate identical wrenches at the corresponding poses result in identical trajectories. Therefore, any trajectory performed with an arbitrary sequence of feasible cable tensions can also be generated with an alternative sequence consisting of WE-optimal tensions generating the same wrenches along the trajectory. With the latter cable tension sequence, the stage cost at each instant along the trajectory is equal to the minimal time-varying cost * . This assertion is rigorously formulated in Lemma 3.

Lemma 3. For each y k ∈ Y, y(•) ∈ S 2n
k and τ τ τ(•) ∈ S m k such that, for all j ∈ N 0 ,

y k+ j = y τ τ τ(•) ( j, y k ) ∈ Y (39a) τ τ τ k+ j ∈ U, (39b)
if Assumption 2 holds, there exists a τ τ τ * (•) ∈ S m k that satisfies

y τ τ τ(•) ( j, y k ) = y τ τ τ * (•) ( j, y k ) = y k+ j , (40a)
τ τ τ wo (x y k+ j , τ τ τ * k+ j ) = τ τ τ * k+ j ∈ U (40b)

(k + j, y k+ j , τ τ τ * k+ j ) = e y (k + j, y k+ j ) = * (k + j, y k+ j ) (40c) for all j ∈ N 0 .

Proof. For y k , y(•) and τ τ τ(•) stated in the lemma, define τ τ τ * (•) ∈ S m k and y * (•) ∈ S 2n k as τ τ τ * k+ j = τ τ τ wo (x y k+ j , τ τ τ k+ j ),

y * k+ j = y τ τ τ * (•) ( j, y k ) (41)
for all j ∈ N 0 . Due to the constraint (15b),

W(x y k ) τ τ τ k = W(x y k ) τ τ τ * k ( 42
)
and, according to [START_REF] Faulwasser | Nonlinear Model Predictive Control for Constrained Output Path Following[END_REF] and ( 9),

y * k+1 = y τ τ τ * (•) (1, y k ) = φ φ φ y (y k , τ τ τ * k ) = y k+1 . (43)
Moreover, ( 41)-( 42) and Lemma 1 imply that

τ τ τ wo (x y k , τ τ τ * k ) = τ τ τ wo (x y k , τ τ τ k ) = τ τ τ * k . (44)
By induction, can obtain

y * k+2 = φ φ φ y (y k+1 , τ τ τ * k+1 ) = y k+2 y * k+3 = φ φ φ y (y k+2 , τ τ τ * k+2 ) = y k+3 . . . (45)
and τ τ τ wo (x y k+1 , τ τ τ * k+1 ) = τ τ τ wo (x y k+1 ,τ τ τ k+1 ) = τ τ τ * k+1 τ τ τ wo (x y k+2 , τ τ τ * k+2 ) = τ τ τ wo (x y k+2 ,τ τ τ k+2 ) = τ τ τ * k+2 . . . [START_REF] Barreto | Design and Implementation of Model-Predictive Control With Friction Compensation on an Omnidirectional Mobile Robot[END_REF] which are equivalent to (40a) and (40b). Equation (40b) and Lemma 2 imply that (k + j, y k+ j , τ τ τ * k+ j ) = e y (k + j, y k+ j ) + τ τ τ * k+ jτ τ τ wo (x y k+ j , τ τ τ * k+ j ) K τ τ τ = e y (k + j, y k+ j ) = * (k + j, y k+ j ), [START_REF] Kayacan | Robust Trajectory Tracking Error Model-Based Predictive Control for Unmanned Ground Vehicles[END_REF] as in (40c).

Finally, the main results on the stability of the NMPC scheme proposed in Section III are summarized in the following theorem.

Theorem 2. Consider the NMPC Algorithm 1 with prediction horizon h p ∈ N and feasible desired trajectory y d (•) ∈ S 2n . Consider also that Assumptions 1 and 2 hold. Then, the nominal closed-loop system (20) is uniformly asymptotically stable on Y provided that h p is sufficiently large.

Proof. The proof consists in analyzing the conditions presented in Theorem 1. First, it is necessary to prove that there exist

α 1 , α 2 ∈ K ∞ such that α 1 ( y -y d (k) ) * (k, y) α 2 ( y -y d (k) ) ∀ k ∈ N and y ∈ Y. (48)
According to Lemma 2, * (k, y) = yy d (k) 2 K y such that the following inequalities hold for all k ∈ N:

K min y y -y d (k) 2 * (k, y) = y -y d (k) 2 K y K max y y -y d (k) 2 . ( 49
)
where K min y and K max y are the minimum and maximum singular value of K y , respectively.

Defining α 1 , α 2 ∈ K ∞ by α 1 (r) = K min y r 2 and α 2 (r) = K max y r 2 , (49) can be written as

α 1 ( y -y d (k) ) * (k, y) α 2 ( y -y d (k) ) ∀ k ∈ N, (50)
which proves [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF].

In addition, ( 27) should be proved. If Assumption 1 holds, then for each y k ∈ Y and k ∈ N, there exist τ τ τ(•) ∈ S m k , real C 1 < ∞ and σ 1 ∈ (0, 1) such that

e y k + j, y τ τ τ(•) ( j, y k ) C 1 σ j 1 e y (k, y k ) (51a) τ τ τ k+ j ∈ U (51b)
for all j ∈ N 0 . Applying Lemma 3 with τ τ τ(•) = τ τ τ(•), the obtained sequence τ τ τ * (•) leads to the same trajectory (as implied by (40a)) and, therefore, (51a) remains valid. In addition, (40c) implies that k + j, y τ τ τ * (•) ( j, y k ), τ τ τ * k+ j = e y (k + j, y τ τ τ * (•) ( j, y k )) (52a) where (51a) and ( 36) were used in order to obtain (52b) and (52c), respectively. Inequality ( 52) proves ( 27), finalizing the stability analysis.

C 1 σ j 1 e y (k, y k ) (52b) = C 1 σ j 1 * (k, y k ), (52c)

V. NUMERICAL SIMULATIONS

The performance of the proposed NMPC scheme is compared through numerical simulations to the LMPC introduced in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF]. Typical pick-and-place trajectories are used and the dimensions of CDPR prototype HRPCable are considered (see Figures 1 and2). The technical characteristics of this CDPR prototype are discussed more in detail in the next section.

Instead of the stage cost (k, y, τ τ τ) as in [START_REF] Oh | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF], the LMPC of [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] is based on the cost function LMPC :

N 0 × Y × U × U → R defined as LMPC (k, y, τ τ τ k , τ τ τ k-1 ) = y -y d (k) 2 K 1 + c 2 τ τ τ k 2 + + c 3 τ τ τ k -τ τ τ k-1 2 , ( 53
)
where c 2 , c 3 are positive real constants, K 1 = diag(k 1 ) with k 1 ∈ R 2n + and τ τ τ kτ τ τ k-1 represents the 2-norm of the variation of the cable tensions between sampling times k -1 and k. More details on LMPC schemes for the position tracking of CDPRs are presented in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], [START_REF] Santos | Model Predictive Control Applied to Large-Dimension Cable-Driven Parallel Robots[END_REF].

Optimal controller gains for both the LMPC from [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] and the NMPC proposed in Section III were obtained through Particle Swarm Optimization (PSO). These optimizations were based on the trajectories depicted in Figures 5 and6. The platform orientation is defined by three Euler angles ψ 1 , ψ 2 , ψ 3 . More precisely, the rotation matrix used in ( 63) is given by

R(ψ ψ ψ) = R z (ψ 3 ) R y (ψ 2 ) R x (ψ 1 ), (54)
with R x , R y , R z : R → R 3×3 the elementary rotation matrices around x, y and z axes, respectively. The cable tension limits are taken as τ τ τ min = 100 N and τ τ τ max = 400 N.

Considering the desired trajectory shown in Figure 5, the simulation of the closed-loop system obtained with the LMPC scheme using a given set of gains {k 1 , c 2 , c 3 } leads to translation errors TE(k) and orientation errors OE(k) for any instant k within the time interval corresponding to the trajectory motion. Accordingly, the optimal LMPC gains were obtained by means of the minimization of a weighted sum of the maximum translation and orientation errors, i.e.:

{k * 1 , c * 2 , c * 3 } = arg min {k 1 , c 2 , c 2 } w t max k TE(k) + w o max k OE(k) ,
for positive constant weights w t and w o . Optimal NMPC gains {k * y , k * τ τ τ } were obtained with an equivalent procedure and using the same weights w t and w o . The sampling period of both control schemes was considered as ∆t = 22 ms. The NMPC prediction horizon was defined as h p = 4 while the LMPC prediction horizon was taken as h p,LMPC = 8. A NMPC prediction horizon smaller than the LMPC one was chosen in order to compensate for the increased complexity of the optimal control problem addressed within the NMPC strategy. These parameter values are in accordance with the experimental tests presented in Section VI.

Figure 6 compares the errors obtained with the LMPC and NMPC schemes using the corresponding optimal control gains. Relatively small errors are obtained with the LMPC scheme since TE(k) < 8 mm and OE(k) < 0.15 degrees for all k in the trajectory time interval. However, the NMPC control strategy proposed in this paper performs much better as it is able to maintain errors virtually equal to zero, namely, TE(k) < 0.03 mm and OE(k) < 2 × 10 -4 degrees for all k in the trajectory time interval. These results are consistent with the discussion presented in Section III. Namely, while the LMPC minimizes the cable tensions in detriment of the tracking errors, the proposed NMPC scheme is able to compute an optimal cable tension distribution without degrading the tracking errors.

In order to better illustrate the contrast between these two control schemes, the values of the maximum translation and orientation errors were evaluated for different controller gains. More precisely, Figure 7 compares the errors obtained with sets of random controller gain values uniformly distributed within the intervals given by 0 for the NMPC, and

.5 × k * y k y 1.5 × k * y 0.5 × k * τ τ τ k τ τ τ 1.5 × k * τ τ τ (55)
     0.5 × k * 1 k 1 1.5 × k * 1 0.5 × c * 2 c 2 1.5 × c * 2 0.5 × c * 3 c 3 1.5 × c * 3 ( 56
)
for the LMPC. The logarithmic plot chart depicted in Figure 7 indicates that the sets of NMPC controller gains {k y , k τ τ τ } within the intervals (55) lead to TE(k) < 0.1 mm and OE(k) < 10 -3 degrees. In contrast, as shown in Figure 7, the 1300 different LMPC controller gains sampled in intervals ( 56) yield max k TE(k) > 7.0 mm and max k OE(k) > 0.12 degrees.

VI. EXPERIMENTAL VALIDATION

The NMPC scheme proposed in Section III was implemented on the HRPCable prototype (Figure 1) leading to the experimental results detailed in the present section. HRP-Cable is a 6-DoF CDPR fully-constrained by eight cables. Its dimensions are shown in Figures 2 and9. The robot is essentially built with Beckhoff components typically used in industry. Each winch is equipped with a servo drive AX5112, a motor AM8061 and gear train AG2210. An industrial PC C6920 equipped with an Intel 2.4GHz i7 core processor communicates with the servo drives through EtherCAT protocol. Cable tensions are measured with load pins Sensy 5300 1T SL positioned in the axes of the routing pulleys. The platform is a cube with length equal to 1 m and total mass m p = 23 kg.

Apart from the mechanical components, the presented experimental results were obtained using the industrial software TwinCAT [START_REF] Beckhoff | TwinCAT 3 -C++ Manual[END_REF]. This software is not compatible with libraries commonly used for matrix manipulation. For this reason, the proposed NMPC scheme was developed from scratch. The environment in which the controller was implemented proves the applicability of the proposed solutions in the industry. Appendix D presents details on the numerical implementation of the control algorithm introduced in Section III. The results presented in this section also use the kinematic model and cable tension control (blocks (b), (c) and (d) in Figure 3) introduced in [START_REF] Santos | A Real-Time Capable Forward Kinematics Algorithm for Cable-Driven Parallel Robots Considering Pulley Kinematics[END_REF] and [START_REF] Santos | A Simple and Efficient Non-Model Based Cable Tension Control[END_REF], respectively.

Similarly to the previous section, the performance of the proposed NMPC scheme is compared to the LMPC introduced in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] using typical pick-and-place trajectories. These control schemes are also compared with respect to their robustness against payload uncertainties.

In addition, an unfeasible desired trajectory is used in order to show that the proposed NMPC scheme can operate on the boundaries of the robot WFW without failure. While the proposed NMPC scheme is able to track a trajectory as close as possible to an unfeasible desired trajectory, a stateof-the-art control scheme leads to a sudden interruption of the robot operation. The results obtained with the proposed NMPC scheme were summarized in a video1 that has been attached to the present paper.

Taking the optimal gains obtained in the previous section as initial estimation, the NMPC control parameters used in the experimental tests presented in this section were tuned by trial and error. Defining

k x , k ẋ ∈ R n + as subvectors of k y such that k y = k T x k T ẋ T
, the used controller parameters can be summarized as follows: As in the previous section, the cycle time of the cable tension control (block (d) in Figure 3) is 2 ms while the one used for the NMPC controller is ∆t = 22 ms. The sampling period of the cable tension control is mainly limited by the response time of the velocity control loop (which is managed at the servo drives) while the cycle time ∆t = 22 ms of the NMPC control loop is limited by the time required to solve the OCP ( 16) in real time.

K y = diag( k T x k T ẋ T ) (57a

A. Pick-and-place trajectories

The desired trajectories representing typical pick-and-place tasks are defined using fifth-degree polynomials based on the CDPR mobile platform poses depicted in Figure 9. Two scenarios of maximum cable tension are considered: (i) τ τ τ max = 400 N and (ii) τ τ τ max = 250 N. These scenarios are addressed respectively in Sections VI-A1 and VI-A2 and both cases consider τ τ τ min = 100 N.

1) Nominal cable tension limits: Figure 8 shows the experimental results obtained using the proposed NMPC scheme taking cable tension limits τ τ τ min = 100 N and τ τ τ max = 400 N. Tracking errors consist of the Translation Errors (TE) and Orientation Errors (OE). Figure 8 also depicts the cable tensions measured using the aforementioned load pins.

The same trajectory was performed with the LMPC proposed in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF]. Figure 10 compares the corresponding results. In accordance with the goals defined in the present paper, the comparative results presented in Figure 10 indicate that the proposed NMPC scheme leads to significantly better tracking errors than the LMPC approach. Indeed, Table I shows that substantial improvements were obtained considering both the RMS of the tracking errors and the maximum errors. These errors are based on the platform pose computed by the forward kinematics (block (b) in Figure 3). y (m) Experiments with τ τ τ max = 250 N. In the charts presenting the platform position and orientation, the dashed and solid lines represent the desired and actual trajectories, respectively. 2) Reduced maximum cable tension: The authors' previous study [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] highlighted that the main contribution of the corresponding LMPC scheme lies in its capability to handle cable tension limits explicitly. Thanks to this characteristic, this LMPC is not prone to the feasibility issues discussed in Section I. In order to verify this property, an unfeasible desired trajectory was used in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF]. The same behavior should be validated for the NMPC scheme proposed in the present paper.

0 3 2 x (m) 1 0 -1 -1 -2 -3 -4
Therefore, a second scenario is considered taking a maximum cable tension reduced to τ τ τ max = 250 N. The same path depicted in Figure 9 is also used in this case with reduced maximum tension. The obtained results are shown in Figure 11. Although the desired trajectory escapes from the WFW defined by this reduced maximum tension, the proposed control scheme is able to track a trajectory as close as possible to the desired one without violating the cable tension limits.

Therefore, as demonstrated in this experiment, the main advantage of the proposed NMPC scheme with respect to state-of-the-art strategies such as [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF] lies in its capability to operate on the robot WFW boundaries without failure. The Experimental results obtained with the control scheme proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] with the reduced maximum cable tension τ τ τ max = 250 N. The robot operation is terminated at t ≈ 10.5 s, when the cable tensions become unfeasible.

control strategies proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]- [START_REF] Alp | Cable suspended robots: design, planning and control[END_REF] combined with some of the tension distribution algorithms in [START_REF] Ueland | Optimal Force Allocation for Overconstrained Cable-Driven Parallel Robots: Continuously Differentiable Solutions With Assessment of Computational Efficiency[END_REF]- [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF] would fail to fulfill this objective. As an example, Figure 12 presents the experimental results obtained using the control scheme proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]. As soon as the desired trajectory reaches unfeasible poses, the robot operation is terminated because the controller requires the tension distribution module to compute an unfeasible wrench.

It is worth noting that the results presented in this section consider the state constraint set Y equal to R 2n . Hence, only the cable tension constraint set U, i.e., the WFW, limits the operation range of the NMPC scheme. This exempts the control designer from the potentially arduous task of defining the set Y. Moreover, the solution of ( 16) is also facilitated by the reduced number of constraints.

In addition, one may note that the velocities of the trajectory depicted in Figure 11 are lower than those of the original trajectory shown in Figure 8. Lower velocities were used in order to reduce platform oscillations when the desired trajectory escapes the WFW. Alternatively, results of [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF] show that the addition to the stage cost of a term τ τ τ kτ τ τ k-1 penalizing cable tension variations is an efficient method to reduce these oscillations. Nevertheless, nontrivial modifications of the stability analysis of Section IV would be necessary in this case. Accordingly, for the sake of consistency, the results presented in this section do not apply a penalization of τ τ τ kτ τ τ k-1 .

B. Robustness against payload uncertainties

CDPRs are often used in order to handle relatively heavy payloads. Typically, the weight of such payloads are not precisely known. Therefore, the tracking accuracy under the incidence of payload uncertainties is a crucial aspect to be considered in the analysis of CDPR position tracking control. Experimental results obtained with the NMPC scheme in the experiment evaluating the robustness against payload uncertainties. Accordingly, in addition to the presented pick-and-place trajectories, the robustness against payload uncertainties of the proposed NMPC scheme was evaluated using the experiment illustrated in Figure 13. An additional payload of 15.1 kg was used, representing 65.6% of the mobile platform mass.

The obtained results are depicted in Figure 14. The same procedure was performed using the LMPC proposed in [START_REF] Santos | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF], leading to the comparison of Figure 15. Table II shows that, as for the results presented in Section VI-A1, substantially better tracking errors were obtained with the NMPC scheme.

C. High-speed trajectory

As discussed in Section III, the proposed NMPC scheme considers that the cable tension control (block (d) in Figure 3) is able to maintain null errors all along the robot motion. This assumption combined to the long sampling period (22 ms) of the NMPC control loop could deteriorate the performance obtained for trajectories presenting high velocities and accelerations.

In order to evaluate whether the NMPC scheme is vulnerable to this issue, the trajectory depicted in Figure 16 was performed. The trajectory consists of two back and forth straightline motions along the x direction with extreme poses 1.66 m away from each other. Maximum velocities and accelerations were equal to 0.915 m/s and 1.18 m/s 2 , respectively. Regarding the positioning errors, the NMPC scheme led to maximum TE and OE equal to 2.06 mm and 0.138 • , respectively. In contrast, the LMPC scheme results in maximum TE and OE of 7.32 mm and 0.210 • . Hence, the NMPC scheme still performs better than the LMPC one in this experiment with velocities and accelerations larger than those in Sections VI-A1 and VI-A2.

Trajectories with higher velocities and accelerations might be successfully tracked by means of the proposed NMPC, but were not tested due to safety reasons since the HRPCable prototype was designed to handle heavy payloads at relatively limited velocities and accelerations. Nevertheless, the applicability of the proposed NMPC scheme to prototypes designed for high velocities and accelerations (such as [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF], [START_REF] Kawamura | Development of an ultrahigh speed robot FALCON using wire drive system[END_REF]) would require further investigations. This issue is mainly related to the computing time required to solve the OCP [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF]. On the one hand, a typical NMPC implementation would use numerical solvers implemented by experts (as discussed in [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF]). This would reduce the computation time in comparison to the controller used in HRPCable, which was implemented from scratch. On the other hand, higher velocities and accelerations may lead to faster changes in the dynamic model [START_REF] Alikhani | Sliding Mode Control of a Cable-driven Robot via Double-Integrator Sliding Surface[END_REF]. The solution of the optimization problem ( 16) may then demand a larger number of iterations of the numerical solver.

VII. CONCLUSIONS

This paper introduced an NMPC scheme for the position tracking of CDPRs. The OCP considered in the NMPC formulation handles explicitly the cable tensions and their limits. Accordingly, a substantial advantage of the proposed control scheme over typical state-of-the-art control strategies lies in its ability to operate on WFW boundary without failing. The proposed OCP formulation minimizes both the tracking error and the Wrench Equivalent Optimality (WEO). The latter is a non-negative measure able to evaluate whether the wrench generated by a given set of cable tensions can be generated by an alternative set with smaller 2-norm. The minimization of both the tracking error and the WEO is shown to be consistent, in the sense that a trajectory with minimal tracking error can be generated by cable tensions with minimal WEO. As a result, the stability of the closed-loop system can be analyzed and sufficient conditions for uniform asymptotic stability were deduced using tools for the analysis of NMPC schemes without terminal constraints and costs. Moreover, the proposed control scheme was validated experimentally on a 6-DoF CDPR driven by eight cables. In comparison to a previously introduced LMPC scheme, the tracking accuracy was substantially reduced.

While the experimental results showed that the proposed control scheme does not fail even in the presence of an unfeasible desired trajectory, the stability analysis presented in this paper is valid only for feasible trajectories. Future works may thus deal with the extension of the stability analysis of the proposed NMPC scheme to the case of unfeasible desired trajectories, aiming at obtaining the conditions under which the NMPC scheme remains stable while following a feasible trajectory as close as possible to the unfeasible desired one. Furthermore, the studied dynamic model considers that the inner control loop managing the cable tensions is able to generate the set of desired cable tensions instantly. Future works may consider, instead of [START_REF] Alikhani | Sliding Mode Control of a Cable-driven Robot via Double-Integrator Sliding Surface[END_REF], an augmented dynamical model taking into account the characteristics of this inner control loop. This approach may notably improve the vibration attenuation capabilities of the proposed NMPC scheme. Besides, a robustness analysis could evaluate the influence on the stability of the closed-loop system of numerical and modeling errors involved in [START_REF] Alikhani | Sliding Mode Control of a Cable-driven Robot via Double-Integrator Sliding Surface[END_REF]. Such an analysis may also consider the presence of measurement errors and disturbances.

APPENDIX A DETAILS ON THE DYNAMIC MODELS

Neglecting cable sagging, the net wrench f(x x x, τ τ τ c ) applied on the platform by all the cable forces is

f(x x x, τ τ τ c ) = W(x x x) τ τ τ c , (58)
where the wrench matrix representing this linear mapping is given by W(x x x) = -J T (x x x) with J(x x x) the Jacobian matrix (see [START_REF] Bruckmann | Wire Robots Part I: Kinematics, Analysis & Design[END_REF]). Additionally, it is well-known that the angular velocity ω ω ω is not equal to the time derivative of the orientation vector ψ ψ ψ. Hence, the following relations need to be introduced:

ṗ p p T ω ω ω T T = S(x x x) ẋ x x p p p T ω ω ω T T = S(x x x) ẍ x x + Ṡ(x x x, ẋ x x) ẋ x x, (59)
with S(x x x) a square matrix mapping ẋ

x x = ṗ p p T ψ ψ ψ T T into ẋ x x T ω ω ω T T
, and Ṡ(x x x, ẋ x x) is the component-wise time derivative of S(x x x) (see for instance [START_REF] Hughes | Spacecraft attitude dynamics[END_REF]Section 2.3] or [START_REF] Siciliano | Robotics: Modelling, Planning and Control[END_REF]Section 3.6]).

The mobile platform equations of motion can be obtained with the Newton-Euler formalism, leading to the following dynamic model when the cable distributed mass is neglected:

M (x x x) Ṡ(ẋ x x, x x x) ẋ x x + S(x x x) ẍ x x + C (ẋ x x, x x x) S(x x x)ẋ x x = g(x x x) + W(x x x) τ τ τ c , (60)
where matrices M and C are given by

M (x x x) = m p I -m p ĉ(x x x) m p ĉ(x x x) H(x x x) and (61)
C (x x x, ẋ x x) ẋ x x = m p ω ω ω ω ω ω c(x x x) ω ω ω H(x x x) ω ω ω . (62)
The scalar m p is the platform mass and I is an identity matrix with dimensions suitably chosen. For any pair of vectors v 1 , v 2 ∈ R 3 , the matrix v1 ∈ R 3×3 represents the skew symmetric matrix such that v1 v 2 = v 1 × v 2 , with × being the cross product. Denoting as c p the position vector of the platform center of mass written in the coordinate system O p attached to the platform (with O p shown in Figure 2), the corresponding vector in the global coordinate system O r is computed from c p using the rotation matrix R(ψ ψ ψ) as follows:

c(x x x) = R(ψ ψ ψ)c p = c x c y c z T . (63)
The angular velocity ω ω ω is obtained from x x x and ẋ x x in accordance with [START_REF] Merlet | A generic numerical continuation scheme for solving the direct kinematics of cable-driven parallel robot with deformable cables[END_REF]. The matrix H is defined as H(x x x) = R(ψ ψ ψ) I G R(ψ ψ ψ) T + m p ĉ(x x x) ĉ(x x x) T where I G is the platform inertia matrix with respect to the coordinate system O p . The vector of gravitational forces is g(x x x) = m p g 0 0 -1 -c y c x 0 T , with g the gravitational acceleration. The z-axis of the fixed reference frame is assumed to be vertical directed upward. Additionally, denoting M(x x x) = M (x x x) S(x x x) and

C(x x x, ẋ x x) = C (x x x, ẋ x x) S(x x x) + M (x x x) Ṡ(x x x, ẋ x x), (64)
the dynamic system (60) can be rewritten as

M(x x x) ẍ x x + C(x x x, ẋ x x) ẋ x x = g(x x x) + W(x x x) τ τ τ c . (65)
In order to define a discrete-time dynamic model based on (1), variables x k+1 and ẋk+1 should be computed based on known x k , ẋk and τ τ τ k . In accordance with (4), the cable tension vector τ τ τ k is considered as constant during a controller cycle period ∆t. Therefore, the solutions x x x : R → R n and ẋ

x x : R → R n of the Initial Value Problem (IVP) given by ẍ

x x = M(x x x) -1 g(x x x) + W(x x x) τ τ τ k -C(x x x, ẋ x x) ẋ x x , (66a) x x x(t k ) = x k and ẋ x x(t k ) = ẋk , (66b)
lead to a transition mapping φ φ φ

x : R n × R n × R m → R 2n such that x x x(t k + ∆t) ẋ x x(t k + ∆t) = φ φ φ x x k , ẋk , τ τ τ k . (67)
Using a more compact notation, the vector of the next state (pose and velocity)

y k+1 = x x x(t k + ∆t) T ẋ x x(t k + ∆t) T T (68)
obtained with initial conditions

x x x(t k ) T ẋ x x(t k ) T T = y k = x T k ẋT k T , (69)
and application of constant cable tensions τ τ τ k ∈ R m is described with the transition mapping φ φ φ y : R 2n × R m → R 2n according to y k+1 = φ φ φ y (y k , τ τ τ k ).

In general, the exact closed-form solution of the nonlinear IVP (66) cannot be determined. This IVP is solved with some numerical method (e.g. Euler or Runge-Kutta methods). It is also known that a numerical solution of (66) can be computed with an arbitrarily small error (according to the convergence theory of ordinary differential equations [START_REF] Deuflhard | Scientific Computing with Ordinary Differential Equations[END_REF]Chapter 4]). In the present work, the numerical solution of the IVP (66) using the Euler method is considered. Accordingly, the discretetime form of the state-space representation of the dynamic model ( 1) is approximated by [START_REF] Williams | Planar translational cabledirect-driven robots[END_REF]. This rather simple solution is discussed in the remainder of this appendix.

Consider an instant t k = k ∆t, pose x k ∈ R n , velocity ẋk ∈ R n and cable tensions τ τ τ k ∈ R m . The platform acceleration at the instant t k is given by ẍ

x x(t k ) = M(x k ) -1 W(x k ) τ τ τ k + g(x k ) -C(x k , ẋk ) ẋk = ẍk (71)
The Euler integration method [71, Chapter 4] assumes the following approximation:

x x x(t k + k ∆t) = x x x(t k ) + ẋ x x(t k ) ∆t ẋ x x(t k + k ∆t) = ẋ x x(t k ) + ẍ x x(t k ) ∆t (72)
or, using the discrete counterparts of x x x, ẋ x x and ẍ x x,

x k+1 = x k + ẋk ∆t ẋk+1 = ẋk + ẍk ∆t. (73)
The application of ( 71) in (73) leads to [START_REF] Williams | Planar translational cabledirect-driven robots[END_REF].

APPENDIX B ON THE PROOF OF THEOREM 1

The results presented in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 6] are based on upper bounds of the cost functional. More precisely, these upper bounds are studied with respect to the optimal value function, which is defined in the following. For a K ∈ N, the optimal value function V K : N × Y → R is defined for an instant k and state y as follows:

V K (k, y) = inf τ τ τ(•)∈U K k J K (k, y, τ τ τ(•)) , s. t. y τ τ τ(•) ( j, y) ∈ Y ∀ j ∈ N 1,K (74)
with the cost functional J K defined in [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF]. The scalar V K (k, y) is also known as the cost-to-go for given instant k and state y in the context of dynamic programming [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. Note that, if the following optimal control problem possesses a global minimum

τ τ τ * (•) = arg min τ τ τ(•)∈U K k K ∑ j=0 k + j, y τ τ τ(•) ( j, y) , τ τ τ k+ j , s. t. y τ τ τ(•) j, y ∈ Y ∀ j ∈ N 1,K (75)
then τ τ τ * (•) satisfies

J K (k, y, τ τ τ * (•)) = V K (k, y). (76)
In order to deduce Theorem 1, let us first consider a simplified case with a constant desired state

y d (k) = y d ∀ k ∈ N (77)
with a constant vector y d ∈ Y. In this case, the closed-loop system is time invariant, the proposed NMPC scheme obtained with Algorithm 1 may be written as in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Algorithm 3.1] and the [43, Theorem 6.24] can be directly applied as it is. Moreover, the value function V K is constant with respect to time and may be simplified to

V K (y) = V K (k, y) ∀ k ∈ N. (78)
Similarly, the stage cost and minimal stage cost can be simplified according to

c (y, τ τ τ) = (k, y, τ τ τ) * c (y) = * (k, y) ∀ k ∈ N. (79)
The conditions for uniform asymptotic stability stated in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Theorem 6.24] addressing time-invariant systems are summarized as follows.

(i) There exist α 1 , α 2 ∈ K ∞ such that for every y ∈ Y the following inequalities hold:

α 1 ( y -y d ) * c (y) α 2 ( y -y d ); (80a)
(ii) ([43, Assumption 6.3] with linear B K ) For each y ∈ Y and K ∈ N, there exists γ K ∈ R with γ K < ∞ so that the optimal cost functional V K (y) satisfies

V K (y) γ K * c (y). (80b)
Additionally, [START_REF] Grüne | Nonlinear model predictive control[END_REF]Lemma 6.6] shows that, if, for each y ∈ Y and N ∈ N, there exist τ τ τ(•) ∈ U N 0 and real constants C 1, σ ∈ (0, 1) such that

c y τ τ τ(•) (n, y), τ τ τ n C σ n * c (y) ∀ n ∈ N 1,N-1 , (81)
then condition (80b) is satisfied with

γ K = C 1 -σ K 1 -σ (82) and sup K γ K < ∞ ∀ K ∈ N.
Condition (81) is referred to as exponential controllability. Under this condition, there exists τ τ τ(•) that generates a trajectory with stage cost exponentially converging to zero in time.

Therefore, the conditions (80) are equivalently given as (80a) and (81). Considering (77) (and (79), consequently), the conditions of Theorem 1 are simplified to: (a) There exist α 1 , α 2 ∈ K L ∞ such that where e e e 0 = e e e(0), ė e e 0 = ė e e(0). The expressions of the error e e e and its time derivative satisfy e e e(t) = e -t/2 a 1 sin(ω t) + a 2 cos(ω t) e -1/2 t ( a 1 + a 2 ) (90)

and ė e e(t) = e -t/2 a 3 (sin(ω t) + a 4 cos(ω t))

e -1/2 t ( a 3 + a 4 ),

with a 3 = a 1 /2ω a 2 and a 4 = a 2 /2 + ω a 1 .

A tedious but straightforward computation using (89) leads to the bound e e e(t) 2 + ė e e(t) 2 c ey e -t ( e e e 0 2 + ė e e 0 2 ), with c ey = (16 + 2 √ 3)/3. Consider k x , k ẋ ∈ R n + , subvectors of k y (as in Section VI), and matrices K x , K ẋ ∈ R n×n diagonal matrices with K x = diag(k x ), K ẋ = diag(k ẋ). Denoting also

K max y = max i (k y,i ) and K min y = min i (k y,i ) , (92)
the following inequalities hold:

e e e(t) which is the continuous counterpart of (29a) with C 1 = c ey K max y /K min y < ∞ and σ 1 = e -1 < 1.

APPENDIX D NUMERICAL SOLUTION OF THE OCP [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF] As discussed in Section VI, the presented experimental results were obtained with a control scheme implemented from scratch. Whereas the implementation of the basic mathematical operations are not noteworthy in this paper, the optimization techniques used in the solution of the OCP (16) deserves special attention. Since [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF] is solved in real time, the implemented methods should combine good numerical accuracy and short computation time. The present section briefly discusses the optimization techniques used.

Besides the optimization problem corresponding to (16) itself, the solution of this OCP also includes the computation of the WEO, in accordance with [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF]. More precisely, the OCP obtained considering explicitly the upper [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF] and lower [START_REF] Chellal | Model identification and vision-based H-infinity position control of 6-DoF cable-driven parallel robots[END_REF] level optimization problems can be expressed as the bilevel minimization given by min τ τ τ(•)∈U W(x k+ j ) τ τ τ k+ j = W(x k+ j ) τ τ τ k+ j ,

τ τ τ ∈ U ∀ j ∈ N 0,h p -1
The present appendix discusses the numerical solution of both lower and upper level problems. Accordingly, Section D-A addresses the solution of ( 15) and Section D-B briefly describes the optimization strategy used in the solution of [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF].

A. Computation of the WEO

The optimization problem ( 15) is a convex inequality constrained QP problem with non-empty feasible region (as discussed in the proof of Lemma 1). The results presented in Section VI were obtained with a control scheme that solves this optimization problem with an active-set method for convex QPs [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5]. More in detail, this method can determine the exact solution of the problem min τ τ τ τ τ τ T τ τ τ (93a) s. t. A eq τ τ τ = b eq (93b)

A ineq τ τ τ b ineq (93c)
within a finite number of iterations. Matrices A ineq , A eq and vectors b ineq , b eq can be chosen so that constraints (93b)-(93c) are equivalent to (15b)-(15c), respectively. One may notably note that A eq = W(x). Each iteration consists of the solution of an equality constrained QP in which, besides (93b), some of the inequalities (93c) are taken as equalities and the remaining ones are neglected. More precisely, for an iteration j ∈ N, the equality constrained QP

with A ineq, j and b ineq, j composed of some of the rows of A ineq and b ineq , respectively. The solution of (94) can be obtained finding the vectors of cable tensions τ τ τ and of Lagrange multipliers µ µ µ that satisfy the first-order optimality conditions (KKT conditions)

τ τ τ + A T j µ µ µ = 0 A j τ τ τ = b j . ( 96
)
The system of equations (96) can be solved straightforwardly applying the QR factorization [74, Section 5.2] of the matrix A T j . This method can be particularly efficient since the QR factorization of A T j+1 obtained in the iteration j + 1 may be based on the factorization of A T j performed in the previous iteration. As detailed in [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5], matrices A j+1 and A T j differ at most of one row, and this similarity can be used in order to update the QR factorizations of successive iterations. Furthermore, one may note that the QR factorization of A T j is closely related to the computation of the null space of W(x), which is typically used in the existing state-of-the-art cable tension distribution algorithms.

with the case addressed in Section VI, the state constraint set is considered as Y = R 2n . Therefore, the QP inequality constraints consist uniquely of the cable tension limits, i.e.: τ τ τ min τ τ τ k+ j τ τ τ max , ∀ j ∈ N 0,h p -1 .

(101)

One may note that the constraints (101) lead to a box constrained (or bound-constrained) QP. The gradient projection method [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.7] is well known as an appropriate method for such QP problems. As illustrated in Figure 17, the proposed strategy couples the gradient projection with the conjugate-gradient method [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 5] in order to obtain faster convergence. The implemented SQP algorithm can be classified as a nonlinear gradient projection method, briefly described in [START_REF] Nocedal | Numerical optimization[END_REF]Section 18.6].

Fig. 1 .

 1 Fig. 1. HRPCable prototype, a 6-DoF fully-constrained CDPR.

Fig. 2 .

 2 Fig. 2. CAD view of the HRPCable prototype.

Fig. 3 .

 3 Fig. 3. Block diagram of a position tracking control scheme integrating the proposed NMPC.

Definition 1 .

 1 Consider the sequences x d (•), ẋd (•) ∈ S n of desired poses and velocities, respectively, and y d (•) ∈ S 2n , with y d (k) = x d (k) T ẋd (k) T T , the desired trajectory. The sequence y d (•) is called feasible if there exists a sequence of nominal desired cable tensions τ τ τ nd (•) ∈ S m such that, for y 0 = y d (0),

Definition 2 .

 2 Consider the NMPC Algorithm 1 with prediction horizon h p ∈ N and a feasible sequence of desired states y d (•) ∈ S 2n . The resulting closed-loop system

Fig. 4 .

 4 Fig. 4. An illustrative example of a function β ex ∈ K L given by β ex (r,t) = r 3 /t.

Assumption 1 .

 1 Consider the error e y : N × Y → R given by e y (k, y) = yy d (k) 2 K y .

Assumption 2 .

 2 Consider y, y + ∈ Y and τ τ τ a , τ τ τ b ∈ U. If these vectors satisfy y + = φ φ φ y (y, τ τ τ a ) and

Fig. 5 .

 5 Fig. 5. Visited poses for the pick-and-place trajectories performed in the numerical simulations.

Fig. 6 .

 6 Fig. 6. Numerical results comparing LMPC and NMPC schemes with optimal gains.

Fig. 7 .

 7 Fig. 7. Monte Carlo simulation comparing LMPC and NMPC schemes for different controller gains.

K

  τ τ τ = 6 × 10 -7 I and h p = 4. (57d)

Fig. 8 .

 8 Fig. 8. Experimental results for τ τ τ max = 400 N. In the charts presenting the cable tensions, the dashed and solid lines represent the desired and measured tensions, respectively.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Visited poses for the pick-and-place trajectories performed in the experimental tests.

  Fig. 11.Experiments with τ τ τ max = 250 N. In the charts presenting the platform position and orientation, the dashed and solid lines represent the desired and actual trajectories, respectively.

  Fig. 12.Experimental results obtained with the control scheme proposed in[START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] with the reduced maximum cable tension τ τ τ max = 250 N. The robot operation is terminated at t ≈ 10.5 s, when the cable tensions become unfeasible.

Fig. 13 .

 13 Fig. 13. Illustration of the robustness test.

Fig. 15 .

 15 Fig. 15. Tracking errors of the LMPC and NMPC schemes in the experiment testing the robustness against payload uncertainties.

Fig. 16 .

 16 Fig. 16. High-speed trajectory results.

α 1 (

 1 yy d ) * c (y) α 2 ( yy d ) ∀ y ∈ Y. (83) (b) For all y ∈ Y, there exist feasible τ τ τ e (•) ∈ S m , real C < ∞ and σ ∈ (0, 1) satisfyingc y τ τ τ e (•) ( j, y), τ τ τ e j C σ j * c (y) , ∀ j ∈ N.(84)These conditions are equivalent to (80a) and (81). The time-varying case addressed in Theorem 1 is obtained considering [43, Assumptions 6.29 and 6.30], indicating that the results above can be straightforwardly extended to timevarying systems.APPENDIX C EXAMPLE OF TRIVIAL CONTROL LAW SATISFYINGASSUMPTION 1 Consider, for instance, a continuous-time control scheme that applies the wrenchf = M(x x x) ẍ x x d + K P (x x x dx x x) + K D (ẋ x x dx x x) + + C(x x x, ẋ x x) ẋ x x -g(x x x)(85)with diagonal gain matrices K P and K D . According to (1), the nominal closed-loop system dynamics is given by ẍx x = ẍ x x d + K P (x x x dx x x) + K D (ẋ x x dx x x)(86)or, written in terms of the error e e e(t) = x x x d (t)x x x(t), ë e e = -K D ė e e -K P e e e.(87)Choosing, for instance, K P = K D = I, the solution to the differential equation (87) can be written as e e e(t) = e -t/2 (a 1 sin (ω t) + a 2 cos (ω t)) ,(88)with ω = √ 3/2 and constant vectors a 1 , a 2 ∈ R n given by a 1 = e e e 0

y

  k+ j -y d (k + j) 2 K y + τ τ τ k+ jτ τ τ wo,k+ j 2 K τ τ τ s. t. y k+ j = x T k+ j ẋT k+ j T = y τ τ τ(•) j, y k ∈ Y ∀ j ∈ N 0,h pτ τ τ wo,k+ j = arg min τ τ τ τ τ τ 2 :

  T τ τ τ s. t. A j τ τ τ = b j (94)is considered, where matrix A j and vector b j are given byA j =A eq A ineq, j and b j = b eq b ineq, j

TABLE I SUMMARY

 I

OF COMPARATIVE ERRORS BETWEEN LMPC AND NMPC: ROOT MEAN SQUARE (RMS), STANDARD DEVIATION (SD) AND MAXIMUM ERRORS. RMS (SD) Maximum error TE OE TE OE LMPC 4.28 (2.17) mm 0.15 (0.07) • 7.47 mm 0.27 • NMPC 1.52 (0.76) mm 0.08 (0.04) • 2.34 mm 0.14 • Improvement 64.5 (64.8) % 44.6 (42.1) % 68.6% 48.1%

TABLE II SUMMARY

 II OF COMPARATIVE ERRORS BETWEEN LMPC AND NMPCREGARDING THE ROBUSTNESS AGAINST PAYLOAD UNCERTAINTIES.

		RMS (SD)		Maximum error
		TE	OE	TE	OE
	LMPC	4.66 (3.46) mm 0.52 (0.40) • 12.37 mm 1.05 •
	NMPC	3.92 (2.82) mm 0.31 (0.23) • 7.83 mm 0.63 •
	Improvement 16.0 (18.3) % 40.7 (42.4) %	36.7%	40.0%

HD video available in https://youtu.be/NIb_XYSVv84

This work was supported by the European Union's H2020 Program (H2020/2014-2020) under the grant agreement No. 732513 (HEPHAESTUS project). The authors also (greatly) acknowledge support of the European Union through FEDER grant n • 49793. Cables Platform Winches

B. Overall solution of (16)

The implemented control scheme solves the OCP ( 16) applying recursive elimination (in accordance with [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 12.1]). Therefore, the cost functional is computed with states y τ τ τ(•) (•, y k ) obtained by applying the transition mapping (8) recursively, as in [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF]. Alternative methods such as the full discretization and multiple shooting may also be used at the cost of increased dimension of the optimization problem. Details on these different approaches are discussed in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 12.1] and [START_REF] Rawlings | Model predictive control: Theory, Computation, and Design[END_REF]Section 8.5].

The application of the recursive elimination leads to a Nonlinear Programming Problem (NPP) corresponding to the OCP [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF]. Sequential Quadratic Programming (SQP) and interiorpoint methods are often used to solve NPPs. The applicability of such algorithms in the implementation of NMPC controllers is studied in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 12.3] and [START_REF] Rawlings | Model predictive control: Theory, Computation, and Design[END_REF]Section 8.7].

The use of an SQP algorithm is considered appropriate for the solution of the obtained NPP since the proposed OCP possesses a quadratic cost functional. Nocedal and Wright discuss the details of SQP methods in [73, Chapter 18]. These are iterative algorithms in which, at each iteration, the original nonlinear problem is approximated as a QP problem. A standard QP strategy may thus be used in the solution of each of these sub-problems.

In order to avoid the expensive computation of the second derivatives of the nonlinear cost function, a quasi-Newton BFGS approximation [73, Section 6.1] was used in the present work. This method estimates the Hessian of the cost function based on its first derivatives. While BFGS could be applied directly to the NPP obtained from [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF], it has been observed that better performance is obtained with the optimization strategy summarized in Figure 17. Let τ τ τ (•) ∈ U h p k be an initial guess for the optimal sequence of cable tensions at an instant k ∈ N and state y k ∈ Y. The sequence τ τ τ d (•) is defined as the WE-optimal tensions of τ τ τ (•) given by

with y (•) the sequence of predicted states within the prediction horizon, obtained according to

The WE-optimal tensions τ τ τ d (•) in ( 97), which are computed with the method presented in Section D-A, are taken as reference to the cable tension optimization. Accordingly, one may define a simplified stage cost τ τ τ d (•) : N×Y×U → R given by τ

The NPP corresponding to the stage cost τ τ τ d (•) can be written as min

A standard SQP algorithm is used in order to solve this NPP, resulting in an updated τ τ τ (•) solution of (99) (as depicted in Figure 17). This procedure is applied iteratively for NPPs with updated τ τ τ d (•), obtained according to (97). The overall algorithm is terminated once the difference between sequences τ τ τ d (•) and τ τ τ (•) is less than a given threshold. More precisely, the optimal sequence of cable tensions τ τ τ * d (•) is taken as τ τ τ (•) when the following convergence condition is satisfied:

with a real positive ε, representing the desired convergence threshold. The optimal τ τ τ * d (•) is used as initial guess τ τ τ (•) in the next sampling time k + 1.

The implemented SQP algorithm applies the quasi-Newton BFGS estimation of the Hessian of (99). This estimation results in an inequality constrained QP problem. In accordance