Obtenir De

Docteur Le Grade De

Université De L'

Montpellier De

De L'université De Montpellier

Présentée Par

João Cavalcanti

Santos Le

Stéphane Caro

Jacques Gangloff

Ionela Prodan

Ben- Jamin Navarro

Robin Passama

Arnaud Meline

Jean Baptiste Izard

Jose Gorrotxategi

Ibon Christoph Martin

Eskudero Gurrutxaga

Kepa Iturralde

Marc Fabritius

Cécile Lukasik

Gwladys Imbert

Anne-Laure Irastorza

Najima Maatoug

Vinayak

Karine Thibault

Maxime Virginie Lucas

Maxence Benjamin Guillaume

THÈSE POUR OBTENIR LE GRADE

saturation (A.8) may reduce the necessary number of iterations. . . . . . . . . .

List of Figures

1.8

The NIST ROBOCRANE, controlled with a Joint Space control [START_REF] Albus | The NIST robocrane[END_REF]. . . . . 1.9 A Joint Space Control scheme proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]. The main variables are:

The desired and actual pose x d and x, desired and measured motor positions q d and q, Jacobian matrix J, desired cable tensions τ d , dynamic motor friction constant F v , motor inertias I m , motor currents i, a diagonal constant matrix C T and PD gains K p and K v . The blocks TD and IK stand for the Tension Distribution algorithm and Inverse Kinematics, respectively. Time derivatives are denoted using the dot notation. . . . .

1.10

The Operational Space Control proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]. The variables that were not introduced in Figure 1.9 are: The pseudo-inverse jacobian matrix J + , static motor friction constant F s , tracking error of the pose ∆x and PID gains on operational space K p , K i and K v . . . . . . . . . . . . . . . . . . .

General Introduction

Cable-Driven Parallel Robots (CDPRs) represent a subclass of parallel kinematic manipulators in which rigid links are replaced by cables. They consist mainly of a mobile platform driven by cables, which transmit the forces generated by winches. Accordingly, the platform motion is managed controlling the winch motors.

Since a large length of cable may be wound on the drum of each winch, CDPRs can have a large workspace. Furthermore, thanks to the reduced inertia of the cables, CDPRs can attain high velocities and accelerations. In addition, cables are able to transmit large forces, which make possible the handling of heavy payloads. Among other advantages, these characteristics motivated studies applying CDPRs for high speed manipulation, simulations in virtual reality, handling of heavy payloads in large workspaces, etc. As a matter of fact, the design, control and analysis of CDPRs are well established research subjects [START_REF] Pott | Cable-Driven Parallel Robots: Theory and Application[END_REF].

The work presented in this thesis has been made in the framework of the european H2020 project called Hephaestus [START_REF]Hephaestus Project[END_REF]. The main goals of this project include the design and control of a CDPR able to install the Curtain Wall Modules (CWMs) of a building facade. Inspired by the specificities of the Hephaestus project, this thesis contributes to the development of advanced control techniques.

Clearly, the mounting of CWMs demands high positioning precision (10 mm). In addition, safety and disturbance rejection capabilities are major concerns considering that the robot should work in a construction site. Moreover, the control scheme should be compatible with a specific set of industrial hardware and software. Finally, due to design constraints, the robot needs to operate close to the boundaries delimiting the safe conditions of the system operation.

Considering this context, it was noted that the state-of-the-art control schemes would lead to a crucial issue. Namely, the operation close to system constraints with a significant influence of disturbances may lead to undesired responses. This is due to the fact that the existing strategies do not consider system constraints as an integral part of the main control scheme. In order to overcome this issue, this thesis introduces Model Predictive Control (MPC) schemes (linear and nonlinear) able to explicitly handle the system constraints within the formulation of the controller. Numerical and experimental results validate the proposed strategies. Moreover, the stability of a nonlinear MPC (NMPC) scheme is also analyzed.

Before applying the proposed MPC schemes, several issues related to the experimental implementation had to be solved. Some of them are addressed in this manuscript. In particular, a chapter is dedicated to contributions related to the cable tension control and forward kinematics of CDPRs. In addition, the programming solutions proposed to the development of the control schemes in an industrial software are briefly examined.

Organization of the Thesis

A brief description of the chapters and appendices of this thesis are presented in the following.

• Chapter 1 presents the context of this thesis and review the pertinent state-ofthe-art on the control of CDPRs. Considering this context, an issue raised by the existing control schemes is identified and discussed. It is noted that a control scheme suitable for the project Hephaestus should consider the system constraints as an integral part of the main controller. Some details on MPC are reviewed, demonstrating that this control architecture may be suitable in the context of this thesis;

• Chapter 2 introduces the implemented kinematic model and the control of cable tensions. A forward kinematic algorithm considering the pulley geometry is proposed. These results were published in (I). Regarding the cable tension control, two different approaches are considered. Numerical simulations and experimental tests indicate that a velocity-based scheme is suitable for the studied systems;

• Chapter 3 introduces a linear MPC for the position tracking control of CDPRs.

The strategy is validated numerically and experimentally, indicating that the proposed MPC controller overcomes the issues described in Chapter 1. Simulation results were published in (III) and experimental ones in (II);

• Chapter 4 proposes a nonlinear MPC scheme. In contrast to the linear MPC, this approach leads to satisfying results even with an increased incidence of nonlinearities. The stability of the nominal closed-loop system is analyzed. The comparison between the linear and nonlinear MPC is performed through realtime experimental tests and numerical simulations;

• Appendix A briefly discusses the development of a C++ numerical library compatible with the industrial software TwinCAT. Thanks to this library, complex matrix operations can be programmed in an intuitive manner, facilitating the implementation of the proposed MPC schemes;

• Appendix B presents a closed-form tension distribution algorithm able to minimize the infinity-norm of the cable tensions vector. The optimality of the algorithm output is proved and experimental results are discussed. The proposed algorithm is based on the results published in (IV).

• Appendix C addresses the stability criterion that is used in the analysis of the proposed NMPC scheme. The results introduced in [START_REF] Grüne | Nonlinear model predictive control[END_REF] that are used Chapter 4 are discussed more in detail.

List of Publications -International Conference Proceedings

Introduction and State of the Art

Cable-Driven Parallel Robots (CDPRs) present several particular advantages when compared to rigid link robots. Section 1.1 discusses the reasons that lead to these advantages and exemplify some of the various applications of CDPRs that take advantage of them. Section 1.2 then introduces the main goals of the project Hephaestus that form the context of this thesis. The main models of CDPRs proposed in the literature are briefly presented in Section 1.3. Notably, this section examines the various assumptions that may be taken and the main consequences on the resulting models. A summary of the state of the art on position tracking of CDPRs is introduced in Section 1.4. The vast majority of the studies presented in this section are coupled to some tension distribution algorithm. This subject is addressed in Section 1.5. Section 1.6 presents existing cable tension control schemes. The very beginning of Section 1.7 highlights an important limitation of the coupling between schemes presented in Section 1.4 and algorithms of Section 1.5. Model Predictive Control (MPC) is then proposed as a solution to overcome this drawback. Section 1.7 presents details on this type of feedback control strategies. Finally, Section 1.8 summarizes the main contributions presented in this thesis.

Cable-Driven Parallel Robots

The main components of a Cable-Driven Parallel Robot (CDPR) are cables, winches and a mobile platform. Each cable has one end attached to the platform and the other one wound on a drum. Each one of these drums is coupled to a motor, composing the winches. Figure 1.1 depicts a CDPR and highlights these main elements. Altogether, the cables transmit the torques applied on the drums so that the platform motion may be driven by the actuation of the motors.

CDPR characteristics and main applications

CDPRs have several advantages compared to rigid link robots. Since the rigid links are replaced by cables that may present a large length wound on the drums, the mobile platform may be displaced over a large workspace. A few examples among the numerous robots that take advantage of this feature are the well known Skycam [START_REF] Lawrence | Skycam: an Aerial Robotic Camera System[END_REF] and large-scale radio telescopes [START_REF] Yao | Dimensional optimization design of the four-cable-driven parallel manipulator in FAST[END_REF][START_REF] Lambert | Implementation of an aerostat positioning system with cable control[END_REF]. Figure 1.2 (a) depicts the large dimension telescope FAST installed in China. Moreover, the use of cables may be very efficient to handle heavy payloads, rendering the application of CDPRs favorable for the handling of heavy objects [3,[START_REF] Albus | The NIST robocrane[END_REF], as the CoGiRo prototype shown in Figure 1.

(b).

The combination of a large workspace with the reduced visual interference obtained using cables instead of rigid links makes cable-driven parallel mechanisms a pertinent solution for haptic devices [4,[START_REF] Sato | Development of string-based force display: SPIDAR[END_REF][START_REF] Martin | Cable tension control and analysis of reel transparency for 6-dof haptic foot platform on a cable-driven locomotion interface[END_REF][START_REF] Martin | Hybrid control with multicontact interactions for 6 DOF haptic foot platform on a cable-driven locomotion interface[END_REF]. Rehabilitation is a typical application of such solutions. The haptic device Inca, commercialized by Haption, is shown in Figure 1.3.

As a subclass of parallel robots, CDPRs typically present their actuators fixed to the inertial base. In addition, the inertia of the moving links may be drastically reduced using cables instead of rigid links. As a result, several works focused on the capability of CDPRs to reach high velocities and accelerations [START_REF] Kawamura | Development of an ultrahigh speed robot FALCON using wire drive system[END_REF][START_REF] Dekker | Design and testing of an ultra-high-speed cable robot[END_REF][START_REF] Zhang | Optimization and implementation of a high-speed 3-DOFs translational cable-driven parallel robot[END_REF][START_REF] Hiller | Design, analysis and realization of tendon-based parallel manipulators[END_REF]. For instance, the FALCON prototype is able to generate 13 m/s of maximal velocity and 43 g of maximal acceleration [START_REF] Kawamura | Development of an ultrahigh speed robot FALCON using wire drive system[END_REF]. Similarly, motion simulators may take advantage of the high dynamic capabilities obtained with the architecture of CDPRs [START_REF] Miermeister | The CableRobot simulator large scale motion platform based on Cable Robot technology[END_REF] (see Figure 1.4).

In contrast to these advantages, CDPRs also have important drawbacks. The cable elasticity may introduce undesired vibrations [START_REF] Baklouti | Dynamic and oscillatory motions of cable-driven parallel robots based on a nonlinear cable tension model[END_REF][START_REF] Baklouti | Vibration reduction of Cable-Driven Parallel Robots through elasto-dynamic model-based control[END_REF][START_REF] Weber | Active vibration canceling of a cable-driven parallel robot in modal space[END_REF][START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF]. Similarly, substantial degradation of the positioning precision may be obtained when neglecting cable elongation [START_REF] Schmidt | Increase of position accuracy for cabledriven parallel robots using a model for elongation of plastic fiber ropes[END_REF]. Moreover, the cable tensions should be kept sufficiently high in order to avoid the presence of slack cables, which not only may lead to the malfunctioning of the winches but also degrades positioning accuracy.

CDPR types

Whereas rigid links may transmit both compression and tensile forces, cables can only operate in tension. Considering this limitation, the Wrench Closure Workspace (WCW) defined in [START_REF] Gouttefarde | Analysis of the wrench-closure workspace of planar parallel cable-driven mechanisms[END_REF] is "the set of poses of the platform for which any wrench can be gener- ated at the platform by tightening the cables". Here, the term pose refers to the combination of position and orientation. A mobile platform pose within the WCW is said to be fully-constrained.

It is interesting to note that many CDPRs present nonexistent WCW. This is the case of the suspended CDPRs, such as the ones depicted in Figures 1.2 (a) and (b). Since every cable force applied on the platform is directed upwards, this type of CDPR can not generate a force directed downwards. Therefore, suspended CDPRs rely on the gravitational forces in order to constraint the platform.

In contrast, if the WCW represents a substantial part of the robot workspace, the CDPR is classified as fully-constrained. While suspended CDPRs may have superior payload capacity, fully constrained CDPRs typically present better precision, stiffness and dynamic capabilities. Robots depicted in Figures 1.3, 1.4, and 1.5 are examples of fully-constrained CDPRs. It is worth noting that a fully-constrained CDPR necessarily presents a number of cables greater than the number of Degrees of Freedom (DoF) of the mobile platform. This actuation redundancy yields some particularities on the control. These particularities play an important role on the main results presented in this thesis. Sections 1.5 and 1.7 discuss these particularities more in detail. At any rate, even though the studies presented in this thesis are focused mainly on redundant CDPRs, they may be directly applied for non-redundant and suspended CDPRs.

Conversely, the platform pose cannot be determined geometrically for CDPRs with less cables than DoF. In addition to the kinematic model, a static or dynamic model should be considered in order to determine the platform pose. Accordingly, such robots are classified as under-constrained CDPRs. This class of CDPRs was studied, for instance, in [START_REF] Carricato | Stability analysis of underconstrained cable-driven parallel robots[END_REF][START_REF] Abbasnejad | Direct Geometrico-static Problem of Underconstrained Cable-Driven Parallel Robots With n Cables[END_REF]. In spite of the interesting research subjects devoted to underconstrained CDPRs, this thesis does not deal with this kind of robot.

Context of the Thesis

The contributions of this thesis were developed within the framework of the European Horizon 2020 project Hephaestus [START_REF]Hephaestus Project[END_REF]. The main goal of this project is to develop a robotic solution for the installation of Curtain Wall Modules (CWM) on building facades. Figure 1.5 presents the Hephaestus CDPR prototype and illustrates its main application. Indeed, one may note that the choice of using a CDPR for this application is in accordance with the advantages introduced in the previous section, notably the possibility to reach large workspaces (the size of a building facade) and the capability to handle heavy payloads (approximately 1300 kg in this specific case).

The positioning of the CWMs is performed with the aid of suction cups installed on the mobile platform. Furthermore, the prototype also assists the installation of the brackets on the building slabs. This task includes the drilling of the slabs, positioning and mounting of the brackets. In spite of the interesting challenges involved in the bracket installation, the contributions presented in the present manuscript are not dependent on technical details related to them. For this reason, among the extensive technical data involved in the whole project, just a few relevant inputs are detailed in Table 1.1.

The main contributions of LIRMM in the Hephaestus project are (i) the geometric design optimization of the CDPR and (ii) the development of advanced control techniques. Since most of the remaining tasks of the project depend on the outputs of (i), the results related to this task were delivered in the first half of the project. These results are presented in [START_REF] Hussein | Geometric Optimization of a Large Scale CDPR Operating on a Building Facade[END_REF] and my contributions in this task were minor. Therefore, this thesis is mainly focused on the development of CDPR control strategies.

The time length in which the Hephaestus prototype was fully built and available for tests is short in comparison to the duration of the project. Moreover, this available time needs to be shared between the partners of the consortium who participate to the experimental tests. As a result, knowing that the experimental activities were crucial and might reveal non-predicted difficulties, a secondary experimental set-up was used. This set-up is a CDPR prototype, HRPCable, with characteristics rather similar to the ones of the Hephaestus prototype. The HRPCable is shown in Figure 1.1 and Table 1.1 summarizes its main features. Indeed, most of the experimental results presented within this thesis were performed with this prototype. Before addressing the modeling and control of CDPRs, the main concerns related to them should be summarized considering the specificities of the Hephaestus project:

• A crucial challenge is the mounting of the CWMs on the brackets. High precision is necessary and a versatile control strategy may be needed, since the assembly is performed in four brackets simultaneously (see the illustration in Figure 1.6). The position of the brackets are not necessarily known with the precision required to the assembly; • Since the robot should work outdoors, in a non-controlled environment, the proposed control scheme should possess external disturbance rejection capabilities; • Given that the robot is intended to operate on a construction sites (with the possible presence of non-trained personnel) and the payload of the machine is large, safety is a critical concern; • Due to the reduced available area in front of the building facade, the nominal designed workspace presents a small margin dedicated to uncertainties. In other words, some poses necessary to the accomplishment of the desired tasks are close to the nominal workspace limits. Therefore, the proposed control scheme should be able to safely handle the robot operation close to system constraints; • Considering that the robot is intended to operate on a construction sites (with the possible presence of non-trained personnel) and the payload of the machine is heavy, safety is a critical concern; • Knowing that the machine is intended to be commercialized, the used hardware and software should be compatible with the ones commonly used in the industry. Since the commonly available industrial softwares are not compatible with advanced mathematical programming libraries, this last aspect leads to remarkable constraints related to the implementation of advanced control schemes.

Modeling

Most of the control strategies discussed in Section 1.4 are model-based control schemes. Indeed, the design of these control strategies takes advantage of the wellknown modeling of CDPRs. For this reason, it is interesting to present some details on the kinematic and dynamic models of CDPRs before addressing the control design itself.

The kinematic model of a CDPR correlates the pose and velocity of the mobile platform to the motor positions and velocities. The model relating a set of joint positions to a known pose of the platform is the Inverse Kinematic (IK) model. Conversely, the calculation of a platform pose for a known set of motor positions is the Forward Kinematic (FK) model. The different mathematical models used to implement et IK and FK may be classified with respect to their assumptions, leading to different levels of precision.

One may neglect the cable elasticity and mass, considering that the profile of each cable is a straight line linking an attachment point to the corresponding drawing point. Even though these may seem rather naive assumptions, several important studies were based on the corresponding model [START_REF] Ryeok | 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][START_REF] Pusey | Design and workspace analysis of a 6-6 cable-suspended parallel robot[END_REF]. On the other hand, several approaches were proposed seeking more precise models, as discussed in the following.

Typically, previous works search for a compromise between model simplicity and precision in the context of each specific study. Firstly, one may address the problem that The influence of these changes, which is often referred to as pulley kinematics, may be significant [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF] and may be taken into account in the control of CDPRs [START_REF] Schmidt | Implementing Extended Kinematics of a Cable-Driven Parallel Robot in Real-Time[END_REF][START_REF] Cavalcanti | Redundancy Resolution integrated Model Predictive Control of CDPRs : Concept , Implementation and Experiments[END_REF].

Moreover, the fact that in reality the cables are elastic and do not present a straight line profile may be also studied. In cases where the cable tensions are low and the lengths are large, the influence of the cable distributed mass may be non negligible. Therefore, as an example, it is not surprising that this phenomenon was considered in the design of the CDPR FAST illustrated in Figure 1.2 (a). Yao et al. obtained in [START_REF] Yao | Dimensional optimization design of the four-cable-driven parallel manipulator in FAST[END_REF] a good compromise between precision and computing time approximating the cable profile with a parabola. A model proposed in [START_REF] Max | Cable structures[END_REF] (usually referred to as Irvine model) is often used to obtain more precise results [START_REF] Kozak | Static analysis of cable-driven manipulators with non-negligible cable mass[END_REF][START_REF] Gouttefarde | Kinetostatic analysis of cable-driven parallel robots with consideration of sagging and pulleys[END_REF][START_REF] Hussein | Static Modeling of Sagging Cables with Flexural Rigidity and Shear Forces[END_REF][START_REF] Nguyen | On the simplifications of cable model in static analysis of large-dimension cable-driven parallel robots[END_REF][START_REF] Merlet | Some properties of the Irvine cable model and their use for the kinematic analysis of cable-driven parallel robots[END_REF]. In contrast, if cable tensions are sufficiently high and cable specific mass is relatively low, cables may also be considered as springs [START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF][START_REF] Mohammad | Dynamic modeling and control of parallel robots with elastic cables: Singular perturbation approach[END_REF][START_REF] Jamshidifar | Vibration Decoupled Modeling and Robust Control of Redundant Cable-Driven Parallel Robots[END_REF].

At any rate, one may note that if a model considering the cable elasticity and/or distributed mass is used, the static and dynamic loads applied on the platform result in different outputs of the kinematic model. In other words, the pose does not depend only on cable lengths, but also on the loads applied on the platform. The kinematic model is, therefore, coupled to the dynamic model. Whereas cable sagging and elasticity may be very difficult to integrate in the model of the robot, the platform motion for known cable forces may be modeled with a rather simple approach. As a matter of fact, for a known set of cable forces, the platform motion can be modeled using the Euler's equations of motion [START_REF] Goldstein | Classical mechanics[END_REF]Section 5.5].

Apart from the modeling of the cables, pulley kinematics and platform motion, a particular attention was dedicated to the forward kinematics. As for parallel robots in general, the estimation of the platform pose for given positions of the actuators is not a trivial task. Whereas this issue may be avoided with the application of an external measurement system capturing the pose of the platform. This can be done, for instance, with the motion-capture system Bonita developed by Vicon [START_REF] Chellal | Model identification and visionbased H-infinity position control of 6-DoF cable-driven parallel robots[END_REF], with a multi-cameras setup [START_REF] Dallej | Modeling and vision-based control of large-dimension cable-driven parallel robots using a multiple-camera setup[END_REF] or with a single camera embedded on the mobile platform [START_REF] Zake | Vision-Based Control and Stability Analysis of a Cable-Driven Parallel Robot[END_REF]. Nevertheless, the implementation of such measurement systems within uncontrolled environments may lead to various issues related to occlusions and lighting variations. Furthermore, these solutions also depend on additional sensors, which may be costly.

Therefore, several studies sought a reliable and accurate algorithm able to compute the CDPR platform pose for a given set of winch positions, i.e., Forward Kinematic (FK) Algorithms. The FK problem may be solved with interval analysis [START_REF] Pierre | Solving the forward kinematics of a Gough-type parallel manipulator with interval analysis[END_REF][START_REF] Berti | Solving the direct geometrico-static problem of underconstrained cable-driven parallel robots by interval analysis[END_REF]. This approach allows one to find all the possible FK problem solutions, taking into account that some cables may be slack. Even though this method is considered mathematically complex, real-time capability has been presented [START_REF] Pierre | On the real-time calculation of the forward kinematics of suspended cable-driven parallel robots[END_REF]. Additionally, iterative numerical schemes constitute an alternative typical solution [START_REF] Bruckmann | Wire Robots Part I: Kinematics, Analysis & Design[END_REF][START_REF] Miermeister | Differential kinematics for calibration, system investigation, and force based forward kinematics of cable-driven parallel robots[END_REF][START_REF] Pott | An algorithm for real-time forward kinematics of cable-driven parallel robots[END_REF]. In general, the FK problem is formulated as a nonlinear optimization problem. The iterative algorithm (usually, Levenberg-Marquadt) delivers a platform pose that locally minimizes modeling errors, taking the winch positions as input. To this end, a differential kinematic model is necessary.

Position Tracking Control

As discussed in the previous sections, this thesis is focused on advanced control techniques for CDPRs. The present section introduces the state-of-the-art of the studies addressing this topic. The main subject addressed in this thesis is the problem of position tracking control. The fundamental input of this problem is a desired trajectory. The control scheme is responsible for tracking this desired trajectory, i.e. generate a motion so that, at each instant of the movement, the measured (or estimated) actual pose is as close as possible to the desired one.

The approaches presented in the literature of position tracking control for CDPRs may be organized as follows. A commonly used classification separates the control schemes in two groups, namely, Joint Space Control (JSC) and Operational Space Control (OSC). In this thesis, a control strategy is classified as an OSC scheme if it measures (or estimates) the mobile platform pose error. Otherwise, the strategy is classified as JSC. Loosely speaking, an OSC strategy presents a closed loop which handles the platform pose error.

Joint space control

Most probably, the simpler schemes to first operate a CDPR were JSC strategies. Consider an actuation in which the motor angular positions are controlled taking the values obtained using the IK model for each pose of the desired trajectory. Such control scheme is commonly referred to as kinematic control and was the control scheme used to control the NIST ROBOCRANE [START_REF] Albus | The NIST robocrane[END_REF] shown in the Figure 1.8.

As the reader may note, the JSC approach may certainly lead to a simple CDPR control strategy. Nevertheless, more elaborate strategies may be derived with the same

J + + + + + + + + + + d q i d I m F v xd xd K p K v q q d C T d Figure 1
.9 -A Joint Space Control scheme proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]. The main variables are: The desired and actual pose x d and x, desired and measured motor positions q d and q, Jacobian matrix J, desired cable tensions τ d , dynamic motor friction constant F v , motor inertias I m , motor currents i, a diagonal constant matrix C T and PD gains K p and K v . The blocks TD and IK stand for the Tension Distribution algorithm and Inverse Kinematics, respectively. Time derivatives are denoted using the dot notation.

basis. Kawamura et al. propose a similar scheme adding the compensation of the gravitational effects and the management of internal forces [START_REF] Kawamura | High-speed manipulation by using parallel wire-driven robots[END_REF]. Fang et al. propose in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF] the control strategy of Figure 1.9. The closed-loop in this block diagram is exclusively dedicated to the joint space and the Inverse Dynamics block delivers the feedforward wrench which generates the desired motion.

These control strategies implement uncoupled feedback gains, i.e. the feedback gain matrices are diagonal. Alternatively, the synchronous tracking control approach introduced in [START_REF] Sun | Synchronous tracking control of parallel manipulators using cross-coupling approach[END_REF] may be used in order to improve the coordination between the cable displacements [START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF]. In practice, this strategy is based on gain matrices that are not diagonal.

+ + + + + + + + + + d q i d I m F v F s xd xd K p K i C T K v q q d d
Figure 1.10 -The Operational Space Control proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]. The variables that were not introduced in Figure 1.9 are: The pseudo-inverse jacobian matrix J + , static motor friction constant F s , tracking error of the pose ∆x and PID gains on operational space K p , K i and K v .

Kraus et al. [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF][START_REF] Kraus | Force control of cable-driven parallel robots[END_REF] propose a control scheme in which the cable tensions are controlled using the measurement of force sensors. As in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF], the desired wrench to be applied on the platform (and the corresponding cable tensions) are computed based on the desired trajectory without the estimation of the error of the platform pose.

Operational space control

As a matter of fact, the majority of the advanced control techniques for the position tracking of CDPRs belongs to the class of OSC. The strategy implemented by Lamaury and Gouttefarde in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] and presented in Figure 1.10 may be taken as an example to introduce the rationale commonly used on OSC. Comparing Figures 1.9 and 1.10, one may note that the schemes proposed in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF] and [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] are rather similar. The main difference lies on the fact that instead of applying the closed-loop strategy directly on the difference between desired and measured motor positions, the PID correction in Figure 1.10 takes as input the estimated error of the platform pose (an estimation obtained as ∆x = J + (q dq)).

Inasmuch as a large amount of studies present a very similar approach, it is worthwhile to highlight some details involved in its motivations. One may note that the feedforward terms compensate most of the nonlinearities of the system, rendering the proposed approach closely related to the widely used computed torque control [START_REF] Richard | Robot manipulators: mathematics, programming, and control: the computer control of robot manipulators[END_REF] or the feedback linearization method [START_REF] Hassan | Nonlinear control[END_REF]Section 9.3]. Indeed, most of the control strategies discussed in this section presents a very similar rationale. The corresponding schemes are alike the one in Figure 1.10, changing notably the highlighted region, which computes the wrench representing the feedback correction for a given error of the platform pose. Note that full feedback linearization may be easily applied for nonredundant CDPRs, as proposed in [START_REF] Habibnejad Korayem | Dynamic load carrying capacity of flexible cable suspended robot: Robust feedback linearization control approach[END_REF]. Nevertheless, the presence of a number of cables greater than the number of degrees of freedom may prevent the application of such method for redundant CDPRs.

Whereas the nominal stability of these approaches can be straightforwardly proven given the simplicity of the nominal closed-loop system, interesting results were ob-tained regarding robustness. The robust stability considering structured and unstructured uncertainties of a rather similar scheme was proved in [START_REF] Mohammad | Robust PID control of fullyconstrained cable driven parallel robots[END_REF] using Lyapunov theory. Similarly, Korayem et al. analyzed the robustness regarding the flexibility of the joints of a fully feedback linearized control scheme [START_REF] Habibnejad Korayem | Dynamic load carrying capacity of flexible cable suspended robot: Robust feedback linearization control approach[END_REF]. Alternatively to the direct Lyapunov theory, Khosravi and Taghirad addressed in [START_REF] Mohammad | Dynamic modeling and control of parallel robots with elastic cables: Singular perturbation approach[END_REF] the stability of a similar scheme using singular perturbation theory [START_REF] Kokotovic | Singular perturbation methods in control: analysis and design[END_REF]. A simple and elegant solution was proposed to extend the stability results when the elasticity of the cables are considered.

Yet in the field of linear control systems, more advanced techniques may be implemented. The Linear Quadratic Regulator (LQR) method was used in [START_REF] Habibnejad | Maximum DLCC of spatial cable robot for a predefined trajectory within the workspace using closed loop optimal control approach[END_REF] considering the feedback linearized system. As for the study in [START_REF] Habibnejad Korayem | Dynamic load carrying capacity of flexible cable suspended robot: Robust feedback linearization control approach[END_REF], a non-redundant robot is necessary to obtain a linear system. The LQR method delivers the optimal gain matrices corresponding to a deterministic infinite horizon optimal control problem applied for linear systems [START_REF] Skogestad | Multivariable feedback control: analysis and design[END_REF]Section 9.2.1]. Due to the fact that its theoretical basis relies on the deterministic model of the system, the LQR may be sensitive to disturbances. To deal with this issue, the H ∞ (H-infinity) robust control [START_REF] Skogestad | Multivariable feedback control: analysis and design[END_REF]Section 9.3.4] method may be used. This possibility was investigated in [START_REF] Chellal | Model identification and visionbased 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] using the Inca robot (as the one in Figure 1.3). Motivated by the fact that all the aforementioned studies present control schemes that are based on the theory of linear systems, such methods are denoted in this section as linear-based methods.

Clearly, methods developed with theories that are not based on the linearization of the studied system may be used. These approaches are classified in this section as nonlinear-based methods. Firstly, Sliding Mode Control (SMC) received particular attention in the control of CDPRs. The main advantages of this nonlinear control method are finite time convergence and robustness to modeling uncertainties and disturbances. Moreover, the implementation of SMCs are often simpler than other advanced nonlinear control techniques.

Oh and Agrawal [START_REF] Ryeok | Nonlinear sliding mode control and feasible workspace analysis for a cable suspended robot with input constraints and disturbances[END_REF] implemented the theory proposed in [START_REF] Slotine | Sliding controller design for non-linear systems[END_REF] in order to estimate the feasible workspace of CDPRs. Subsequently, improved versions were proposed, e.g. 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 | Robust adaptive sliding mode control of a redundant cable driven parallel robot[END_REF]. Nevertheless, chattering is a well-known issue commonly addressed in studies implementing SMCs. Even though numerous methods were proposed in order to reduce this difficulty, chattering is still observed in recent results [START_REF] Schenk | Robust adaptive sliding mode control of a redundant cable driven parallel robot[END_REF]. Besides SMC, other nonlinear feedback control approaches were proposed, such as fuzzy control [START_REF] Zi | Dynamic modeling and active control of a cable-suspended parallel robot[END_REF], neural networks [START_REF] Jabbari | Adaptive neural network control of cable-driven parallel robots with input saturation[END_REF] and the intricate method proposed in [START_REF] Li | Tracking Control of Fully-Constrained Cable-Driven Parallel Robots using Adaptive Dynamic Programming[END_REF] involving dynamic programming with an actorcritic structure.

Inspired by the fact that the applications of CDPRs are frequently subjected to nonnegligible uncertainties (in particular regarding the payload and the friction of the winches), a lot of effort has been devoted to adaptive strategies. Most of them are based on the adaptation law introduced in [START_REF] Weiping | On The Adaptive Control of Robot Manipulators[END_REF], also proposed by Slotine. An adaptive version of [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF] is presented in [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF], resulting in a better tracking performance than its predecessor. As a matter of fact, many of the aforementioned nonlinear control schemes are classified as adaptive strategies [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF][START_REF] Schenk | Robust adaptive sliding mode control of a redundant cable driven parallel robot[END_REF][START_REF] Jabbari | Adaptive neural network control of cable-driven parallel robots with input saturation[END_REF][START_REF] Li | Tracking Control of Fully-Constrained Cable-Driven Parallel Robots using Adaptive Dynamic Programming[END_REF]. Figure 1.11 summarizes the classifications discussed in the present subsection.

One may note that most of the schemes discussed in the present section define,

Joint Space

Position Tracking of CDPRs

State of the Art [72]

[74]

[75]

[73]

[69]

• Neural Networks Asl

• Dynamic Programming Li [7]

[59]

[6]

[64]

[8]

[47], [START_REF] Laroche | A preliminary study for H infinity control of parallel cable-driven manipulators[END_REF] [63]

[66] at some point, a desired wrench to be applied on the platform. This wrench should be generated by the combination of the forces applied by the cables on the mobile platform. For a non-redundant CDPR, i.e. having a number of cables equal to the number of DoF, the relation between these sets of variables may be straightforwardly obtained. Nevertheless, for a redundantly actuated CDPR, this subject merits special attention. The next subsection is devoted to this issue.

Tension Distribution Algorithms

Different kinds of redundancies may be defined in robotics [START_REF] Müller | On the terminology and geometric aspects of redundant parallel manipulators[END_REF][START_REF] Sahin | Clarifying the definition of redundancy as used in robotics[END_REF][START_REF] Pierre | Parallel robots[END_REF]. In the context of this thesis, a specific type is of particular interest, namely, actuation redundancy. Although the definition of this class of redundancy may be the subject of discussions [START_REF] Müller | On the terminology and geometric aspects of redundant parallel manipulators[END_REF], a simple definition may be synthesized for the CDPRs studied in this thesis. The classification of redundancy for reconfigurable CDPRs [START_REF] Gagliardini | Discrete reconfiguration planning for Cable-Driven Parallel Robots[END_REF][START_REF] Abdolshah | Optimizing Stiffness and Dexterity of Planar Adaptive Cable-Driven Parallel Robots[END_REF][START_REF] Gagliardini | A reconfiguration strategy for Reconfigurable Cable-Driven Parallel Robots[END_REF][START_REF] Nguyen | On the analysis of large-dimension reconfigurable suspended cable-driven parallel robots[END_REF] or robots presenting articulated mobile platform [START_REF] Miermeister | Design of Cable-Driven Parallel Robots with Multiple Platforms and Endless Rotating Axes[END_REF][START_REF] Reichenbach | Static Analysis of a Two-Platform Planar Cable-Driven Parallel Robot with Unlimited Rotation[END_REF] may render necessary a detailed discussion. Nevertheless, the studies presented here are focused on non-reconfigurable robots with a platform consisting of a unique non-articulated body. In this case, a CDPR having a number of cables (m) greater than the number of DoF of the mobile platform (n) is classified as redundantly actuated. Most of the results presented in this thesis takes as motivation some particularities of this class of CDPRs.

For the majority of commonly used modeling approaches, for a given platform pose, one may define a constant linear mapping between the vector of cable tensions τ and the resulting wrench f applied on the platform [START_REF] Bosscher | Wrench-based analysis of cable-driven robots[END_REF]. Here, each element τ i of the vector τ is the euclidean norm of the (vectorial) force applied on the platform by the corresponding cable. More precisely, there exists a matrix W, named wrench matrix, so that f = W τ . Therefore, if m = n, a set of cable tensions able to generate a desired wrench f d may be trivially obtained inverting such (bijective) linear mapping (if the considered pose is not singular, leading to an invertible square matrix W). Otherwise, non-trivial problems emerges if m = n. Robots with m < n are referred to as under-constrained CDPRs [START_REF] Carricato | Stability analysis of underconstrained cable-driven parallel robots[END_REF][START_REF] Abbasnejad | Direct Geometrico-static Problem of Underconstrained Cable-Driven Parallel Robots With n Cables[END_REF] and their study is out of the scope of this thesis. In contrast, if m > n, infinitely many sets of cable tensions may generate the same wrench. The algorithms responsible for the definition of the most pertinent set of cable tensions for a given desired wrench are called Tension Distribution Algorithms (TDA).

In general, TDAs should comply with two sets of constraints. Clearly, one set of constraints result from the fact that the cable tensions should generate the desired wrench, i.e. vector τ should satisfy f d = W τ . In addition, cable tensions should be within a positive interval. This positive interval is delimited by sets of maximal and minimal cable tensions, τ max and τ min , respectively. The minimal cable tensions are considered in order to avoid excessively low tensions leading to the malfunctioning of the winches and to slack cables, which may result in deteriorated precision and stiffness. Maximal cable tensions are defined according to the maximal load with under which every mechanical part operates safely.

To sum up, the goal of a TDA is to define a pertinent set of cable tensions τ so that f d = W τ and τ min τ τ max , taking as input the linear mapping W, the desired wrench f d and the cable tension limits τ min and τ max . Accordingly, a TDA can be formulated as an optimization problem

τ * = arg min τ f (τ ) , s. t. f d = W τ τ min τ τ max (1.1)
with a given cost function f : R m → R, so that the essential difference among them is the definition of f .

A natural choice for f is some norm of the vector τ . Using the typical definition of vector norms as in [START_REF] Gene | Matrix computations[END_REF]Section 2.2]

with τ p = (τ p 1 + τ p 2 + • • • + τ p m ) 1/p
, several studies address the TDA formulated as (1.1) for f (τ ) = τ p with p some positive non-null integer.

Using the 1-norm with p = 1, the corresponding optimization problem (1.1) is classified as a Linear Programming (LP) problem [START_REF] Bertsimas | Introduction to linear optimization[END_REF]. Snyman and Hay simplified the general methodology implemented to the solution of LP problems for the specific constraints of TDAs [START_REF] Snyman | Analysis and Optimization of a Planar Tendon-Driven Parallel Manipulator[END_REF]. Similarly, an explicit expression of the vector of cable tensions minimizing τ 1 = τ i is deduced in [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF]. As an alternative this, one may apply weighing costs c i for each individual cable i according to the specific needs of the CDPR. For instance, a critical issue for the suspended CDPR in [START_REF] Shiang | Optimal force distribution applied to a robotic crane with flexible cables[END_REF] is to keep the longer cables sufficiently tighten. Hence, a cost function f (τ ) = c i τ i with negative c i may be used for the cables that need to be under increased cable tensions. Besides the case of f (τ ) = c i τ i , the LP methodology may be used to minimize the maximum element of the vector τ (also known as minimax problem). As a matter of fact, any minimax problem with linear constraints may be solved through an equivalent LP problem [START_REF] Ahuja | Minimax linear programming problem[END_REF]. Based on this approach, a versatile method is proposed in [START_REF] Per | Rapid computation of optimally safe tension distributions for parallel cable-driven robots[END_REF] to minimize f (τ ) = max τ i . Since max i τ i = lim p→∞ τ p , the notation τ ∞ := max i τ i is commonly used. Moreover, the proposed method is versatile in the sense that it is able to maximize the minimal cable tension. Nevertheless, as highlighted in this study, TDAs based on norms with p = 1 or p = ∞ may lead to discontinuities in the profiles of cables tensions, i.e. infinitesimal changes in f d may lead to finite changes of the τ * obtained with the corresponding TDAs. Clearly, additional constraints may be taken into account in the optimization problem (1.1) so that the variation of cable tension remains limited [START_REF] Shiang | Optimal force distribution applied to a robotic crane with flexible cables[END_REF].

In contrast, methods taking f (τ ) = τ p with 1 < p < ∞ (with strict inequalities) intrinsically lead to continuous profiles of cable tensions. Gosselin presents an overview on this matter in [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF] emphasizing this advantage. On the other hand, assuming that the presence of cable tensions in the neighborhood of the bounds τ min and τ max should be avoided, he points out that it may be preferable to minimize

f (τ ) = τ -τ ∞ with τ = 1 /2 (τ min + τ max ).
It is shown that this approach, compared to the application of other values of p, tends to maximize the difference between the desired cable tensions and the cable tension limits, preventing the incidence of values close to τ min and τ max (see Figure 1.12). Inspired by this advantage, p = 4 is proposed in order to ally continuity of cable tensions and increased difference between τ * and bounds τ min and τ max . Figure 1.13 illustrates that this value of p may be used to this end. Moreover, an explicit expression for the TDA taking m = n + 1 and f (τ ) = ττ 4 is presented in [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF].

In spite of the advantages highlighted in [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF] led by the use of higher values of p, the majority of studies addressing TDAs with the format (1.1) is focused on values of p 2. Besides the above-mentioned works addressing the case in which p = 1, particular attention has been devoted to the problem of TDAs based on the 2-norm of the vector of cable tensions [START_REF] Bruckmann | Calculating force distributions for redundantly actuated tendon-based Stewart platforms[END_REF][START_REF] Bruckmann | A Modular Controller for Redundantly Actuated Tendon-Based Stewart Platforms[END_REF][START_REF] Martin | Determination and management of cable interferences between two 6-DOF foot platforms in a cabledriven locomotion interface[END_REF][START_REF] Hassan | Analysis of bounded cable tensions in cable-actuated parallel manipulators[END_REF][START_REF] Hamid | An analytic-iterative redundancy resolution scheme for cable-driven redundant parallel manipulators[END_REF]. In addition to the fact that such algorithms intrinsically delivers continuous profiles of cable tensions, the resulting optimization problems are strictly convex Quadratic Programming (QP) problems [START_REF] Nocedal | Numerical optimization[END_REF]. The numerical tractability of this class of optimization problem is an important advantage. Accordingly, the well-known numerical tools [START_REF] Demmel | LAPACK: A portable linear algebra library for supercomputers[END_REF][START_REF] Ford | The NAG library 'machine[END_REF] able to solve QP problems were applied to TDAs minimizing the 2-norm of cable tensions [START_REF] Bruckmann | Calculating force distributions for redundantly actuated tendon-based Stewart platforms[END_REF][START_REF] Bruckmann | A Modular Controller for Redundantly Actuated Tendon-Based Stewart Platforms[END_REF]. Alternatively, the Dykstra's Projections Algorithm [START_REF] Hassan | Analysis of bounded cable tensions in cable-actuated parallel manipulators[END_REF] or direct consequences of the Karush-Kuhn-Tucker (KKT) conditions [START_REF] Kuhn | Nonlinear Programming[END_REF] may be used [START_REF] Hamid | An analytic-iterative redundancy resolution scheme for cable-driven redundant parallel manipulators[END_REF]. All these methods [START_REF] Bruckmann | Calculating force distributions for redundantly actuated tendon-based Stewart platforms[END_REF][START_REF] Bruckmann | A Modular Controller for Redundantly Actuated Tendon-Based Stewart Platforms[END_REF][START_REF] Martin | Determination and management of cable interferences between two 6-DOF foot platforms in a cabledriven locomotion interface[END_REF][START_REF] Hassan | Analysis of bounded cable tensions in cable-actuated parallel manipulators[END_REF][START_REF] Hamid | An analytic-iterative redundancy resolution scheme for cable-driven redundant parallel manipulators[END_REF] are iterative and, therefore, yield a performance highly dependent on the initial guess. Fortunately, since TDAs with 2-norm cost functions deliver continuous cable tensions for continuous variations of desired wrenches, consecutive time samplings should deliver solutions close to each other. This evidence justifies the use of the solution obtained in a sampling time as initial guess for the next sampling time (the so-called warm start). Moreover, as discussed in [START_REF] Bruckmann | Calculating force distributions for redundantly actuated tendon-based Stewart platforms[END_REF], gradient-based algorithms such as the gradient projection method may be used with a limited number of iterations, providing feasible suboptimal solutions if realtime constraints hinder the computation of the global optimal solution.

Even though TDAs may be used in different circumstances, such as the design of CDPRs [START_REF] Joshua T Bryson | Optimal Design of Cable-Driven Manipulators Using Particle Swarm Optimization[END_REF], it is reasonable to assert that the main application of such algorithms is the real-time control. Therefore, the computation time is a major concern. Pott et al. proposes in [START_REF] Pott | Closed-form Force Distribution for Parallel Wire Robots[END_REF] a closed-form expression able to minimize ττ 2 = ττ . Besides the short computation time, this method presents an attractive simplicity of implementation. Additionally, in contrast to iterative methods, the closed-form solution demands well defined floating point operations, which may be used to bound the worst case of the computational burden. Nevertheless, even though an improved algorithm was introduced in [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF], these methods do not have a clear proof that a feasible distribution may be found for any feasible desired wrench.

Alternatively to norms of the vector of cable tensions, indirect consequences of the changes in τ may be addressed as well. In particular, changes in the internal forces may result in variations on the stiffness of the CDPR [START_REF] Lim | Optimization of tension distribution for cable-driven manipulators using tension-level index[END_REF][START_REF] Jamshidifar | Kinematically-Constrained Redundant Cable-Driven Parallel Robots: Modeling, Redundancy Analysis, and Stiffness Optimization[END_REF][START_REF] Yu | Simultaneous Trajectory Tracking and Stiffness Control of Cable Actuated Parallel Manipulator[END_REF]. The influence of the cable tensions on the robot stiffness is studied under different approaches [START_REF] Azadi | Antagonistic variable stiffness elements[END_REF][START_REF] Yuan | Static and dynamic stiffness analyses of cable-driven parallel robots with non-negligible cable mass and elasticity[END_REF]. Probably, the most commonly applied approach considers cables as springs with a profile consisting of a straight line [START_REF] Azadi | Antagonistic variable stiffness elements[END_REF]. Using this model, the stiffness along a given direction may be optimized as proposed in [START_REF] Jamshidifar | Kinematically-Constrained Redundant Cable-Driven Parallel Robots: Modeling, Redundancy Analysis, and Stiffness Optimization[END_REF].

Moreover, other methods not based on the formulation (1.1) were also proposed. Many studies are based on the characterization of the set formed by the constraints in (1.1). This is the set of feasible cable tensions that generate the desired wrench. As discussed in [START_REF] Bouchard | On the Ability of a Cable-Driven Robot to Generate a Prescribed Set of Wrenches[END_REF][START_REF] Gouttefarde | Characterization of Parallel Manipulator Available Wrench Set Facets[END_REF], this set is a polytope which is a subset of R m . As any polytope, it may be described by its vertices. Aiming at the avoidance of cable tensions close to tension limits, Mikelsons et al. proposes the determination of the barycenter of such polytope in [START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF]. The method is proved to be continuous, but the computation of the polytope itself is not discussed in depth. Accordingly, Gouttefarde et al. introduces a method in [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n + 2 Cables[END_REF] (with preliminary study published in [START_REF] Lamaury | A Tension Distribution Method with Improved Computational Efficiency[END_REF]) able to compute this set for the specific case of m = n+2. This approach proved to be versatile in the sense that several TDAs may be adopted based on the output delivered by such algorithm. Explicit computation of the following outputs are presented: centroid, weighed barycenter, 1norm and 2-norm optimal solutions. The algorithm finds the values of λ ∈ R 2 so that τ = W + f d + N λ are the vertices of the polytope of the TDA solutions, with N ∈ R m×2 the null space matrix of W. Therefore, the identification of the polytope vertices lying on R 8 is simplified to computations in the 2-dimensional space of λ.

Cable Tension Control

Inasmuch as the output of the TDAs is a set of desired cable tensions, a strategy should be implemented to convert desired cable tensions to a variable that is meaningful to the motor drivers, typically motor torques (or currents), positions or velocities. Such strategies are denoted in this thesis as cable tension control schemes. Compared to the position tracking control, this problem received less attention among the research studies related to the control of CDPRs. Notably, many works presenting exclusively simulations focused on the position tracking control neglect this subject, considering that the desired cable tensions are instantly generated [START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF][START_REF] Schenk | Robust adaptive sliding mode control of a redundant cable driven parallel robot[END_REF].

Maybe due to the fact that the measurement of forces are often very noisy or because force sensors to measure cable tensions were not available, the majority of published experimental results were obtained with open-loop cable tension control schemes. Figure 1.14 presents a very simple scheme using motor torques as control input to the motor drivers 1 . In spite of its simplicity, this strategy has been widely used. Clearly, this is the case for the schemes in Figures 1.9 and 1.10 of works [START_REF] Fang | Motion control of a tendon-based parallel manipulator using optimal tension distribution[END_REF] and [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF], respectively. Additionally, such open-loop control is also used in other studies cited previously, e.g. [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF][START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF].

Alternatively, force sensors may be used enabling the implementation of a feedback control strategy. Kraus et al. used motor positions as control inputs to the motor q sign()

CDPR + + + + + d d q i d I m F v F s C T Figure 1.14 -A
simple open-loop cable tension control scheme [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF].

drivers [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF]. This work presents details on the identification of the dynamic system and satisfying experimental results on position tracking control. With a very similar cable tension control scheme, the same author implemented a hybrid force-position control [START_REF] Kraus | Hybrid Position/Force Control of a Cable-Driven Parallel Robot with Experimental Evaluation[END_REF]. Reichert et al. attains a similar goal using torques as control inputs of the motor drivers [START_REF] Reichert | Robust internal force-based impedance control for cable-driven parallel robots[END_REF]. Using a disturbance observer, a robust internal force-based impedance control is proposed [START_REF] Reichert | Robust internal force-based impedance control for cable-driven parallel robots[END_REF]. A more advanced technique is proposed and tested through numerical simulations in [START_REF] Meunier | Control of an overactuated cable-driven parallel mechanism for a radio telescope application[END_REF] using H ∞ robust control in order to compensate the influence of the cable elasticity.

The use of measured cable tensions as feedback signals may be of particular interest for fully-constrained CDPRs. Consider a suspended CDPR presenting a pose and payload in accordance with the designed nominal workspace and admissible payloads. Excessively high cable tensions necessarily lead to high accelerations. Therefore, safety of suspended CDPRs regarding the limits of cable tension may be indirectly monitored by other aspects, namely, pose, payload, velocity and accelerations. For this reason, some suspended CDPRs do not present force sensors (even though most of the cranes directly monitor cable tension for safety purposes). In contrast, fully-constrained CD-PRs may develop excessively high cable tensions without any other apparent variation on these variables. Therefore, for safety reasons, the use of force sensors monitoring cable tensions of a fully-constrained CDPR is highly recommended. This is the case for the robots studied in this thesis. Accordingly, the implementation of an inexpensive (regarding the cost) cable tension control scheme taking advantage of the presence of the force sensors is pertinent.

Model Predictive Control

According to the studies presented in Section 1.5, the minimal and maximal cable tensions play a major role in the control of CDPRs. However, one may note that the cable tension limits, in spite of their critical importance, are generally taken into account (in the TDA) once the desired wrench is already computed. In other words, the wrench feasibility cannot be directly considered in the corresponding control approaches. In practice, if a non-feasible wrench is computed by any of the OSC modules [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]- [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF], the corresponding TDA will not be able to compute cable tensions, resulting in a sudden interruption of the robot operation. Similarly, the discussed JSC schemes handle the tension limits uniquely in the feedforward terms, as in the scheme in Figure 1.9.

Therefore, one cannot guarantee that the final cable tension is within the feasible limits.

Typically, a non-feasible desired wrench may be necessary in case of excessive disturbances or if the desired pose is out of the actual workspace. Even though the robot workspace can be theoretically computed beforehand, so that the feasible trajectories may be analyzed out of the real-time operation, the proposed workspace computation methods were proven to deliver not more than estimative results. Depending on the considered assumptions and all kind of uncertainties including approximate knowledge of the CDPR geometric parameters, the dimensions of the estimated workspace may change considerably and the computation itself may be troublesome [START_REF] Pierre | On the workspace of suspended cable-driven parallel robots[END_REF]120]. Moreover, changes on the values and location of the payload lead to variations in the workspace dimensions.

As discussed in Section 1.2, the robot designed in the Hephaestus project should operate in a construction site so that the presence of disturbances is a major concern and safety deserves particular attention. In addition, the limited available area in the neighborhood of the building results in a workspace with small margin for uncertainties, i.e. the difference between the workspace predicted in the robot design and the necessary to perform the tasks is very narrow. In this context, the operation of the Hephaestus prototype with a control scheme similar to state of the art approaches [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]- [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF] would be probably subjected to the sudden interruptions described above, which may lead to unacceptable safety issues. In practice, a control scheme able to resolve this problem should take into account the cable tension limits within the computation of the desired wrench.

Model Predictive Control (MPC) is one of the few control methodologies able to explicitly handle system constraints [START_REF] Marian | Predictive control: with constraints[END_REF][START_REF] Camacho | Model predictive control[END_REF]. The main idea of MPC is well described in the early work of Lee and Markus [START_REF] Bruce | Foundations of optimal control theory[END_REF]:

One technique for obtaining a feedback controller synthesis from knowledge of open-loop controllers is to measure the current control process state and then compute very rapidly for the open-loop control function. The first portion of this function is then used during a short time interval, after which a new measurement of the process state is made and a new open-loop control function is computed for this new measurement. The procedure is then repeated.

In other words, an open loop optimal control problem is repeatedly solved based on updated measured states and the first portion of the optimal sequence is applied during one control cycle time. The procedure is repeated at each control cycle. As better discussed in the course of the present section, this optimal control problem may explicitly handle the system constraints.

Although the main idea of MPC is remarkably simple, its implementation relies on the solution of an optimal control problem in real-time. Clearly, considering the numerous possible combinations of system dynamics and strategies to solve the corresponding optimal control problem, the class of MPC schemes includes a variety of control strategies with characteristics remarkably different. This section tries to summarize these main characteristics and evaluate their pertinence in the context of this thesis.

f(t) v(t) x(t) t = 0 x d

Optimal control

According to Vinter [START_REF] Vinter | Optimal control[END_REF], Optimal Control Problems (OCPs) emerged in the 1950s taking the aerospace engineering field of the American and Russian space programs as important source of motivating problems. As an illustrative example, take the very simple system of Figure 1.15. The motion of a body with unitary mass within a onedimensional space is studied in a time span t ∈ [0, t f ]. The body departs from position x(0) = 0 with null initial velocity v(0) = 0 and is actuated by a force, namely, f (t). The motion satisfies the Newton's Second Law v(t) = f (t) and the system is constrained with respect to the applied force and velocity according to |f (t)| f max and v(t) v max for all t ∈ [0, t f ]. A twofold goal is defined: The final position x(t f ) should be as close as possible to a desired position x d and the squared final velocity v(t f ) 2 should be minimized as well. Accordingly, the following OCP may be formulated:

min f x(t f ) -x d 2 + v(t f ) 2 . s. t. x(0) = 0, v(0) = 0 ẋ(t) = v(t), v(t) = f (t) |v(t)| v max , |f (t)| f max , ∀ t ∈ [0, t f ] (1.2)
In spite of its simplicity, this illustrative example introduces the main elements of a general OCP, which may be summarized as follows:

Present example

• A dynamic system v(t) = f (t) • A time span (bounded or unbounded) t ∈ [0, t f ] • System constraints |v(t)| v max , |f (t)| f max • A cost functional x(t f ) -x d 2 + v(t f ) 2
It is worth noting that, as presented in the list above, the analyzed time span may be bounded or unbounded. In contrast to this illustrative example that implements a finite time span, the interval of time t ∈ [t 0 , ∞) for a given initial time t 0 is used in OCPs with an unbounded time span. The theory of infinite horizon optimal control [START_REF] Dean A Carlson | Infinite horizon optimal control: deterministic and stochastic systems[END_REF] addresses this class of OCP. The (LQR) and Linear Quadratic Gaussian (LQG) controllers are well-known methods using this approach.

Even though OCPs based on continuous time systems may be categorized within the class of optimization problems, it is worth highlighting that they present a rather particular characteristic when compared to the classical problems addressed in optimization theory. Instead of taking as argument a finite set of numbers, a continuous OCP take as argument a vector function, namely, the control inputs within the studied time span (f (t) in the example above). This is the reason why the term cost functional (instead of cost function) was used in the list of the main elements of an OCP. As defined in [START_REF] Gelfand | Calculus of variations[END_REF], a functional is a "correspondence which assigns a definite (real) number to each function (or curve) belonging to some class". Briefly speaking, the main goal of a continuous OCP is to minimize (or maximize) a cost functional, which quantifies how good is the obtained behavior. The solution of the problem represents the optimal profile of control inputs.

Once the problem is stated, a method to compute the control inputs as functions of time should be defined. To this end, the variations of the control inputs should be mapped to the variations of the cost functional. This subject is addressed by means of the Variational Calculus, which constitutes the theoretical basis to the solution of classical problems such as the derivation of the seminal Euler-Lagrange equation [START_REF] Gelfand | Calculus of variations[END_REF]Section 1.4]. More precisely, the Euler-Lagrange equations are differential equations that are necessary conditions to the optimality of a functional. In order to cope with the solution of a general optimal control problem (in particular, the system constraints), the results of the Variational Calculus were extended by Pontryagin [START_REF] V G Boltyanskiy | Mathematical theory of optimal processes[END_REF], who introduced the maximum principle.

As an alternative to these methods, the continuous system may be studied with a discrete-time approach. In this case, continuous-time control inputs and states are converted to sequences of control inputs and states. This numerical discretization may be performed using the methods proposed in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 9] and leads to a discrete-time system, which is an approximation of the original continuous system.

As an illustration, considering the example introduced earlier, the functions representing the force f (t) and states s(t) for t ∈ [0, t f ] are converted to finite sequences {f 0 , ..., f n f } and {s 0 , s 1 , ..., s n f +1 }, with s k = x k v k T containing the actual position and velocity, x k and v k respectively. In this discrete-time version (in the sense of [128, Section 2.3]) of the aforementioned example, the sequence of states {s 0 , s 1 , ..., s n f +1 } and the sequence of control inputs {f 0 , f 1 , ..., f n f } satisfy a transition mapping

s k+1 = φ(s k , f k ), ∀ k ∈ {0, 1, ..., n f }, (1.3)
and the OCP (1.2) may be reformulated as

min f 0 ,...,fn f x n f +1 -x d 2 + v 2 n f +1 . s. t. s 0 = 0 s k+1 = φ(s k , f k ) ∀ k ∈ {0, ..., n f } |v k | v max , |f k | f max , ∀ k ∈ {0, ..., n f } (1.4)
This discrete counterpart of the OCP discussed earlier takes as argument a sequence of control inputs, instead of a continuous function. Therefore, the resulting optimization problem may be solved with standard optimization techniques. Furthermore, since the vast majority of control systems used nowadays are digital, the discrete-time approach of OCP led to relevant influences on the current state of the art of the MPC. For this reason, this subject is discussed here more in detail.

The Dynamic Programming Principle (DPP), introduced by Bellman [START_REF] Bellman | Dynamic Programming[END_REF], is able to break some optimization problems in several smaller subproblems in a recursive manner. As illustrated in Figure 1.16, the DPP is useful when decisions should be taken successively. This figure introduces the main elements of a problem addressed by this principle. The possible states are distributed in several layers, and successive layers present their states connected by transitions with corresponding costs. The problem of finding the sequence of transitions linking the initial state to the final state with a minimal accumulated cost (optimal path) may be divided into several subproblems in which a cost-to-go is attributed to each state of a layer. For a given state s, the cost-to-go is the minimal cost to attain the final state departing from the state s.

More precisely, the DPP may be sketched as follows. Denote s i the possible states that can be reached directly from s, and c i the stage cost corresponding to the transition from s to s i , with i ∈ {1, ..., n s } and n s the number of possible states reachable from s. Denote also V (s) the cost-to-go of any state s. The DPP states that

V (s) = min i c i + V (s i ) .
(1.5)

The applicability of the DPP in discrete-time OCPs is evident. However, taking a continuous system such as the dynamic model of a CDPR, one may note that the application of the DPP as it is involves two instances of discretization. Firstly, the system should be discrete-time, i.e. time should be considered as a discrete sequence t ∈ {t 0 , t 1 , ..., t f }, instead of a continuous interval. Secondly, the states and control inputs themselves should be discretized as well. That is to say, instead of considering continuous intervals of admissible cable tensions, platform poses and velocities, one should consider finite discrete sets. As a result, a mesh of possible states and control inputs need to be considered. As already pinpointed by Bellman [129, page 322], this second instance of discretization leads to a severe drawback for systems with high dimensions (this drawback is commonly referred to as the curse of dimensionality). Nevertheless, the DPP plays a very important role as an analysis tool for the performance and stability evaluation of MPC schemes [START_REF] Grüne | Nonlinear model predictive control[END_REF]. Moreover, the theory introduced in [START_REF] Bellman | Dynamic Programming[END_REF] inspired other MPC approaches that do not require the discretization of states and control inputs, such as the Differential Dynamic Programming (DDP) [START_REF] Sen | A quasi-newton differential dynamic programming algorithm for discrete-time optimal control[END_REF] which has been used in recent and advanced applications of MPC [START_REF] Tassa | Control-limited differential dynamic programming[END_REF].

Receding horizon model predictive control

Note that the discussions related to optimal control hitherto do not consider a closed-loop strategy. Recalling the discussion at the beginning of this section, the main idea involved in MPC is to recompute the optimal control inputs as frequently as possible taking the measured states into account, leading to a closed-loop algorithm (see the illustration in Figure 1.17). Consequently, an optimal control problem should be solved in each cycle of the feedback loop. Such as illustrated in Figure 1.17, this concept results in an interesting consequence on the analyzed time span of MPC schemes. First, consider a control scheme that solves repeatedly an OCP with a finite time span. In the context of MPC, the final instant of such time span is called prediction horizon. Obviously, the initial time of the OCP is shifted every time it is computed. Since the analyzed length of time is kept constant, the prediction horizon is also shifted. This feature presented by MPC schemes solving OCP with finite time spans inspires the commonly used term Receding Horizon Model Predictive Control. On the other hand, if infinite time span is used, the horizon is constant (t f = ∞). Therefore, the aforementioned examples LQR and LQG do not present receding horizons. For this reason, even though this class of feedback control repeatedly solves OCPs taking updated measures of states, they are rarely included within the set of MPC schemes. In brief, the commonly used classification indicates that a control scheme is considered as an MPC scheme only if it presents a receding horizon. Accordingly, the definition of MPC used in the present thesis is as in [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]:

Model predictive control (MPC) or receding horizon control (RHC) is a form of control in which the current control action is obtained by solving on-line, at each sampling instant, a finite horizon open-loop optimal control problem, using the current state of the plant as the initial state; the optimization yields an optimal control sequence and the first control in this sequence is applied.

In spite of the fact that earlier works, such as [START_REF] Propoı | Application of linear programming methods for the synthesis of automatic sampled-data systems[END_REF], introduced strategies matching with the aforementioned rationale, the origin of MPC is commonly attributed to Richalet. The work [START_REF] Richalet | Model predictive heuristic control[END_REF] coined the name Model Predictive Heuristic Control to the method in which impulse-response methodology is used to identify and model the controlled system. Promising experimental results were obtained for the control of a steam generator, a distillation column in an oil refinery and a complete chemical plant synthesizing vinyl chloride. Richalet already highlighted the prominent advantage that constraints may be handled explicitly within the feedback loop. Handling constraints is particularly relevant since the optimal operation in industrial processes are frequently attained on the boundary of the set of feasible states and control inputs. However, some strategies which cannot be neglected on the timeline of the development of MPC schemes do not consider system constraints, e.g. the Generalized Predictive Control (GPC) proposed by Clarke [START_REF] Clarke | Generalized predictive control-Part I. The basic algorithm[END_REF]. Since Richalet's work [START_REF] Richalet | Model predictive heuristic control[END_REF], MPC has been the focus of a great deal of attention both in academia and industry. Indeed, as advocated by Maciejowski [START_REF] Marian | Predictive control: with constraints[END_REF], MPC is sometimes known as "the only advanced control technique -that is, more advanced than standard PID -to have had a significant widespread impact on industrial process control". Here, it is interesting to highlight the scope of industrial process control. Besides the evident profitability obtained with a control scheme able to safely run close to constraints, industrial processes often present system dynamics much slower than other typical applications of advanced control strategies, such as robotics. Systems with slow dynamics enable the application of control running with longer cycle times. Therefore, issues related to the numerical solution of OCPs in real-time may be overcome. For these reasons, the vast majority of relevant applications of MPC in the 20th century was focused on industrial process control, or more specifically, to chemical processes [START_REF] Hidalgo | Nonlinear model predictive control of styrene polymerization at unstable operating points[END_REF][START_REF] Eaton | Model-Predictive Control of Chemical Processes[END_REF][START_REF] García | Model predictive control: Theory and practice-A survey[END_REF].

Fortunately, thanks to the improvements on the numerical optimization methods and on the computational capability of the commercial hardware, the applicability of MPC has been extended to systems with fast dynamics. As a matter of comparison, in 1990, Hidalgo and Brosilow used a sampling time of 1 min [START_REF] Hidalgo | Nonlinear model predictive control of styrene polymerization at unstable operating points[END_REF], whereas numerous recent studies apply intricate MPC schemes on robots with sampling time of a few milliseconds, e.g. the trajectory optimization proposed in [START_REF] Neunert | Fast nonlinear Model Predictive Control for unified trajectory optimization and tracking[END_REF] that may be computed within 25 ms. Accordingly, MPC is a mature subject within the field of robotics, e.g. [START_REF] Neunert | Fast nonlinear Model Predictive Control for unified trajectory optimization and tracking[END_REF][START_REF] Ginhoux | Active filtering of physiological motion in robotized surgery using predictive control[END_REF][START_REF] Jacques | High-speed visual servoing of a 6-d.o.f. manipulator using multivariable predictive control[END_REF][START_REF] Vivas | Predictive functional control of a parallel robot[END_REF][START_REF] Heshmati-Alamdari | A self-triggered visual servoing model predictive control scheme for under-actuated underwater robotic vehicles[END_REF][START_REF] Gonzalez | Robust tube-based MPC for constrained mobile robots under slip conditions[END_REF][START_REF] G P Incremona | MPC for Robot Manipulators With Integral Sliding Modes Generation[END_REF][START_REF] Hsieh | Nonlinear model predictive control for wheeled mobile robot in dynamic environment[END_REF][START_REF] Nubert | Safe and Fast Tracking on a Robot Manipulator: Robust MPC and Neural Network Control[END_REF][START_REF] Li | Robust Vision-Based Tube Model Predictive Control of Multiple Mobile Robots for Leader-Follower Formation[END_REF][START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF][START_REF] Ginhoux | Model predictive control to cancel repetitive deformations of organs in robotized laparoscopic surgery[END_REF]. MPC schemes may be categorized according to several classes. The classes that are relevant to the discussions present in this thesis are discussed in the sequel. Firstly, as for the control strategies discussed in Section 1.4, several simplifications can be considered if the system to be controlled is assumed to be linear. In this case, the resulting strategy is classified as linear MPC. Even though linear systems represent a very specific class of system among the broad set of possible applications of MPC, a great part of the most well-known books dedicated to the MPC theory are focused exclusively on linear MPC [START_REF] Marian | Predictive control: with constraints[END_REF][START_REF] Camacho | Model predictive control[END_REF][START_REF] Rawlings | Model predictive control: Theory and design[END_REF]. For this reason, linear MPC is commonly referred to as MPC, shortly. Otherwise, if the control scheme considers a nonlinear system, it is classified as Nonlinear MPC (NMPC) [START_REF] Grüne | Nonlinear model predictive control[END_REF][START_REF] Magni | Nonlinear model predictive control[END_REF].

In the scope of linear MPC, linear systems generally lead to a simple analytical expression for the cost functional optimized in the OCP. Revisiting the example of Figure 1.15, the corresponding transition mapping φ in (1.3) may be written as

s k+1 = A 1 ∆t 0 1 s k + B ∆t 2 /2 ∆t f k . (1.6)
This expression may be recursively applied so that

s 1 = B f 0 s 2 = A B f 0 + B f 1 . . . s n f +1 = n f i=0 A n f -i B f i = B f , (1.7) with B = A n f B A n f -1 B . . . B and f = f 0 . . . f n f T .
As a result, the following explicit expression for the functional in (1.4) may be obtained

x n f +1 -x d 2 + v 2 n f +1 = B f - x d 0 T B f - x d 0 = = f T B T B f -2 x d 0 B f + x 2 d , (1.8)
which is a quadratic function of the vector of forces f . Similarly, the constraints in (1.4) may be written as linear constraints in f . As a result, the OCP (1.4) may be rewritten in the format min

f f T Q f + c T f , s. t. A ineq f b ineq (1.9)
with constant matrices A ineq and Q, and constant vectors b ineq and c, so that the studied OCP is actually a QP problem. Moreover, the fact that the hessian matrix Q = B T B is positive semidefinite demonstrates that the optimization problem is convex.

As in this illustrative example, MPC schemes very often use cost functionals that are quadratic functions of the states and control inputs. For the sake of comparison, among all the previously cited references of MPC applied to robotics [START_REF] Neunert | Fast nonlinear Model Predictive Control for unified trajectory optimization and tracking[END_REF][START_REF] Vivas | Predictive functional control of a parallel robot[END_REF][START_REF] Heshmati-Alamdari | A self-triggered visual servoing model predictive control scheme for under-actuated underwater robotic vehicles[END_REF][START_REF] Gonzalez | Robust tube-based MPC for constrained mobile robots under slip conditions[END_REF][START_REF] G P Incremona | MPC for Robot Manipulators With Integral Sliding Modes Generation[END_REF][START_REF] Hsieh | Nonlinear model predictive control for wheeled mobile robot in dynamic environment[END_REF][START_REF] Nubert | Safe and Fast Tracking on a Robot Manipulator: Robust MPC and Neural Network Control[END_REF][START_REF] Li | Robust Vision-Based Tube Model Predictive Control of Multiple Mobile Robots for Leader-Follower Formation[END_REF][START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF][START_REF] Ginhoux | Model predictive control to cancel repetitive deformations of organs in robotized laparoscopic surgery[END_REF], the work of Worthmann et al. [START_REF] Worthmann | Model Predictive Control of Nonholonomic Mobile Robots Without Stabilizing Constraints and Costs[END_REF] is the only one that does not belong to this class of MPC.

The main characteristic illustrated in the OCP (1.9) is that a cost functional which is a quadratic function of the states and control inputs applied to linear systems may be optimized with standard QP algorithms. Recalling that, as discussed in Section 1.5, numerous software can be used to efficiently solve QP problems (see [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.8]), the numerical solution of such OCPs may be seen as an elementary task. On the other hand, the same kind of cost functional applied to nonlinear systems leads to OCPs demanding nonlinear numerical optimization methods.

Based on this contrast between linear MPC and NMPC, the first control solution proposed in this thesis applies linear MPC (Chapter 3). Numerical simulations and experimental studies proved that the proposed solution yields to satisfying results. Nevertheless, in order to guarantee stability and good performance independent of the level of nonlinearities of the CDPR dynamics, an NMPC scheme is necessary. For this reason, additional studies presented in Chapter 4 address this type of MPC. Both approaches cope with the issue discussed at the very beginning of this section, i.e. the cable tension limits are explicitly handled in the OCP formulated within the MPC scheme. As a result, the tension distribution is integrated within the main feedback control. Moreover, the proposed scheme computes feasible sets of cable tensions even for reference trajectories that escape from the robot workspace.

The majority of studies addressing the stability and robustness of MPC schemes is focused on the so-called tracking MPC (or stabilizing MPC), which penalizes the deviation to a predefined feasible equilibrium or trajectory [START_REF] Rawlings | Model Predictive Control: Theory, Computation, and Design[END_REF]. More precisely, sequences of desired states and desired control inputs should be computed beforehand. The cost functional should penalize the euclidean distance to these desired sequences. Therefore, the cost functional should be null for states and control inputs equal to their desired counterparts, and strictly positive otherwise. The illustrative example of Figure 1.15 can be used to show the typical approach implemented to obtain a standard tracking MPC.

The control scheme considers sequences of desired states {s d,0 , s d,1 , . . . , s d,n f +1 } and forces {f d,0 , . . . , f d,n f } determined beforehand. These sequences represent the desired behavior of the system. The sequence of desired states determines the reference trajectory departing from s d,0 = s 0 and arriving on s d,n f +1 = x d 0 T . The sequence {f d,0 , . . . , f d,n f } is the sequence of forces that generates the desired trajectory. Clearly, both sequences should satisfy the system constraints. For a given discrete instant k and actual state s a , a typical standard tracking MPC could apply the OCP

min {f k ,..., f k+h } k+h i=k s i -s d,i 2 + (f i -f d,i ) 2 (1.10a) s. t. s k = s a (1.10b) s i+1 = φ(s i , f i ) ∀ i ∈ {k, . . . , k + h} (1.10c) |v i | v max , |f i | f max ∀ i ∈ {k, . . . , k + h} (1.10d)
where a positive integer h was used in order to determine the prediction horizon. As described before, the cost functional (1.10a) penalizes deviations with respect to the desired trajectory (with the term s is d,i

2 ) and also with respect to the desired forces (with

(f i -f d,i ) 2
). Moreover, if the actual state is equal to the desired one (s a = s d,k ), the sequence of desired forces within the prediction horizon {f d,k , . . . , f d,k+h } leads to

f i = f d,i and s i = s d,i ∀ i ∈ {k, . . . , k + h}, ⇒ ⇒ k+h i=k s i -s d,i 2 + (f i -f d,i ) 2 = 0 (1.11)
resulting in a null cost functional in (1.10). If any of the elements of {s k , . . . , s k+h } and {f k , f k+h } is not equal to their desired counterparts, the cost functional is strictly positive. This is the desired characteristic described above.

Note, however, that such approach requires the definition beforehand of the desired sequence of control inputs ({f d,0 , . . . , f d,n f } in the example above). The MPC schemes proposed in this thesis perform the tension distribution as an integral part of the main controller. Therefore, the sequence of control inputs (cable tension vectors) is not known beforehand. Accordingly, an alternative approach is proposed in Chapter 4 in order to obtain a standard tracking NMPC with a similar characteristic as in (1.11). The proposed NMPC scheme allows the stability analysis of the corresponding closedloop system. Even though the design and stability analysis of the proposed NMPC is based on the feasibility of the desired trajectory, numerical simulation results show that the controller is able to operate with unfeasible desired trajectories.

The class of tracking NMPC algorithms may be divided between control schemes (i) with and (ii) without terminal conditions. An NMPC scheme with terminal conditions (i) solves OCPs with additional costs and constraints related to the terminal states so that the cost-to-go function can be easily proved to be a Lyapunov function. Accordingly, NMPC schemes with terminal conditions typically lead to simpler stability and performance analysis [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] compared to algorithms without terminal constraints [START_REF] Grüne | NMPC without terminal constraints[END_REF].

Nevertheless, the application of terminal conditions leads to crucial drawbacks in the context of the present thesis. As a matter of fact, "the construction of such stabilizing (terminal) constraints may be challenging and their use may considerably reduce the operating range of the MPC scheme" [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]. The design of the terminal constraints may be particularly troublesome when the disturbances and modeling errors are considered. As illustrated in Figure 1.18, an NMPC of type (i) forces the terminal state in the prediction horizon to be within a terminal constraint set consisting in a neighborhood of the desired trajectory. In the incidence of significant modeling errors and disturbances, feasibility issues may be met. Recalling the considered requirements for the Hephaestus project, one may note that this limitation may lead to safety issues.

Conversely, even though the stability analysis of tracking NMPC schemes without terminal conditions also assumes feasibility of the desired trajectory, controllers of type (ii) are typically less prone to the incidence of unfeasible OCPs. This subject is discussed in detail 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 in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Section 7.4]. Grüne presents an intuitive example in [START_REF] Grüne | NMPC without terminal constraints[END_REF] illustrating the applicability of NMPC with unfeasible desired states.

In addition, terminal constraints may hinder the numerical solution of the OCP [START_REF] Grüne | Nonlinear model predictive control[END_REF] and the implementation of NMPC schemes of type (ii) is typically easier than (i) [START_REF] Altmüller | Performance of NMPC Schemes without Stabilizing Terminal Constraints[END_REF]. Recalling that the numerical implementation of the proposed MPCs should be performed from scratch, these advantages also play an important role. Therefore, the control strategies proposed in this thesis are MPC schemes without terminal conditions. This choice is in accordance with the MPC schemes commonly used in the industry. Indeed, MPC schemes with terminal constraints "are mostly not used in the process in-dustries" [START_REF] Mayne | An apologia for stabilising terminal conditions in model predictive control[END_REF].

Apart from tracking NMPC, it is also interesting to highlight an alternative MPC architecture. Similarly to the control strategies in Chapters 3 and 4 of this thesis, the use of MPC to track unreachable setpoints was also addressed in [START_REF] Rawlings | Unreachable setpoints in model predictive control[END_REF]. A few years latter, the term Economic MPC (EMPC) has been coined referring to MPC algorithms that do not penalize the deviation to a feasible equilibrium (or trajectory) [START_REF] Diehl | A Lyapunov function for economic optimizing model predictive control[END_REF][START_REF] Amrit | Economic optimization using model predictive control with a terminal cost[END_REF]. In other words, MPC schemes that do not match with the definition of standard tracking MPC may be categorized as EMPC. The common use of cost functions evaluating the economic performance of a system inspired the expression "economic". As a matter of fact, as discussed in [START_REF] Grüne | Nonlinear model predictive control[END_REF], the cost function in an EMPC can "model all kinds of quantities, like energy consumption, yield of a substance, income of a firm, etc., which one would like to minimize or maximize". Considering that the ability to track unfeasible desired trajectories is one of the main goal in this manuscript, the control strategies proposed in this thesis may be also classified as EMPC schemes.

As for standard tracking NMPC, the class of EMPC algorithms may be divided between control schemes with [START_REF] Amrit | Economic optimization using model predictive control with a terminal cost[END_REF][START_REF] Angeli | On average performance and stability of economic model predictive control[END_REF][START_REF] Grüne | On non-averaged performance of economic MPC with terminal conditions[END_REF] and without [START_REF] Müller | Economic model predictive control without terminal constraints for optimal periodic behavior[END_REF][START_REF] Grüne | Economic receding horizon control without terminal constraints[END_REF] terminal conditions. Similarly to the strategy described above, an EMPC scheme applying terminal conditions takes advantage of a terminal cost that is a Lyapunov function within the terminal state constraint set. The same disadvantages described before are valid for EMPC schemes with terminal conditions. Stability and performance estimations of this class of NMPC are closely related to the theory of dissipativity of nonlinear systems, which was introduced in [START_REF] Willems | Dissipative dynamical systems part I: General theory[END_REF] and firstly used within the framework of EMPC in [START_REF] Diehl | A Lyapunov function for economic optimizing model predictive control[END_REF]. Whereas the stability conditions of systems with a constant unfeasible reference are well-known [START_REF] Grüne | Economic receding horizon control without terminal constraints[END_REF][START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF], the case for time-varying references is still an active field of research. Most of these studies are focused on periodical reference trajectories [START_REF] Müller | Economic model predictive control without terminal constraints for optimal periodic behavior[END_REF][START_REF] Kohler | On Periodic Dissipativity Notions in Economic Model Predictive Control[END_REF][START_REF] Kohler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]. Some exceptions are the works [START_REF] Dong | Analysis of economic model predictive control with terminal penalty functions on generalized optimal regimes of operation[END_REF][START_REF] Martin | Dissipativity and economic model predictive control for optimal set operation[END_REF], which address the case for general optimal operation regimes.

Independently of the aforementioned classes, MPC schemes involve the numerical solution of OCPs in real-time. On the other hand, the vast majority of control schemes discussed in Section 1.4 defines a precomputed policy. More precisely, studies [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]- [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF] define an explicit expression to the computation of the desired wrenches. The comparison between the complexity involved in these two approaches reveals an important drawback related to MPC schemes. As discussed in the end of Section 1.2, the control schemes developed in this thesis were implemented in a software commonly used in industrial environments. Since such software is not compatible with common programming libraries, the control schemes proposed here were implemented from scratch. This is a minor aspect in the scientific perspective, but crucial from a practical point of view.

Previous studies addressing the implementation of MPC in the control of CDPRs are very few and recent. Katliar et al. proposed a nonlinear MPC for a motion simulator in [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF]. Its performance was investigated through numerical simulations. To the best of our knowledge, no experimental results were published so far. The feasibility of the MPC implementation in real time was addressed taking into account solvers such as HPMPC [START_REF] Frison | High-performance small-scale solvers for linear model predictive control[END_REF] and qpOASES [START_REF] Hans | qpOASES: a parametric active-set algorithm for quadratic programming[END_REF]. However, these solvers are not compatible with common industrial real-time environments. As a motion simulator control, the focus in [START_REF] Katliar | Nonlinear Model Predictive Control of a Cable-Robot-Based Motion Simulator[END_REF] is a set of desired accelerations and velocities. The present thesis is focused on trajectory tracking, which prioritizes positioning accuracy.

The recent work [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF] presents simulations and experimental results of a linear MPC scheme implemented in a hybrid cable-driven robot. However, the studied cabledriven system consists of a 2-DoF platform with planar movements. Only some of the cables are controlled in tension. The remaining cables are controlled in position. The actual tensions of these cables are not considered in the real-time controller and the proposed MPC strategy is used to control the vibration of the system. This thesis present results on trajectory tracking of CDPRs with six DoF, eight cables, and each cable tension is controlled by the MPC scheme. Additionally, the focus of the work [START_REF] Qi | Decoupled modeling and model predictive control of a hybrid cable-driven robot (HCDR)[END_REF] is vibrations attenuation. One may note that both schemes proposed in [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] cannot be classified as position tracking control.

Contributions of the Thesis

Based on the state of the art discussed in Sections 1.4, 1.5, 1.6 and 1.7, MPC proved to be an appropriate strategy to the context introduced in Section 1.2. The state-ofthe-art control strategies do not explicitly consider the cable tension limits in the computation of the feedback correction. Therefore, feasibility issues may be faced when the robot operates close to the system constraints. These risks present particular relevance in this thesis. Being one of the few control strategies able to explicitly handle the system constraints, MPC corresponds with the specific requirements in the Hephaestus project. Therefore, this thesis proposes different position tracking control solutions using MPC. As a matter of fact, the integration of the TDA within the main feedback control computation is a relevant achievement for the control of CDPRs, in general. Details on the design, numerical simulations, implementation, and real-time experimental results related to two MPC schemes are presented.

Chapter 3 introduces a linear MPC scheme, highlighting its main advantage: The control scheme is able to perform the tracking of unfeasible trajectories satisfying the constraints of feasible cable tensions. In other words, the performed trajectory is as close as possible to the desired one without violating the interval of feasible cable tensions. To the best of our knowledge, this controller led to the first experimental results of a CDPR control scheme with this characteristic. Additionally, characteristics related to disturbance rejection capabilities are also discussed.

Acknowledging the nonlinearities intrinsically inherent to the dynamics of a CDPR, Chapter 4 addresses an NMPC. Whereas this control scheme is inspired by the strategy of Chapter 3, it seeks guarantees of stability and performance based on the theoretical studies of tracking NMPC. Consequently, a more intricate scheme is proposed and its stability analysis is performed. Similarly to the linear MPC of Chapter 3, terminal conditions are not used. The output of both MPC schemes is a set of desired feasible cable tensions for each sampling time.

Even though the results presented in Chapters 3 and 4 represent the main contributions of this thesis, the corresponding implementation in experimental set-ups would not have been possible without the development of some other technical solutions. Some of these technical solutions are discussed in Chapter 2.

Since the output of the proposed predictive controllers is a set of desired cable tensions, a strategy should be implemented to convert desired cable tensions to a variable that is meaningful to the motor drivers, e.g. motor torques, positions or velocities. As discussed in Section 1.6, the CDPRs studied in this thesis have force sensors in order to measure the cable tensions. Therefore, even though the majority of published studies address this problem with open-loop cable tension control schemes, measurements obtained with the force sensors may be used in order to guarantee that the actual tensions are as close as possible to the desired ones. Nevertheless, the CDPRs considered in this thesis have servo motors coupled to gear trains with a large reduction ratio. Real-time experiments described in Chapter 2 demonstrate that the influence of the friction introduced by these gears is non-negligible and scarcely repeatable. Inasmuch as every feedback state of the art cable tension control method discussed in Section 1.6 is highly dependent on a precise identification of the system dynamics, a non-model based scheme was developed. Chapter 2 presents the design of this controller and the corresponding experimental results.

In addition, in accordance with the priorities enumerated in Section 1.2, an efficient estimation of the pose of the platform should be used in order to improve the positioning accuracy of the CDPR. As can be seen in Table 1.1, a relatively large cable diameter (12 mm) is necessary due to the the large payload presented in the Hephaestus project. As a result, a large pulley diameter is necessary in order to satisfy the minimal curvature specific for this cable, which increases the influence of the pulley kinematics discussed in Section 1.3. Therefore, a forward kinematics algorithm considering the pulley kinematics is proposed in Chapter 2. The main contribution here consists of explicit expressions of the corresponding differential kinematics. Such expressions facilitate the implementation of an efficient algorithm using QR factorization.

One may note that implementation of the methods discussed in Chapters 2, 3 and 4 demands advanced numerical algebra, e.g. QP, QR factorization, null space computation, among others. Since industrial software and hardware should be used (see Section 1.2), the real-time operation of the programmable logic controller (PLC) is most often not compatible with existing libraries of advanced mathematical operations. Consequently, the necessary operations were programmed from scratch using C++ language compatible with the used industrial software (TwinCAT -Beckhoff). The Appendix A is devoted to the these implementations.

Additionally, Appendix B introduces a TDA able to minimize the infinity norm of the vector of cable tensions. In contrast to most of the algorithms compatible with m > n + 1, the proposed TDA consists in a closed-form expression. The optimality of the proposed algorithm is proved and experimental results are presented. The discontinuity inherent of TDAs minimizing the infinity norm is exemplified.

Kinematic Model

The main goal of this section is to develop a method with which the most probable platform pose and velocity could be estimated for given measured motor positions and velocities. On the one hand, precise estimations are necessary in order to obtain an accurate position tracking control algorithm. On the other hand, since this estimation should be done in real-time, the computation time is a major concern.

Even though Pott showed in [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF] that relevant errors may be obtained assuming fixed proximal anchor points, most of the previous works discussed in Section 1.3 do not address the influence of the pulley kinematics. An exception is the work of Schmidt [START_REF] Lorenz | Modeling Techniques and Reliable Real-Time Implementation of Kinematics for Cable-Driven Parallel Robots using Polymer Fiber Cables[END_REF], in which the FK problem considering the pulley kinematics is solved but the equations are not derived explicitly.

Accordingly, the present section explicitly derives the differential kinematics for CDPRs considering the pulley kinematics and, based on that, proposes a real-time capable Forward Kinematics (FK) algorithm. The algorithm applies iteratively QR decomposition to solve a linearized version of the least squares problem representing the FK. Differential kinematics consists here of a closed-form expression of the jacobian matrix of CDPRs considering the pulley kinematics. This jacobian matrix is used to construct the linearization of the FK problem. Experimental and numerical results address the convergence capabilities of the proposed algorithm.

Following the typical approach for the kinematic modeling of parallel robots, the kinematic constraints are firstly formulated by means of the inverse kinematics.

Inverse kinematic model

A 6-DoF CDPR is considered but the proposed algorithm can be applied to CDPRs with fewer DoF. The platform is driven by m cables. The main variables and parameters of the kinematic model are introduced in Figure 2.1. The position of each cable attachment point B i on the platform is given by the vector b i , with i = 1, ..., m. Each point A i is a fixed point which lies on the circumference of pulley i. Each pulley i is able to rotate around the line passing through A i tangent to the pulley circumference. The constant position of A i is given by vector a i in the fixed reference frame.

The platform pose is given by the vector x = p T ψ T T . The position of the platform is p = p x p y p z T . Its orientation is represented by the vector ψ = α β γ T . The elements of this vector are Euler angles so that b

i = p + R z (γ) R y (β) R x (α) b pi ,
where b pi is the vector of the platform attachment point i written in the coordinate system O p , attached to the platform, and R x , R y and R z are the rotation matrices around x, y and z axis, respectively. The vector of cable lengths is l = l i ... l m T , where l i is the length of cable i and the cable elongations are neglected. The inverse kinematics model gives the motor positions for a given pose of the platform, i.e. compute l for a given x. The remainder of this section recalls the inverse kinematics considering the influence of the pulley geometry. Similar models were notably presented in [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF][START_REF] Bruckmann | Wire Robots Part I: Kinematics, Analysis & Design[END_REF].

The pulley geometry is illustrated in Figures 1.7 and 2.2. One drawing pulley i among the m CDPR pulleys is considered and subscript i is dropped in the remainder of this section. Consider the definition of coordinate systems E 0 = e 0 1 e 0 2 e 0 3 and E 1 = e 1 1 e 1 2 e 1 3 . These coordinate systems are defined with respect to the inertial frame O b . The pulley can rotate about the fixed axis e 0 3 (with angle ρ). E 0 is any fixed coordinate system with e 0 3 aligned with the rotation axis of the pulley. E 1 rotates with the pulley and e 1 3 = e 0 3 . E 1 is thus pose dependent and its orientation is computed in the inverse kinematics.

Define v = b ia i , written in the CDPR fixed reference frame. Vector a i is constant for a given CDPR and b i depends on the platform pose. Vector v can be written in the coordinate system E 0 as v 0 = E T 0 v. The angle ρ shown in Figure 2.2 is obtained as atan2(v 0 2 , v 0 1 ), with v 0 1 and v 0 2 the first and second components of v 0 . Therefore, the coordinate system E 1 is obtained as

E 1 = E 0 R z (ρ)
, where R z is the rotation matrix around the z axis.

The vector v 1 written in the coordinate system

E 1 is v 1 = E T 1 v = v 1 1 0 v 1 3 T
. The component v 1 2 = 0 because the coordinate system is aligned to the pulley orientation so that vector v lies in the plan defined by the vectors e 1 1 and e 1 3 . This can be seen in Figure 2.2.

Analyzing the geometry presented in Figure 2.2, the following system of trigonometric equations is obtained

v 1r = v 1 1 v 1 3 = l f sin φ cos φ + r p 1 -cos φ sin φ , (2.1)
where r p is the primitive radius of the pulley and l f is the length depicted in the Figure 2.2. Eliminating the variable l f , (2.1) leads to

v 1 3 sin φ + (r p -v 1 1 ) cos φ -r p = 0. (2.2)
Trigonometric functions sin φ and cos φ may be written as

sin φ = 2 tan(φ/2) 1 + tan 2 (φ/2) , cos φ = 1 -tan 2 (φ/2) 1 + tan 2 (φ/2) , (2.3) so that (2.2) leads to tan 2 (φ/2) (v 1 1 -2 r p ) + tan(φ/2) (2 v 1 3 ) -v 1 1 = 0. (2.4)
Equation (2.4) can be solved as a quadratic equation in tan(φ/2). Then, using the inverse tangent function, two solutions are obtained in the interval [-π, π], namely,

φ 1 = 2 atan   v 1 1 2 -2v 1 1 r p + v 1 3 2 -v 1 3 v 1 1 -2r p   , ( 2.5
)

φ 2 = 2 atan   -v 1 1 2 -2v 1 1 r p + v 1 3 2 -v 1 3 v 1 1 -2r p   . (2.6)
If v 1 1 > 2r p (typical situation), φ 1 > 0 and φ 2 < 0. In this case the angle respecting the geometry of the problem is φ 1 . Otherwise, φ 2 should be taken.

The length l f may be easily obtained from (2.1) once φ is known. The cable length l i is given by

l i = r p φ + l f + l ci , (2.7)
where l ci is a constant additional length. This constant additional length may, for instance, take into account the total distance between fixed pulleys.

Additionally, the unitary vector u i shown in Figure 2.1 can be computed as

u i = E 1 R y (φ) 0 0 -1 T = -E 1 sin φ 0 cos φ T . (2.8)
This vector is used in the computation of the wrench matrix, as detailed in Section 3.1.1.

Differential kinematics

The jacobian matrix J is used in the solution of the FK problem since it relates δl to δx as δl = J δx, (

where δ(•) denotes the differential variation of (•). This section presents the analytical expression of the matrix J obtained from the equations in Section 2.1.1. As in Section 2.1.1, subscript i is dropped.

The vector v = b i -a i is computed as v = p + R(ψ) b pi -a i , (2.10) with R(ψ) = R z (γ) R y (β) R(α).
Differentiating this expression with respect to the platform pose, the following equation is obtained

δv = I 3 ∂ ∂ψ (R(ψ) b pi ) ∂v ∂x δx, ( 2.11)
where

∂ ∂ψ (R(ψ) b pi ) = =    b piy σ 1 + b piz σ 2 b piz cα c β cγ -b pix cγ s β + b piy c β cγ sα b piz σ 3 -b piy σ 4 -b pix c β sγ -b piy σ 3 -b piz σ 4 b piz cα c β sγ -b pix s β sγ + b piy c β sα sγ b piz σ 3 -b piy σ 2 + b pix c β cγ b piy cα c β -b piz c β sα -b pix c β -b piz cα s β -b piy sα s β 0    , with σ 1 = (s α s γ + c α c γ s β ), σ 2 = (c α s γ -c γ s α s β ), (2.12
)

σ 3 = (c γ s α -c α s β s γ ), σ 4 = (c α c γ + s α s β s γ ), (2.13)
and c θ = cos θ and s θ = sin θ.

Defining v 0r = v 0 1 v 0 2 T
and E 0r = e 0 1 e 0 2 , one may write

v 0r = E T 0r v = l 0 cos ρ sin ρ , ( 2
v 1r = R T zr (ρ) E T 0 v = l f sin φ cos φ + r p 1 -cos φ sin φ . (2.16)
Differentiating this equation, the following relation is obtained

∂v 1r ∂x D Rz (ρ) E T 0 v ∂ρ ∂x + R T zr (ρ) E T 0 ∂v ∂x δx = = r p sin φ cos φ + l f cos φ -sin φ δφ + sin φ cos φ δl f , (2.17)
with

D Rz (ρ) = -sin ρ cos ρ 0 0 0 0 . ( 2.18)
The differential δφ can be isolated in (2.17) leading to the derivative of

φ δφ = 1 l f cos φ -sin φ ∂v 1r ∂x ∂φ ∂x δx. (2.19)
Moreover, this expression can be substituted for δφ in (2.17) and δl f is thereby obtained as

δl f = sin φ cos φ ∂v 1r ∂x -r p ∂φ ∂x ∂l f ∂x δx. (2.20)
Differentiation of (2.7) leads to

δl i = r p ∂φ ∂x + ∂l f ∂x ∂l i ∂x δx. (2.21)
Finally, the computations detailed above can be applied to each kinematic chain i = 1, ..., m, and the derivatives are concatenated to form the expression of the Jacobian matrix

δl =     δl 1 . . . δl m     =     ∂l 1 ∂x . . . ∂lm ∂x     J δx. (2.22)

Forward kinematic model

The FK model consists in computing the platform pose for a given set of winch motor positions. In the case of redundant CDPRs with m > n, n being the number of degrees of freedom of the mobile platform, the kinematic model is overdetermined. In general, due to modeling and measurement inaccuracies, it is impossible to find a platform pose that satisfies exactly the kinematic constraints set by the m kinematic chains. Therefore, the FK may be formulated as the problem of minimizing the error between the measured cable lengths and the cable lengths computed by the IK with (2.7).

Let l m be the set of measured cable lengths and l(x) the cable lengths obtained with the IK model for a given pose x. The error to be minimized is defined as e(x) = l(x) -l m 2 , where (•) 2 denotes the 2-norm of (•). Therefore, for a given l m , the FK algorithm should return the solution x * = arg min x l(x) -l m 2 .

In general, the function e(x) = l(x) -l m 2 possesses several local minima. The proposed FK algorithm takes an initial guess x g and find a x * that locally minimizes the function e(x). The proposed algorithm assumes that x g is sufficiently close to the current platform pose and consists of an iterative scheme. For an iteration k with platform pose x k , the next iteration takes

x k+1 = x k + ∆x k . Considering the approximation ∆ l ≈ J(x k ) ∆x, a reasonable choice for ∆x k is ∆x k = arg min ∆x J(x k ) ∆x -l m -l(x k ) 2 .
(2.23)

The minimization problem (2.23) is a Linear Least Squares problem. The solution of (2.23) might be computed solving

J(x k ) T J(x k ) ∆x = J(x k ) T l m -l(x k ) , ( 2.24)
for ∆x. However, SVD and QR factorizations may solve this linear least-square problem more efficiently [START_REF] Gene | Matrix computations[END_REF][START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF]. The SVD may be preferred because of its greater diagnostic capability in pathological cases. Nevertheless, QR factorization presents faster computing time. The latter being critical in real-time applications, the QR factorization is preferred. Therefore, the Jacobian matrix J is decomposed using the QR factorization such that

J = Q R,
with Q orthogonal and R upper triangular. The solution of (2.23) is then obtained from the back substitution of the system

R ∆x = Q T l(x k ) -l m .
Once ∆x k is obtained, the next pose x k+1 = x k + ∆x k is computed. This procedure is repeated until ∆x 2 < , for the desired tolerance.

Note that the proposed algorithm does not use numerical damping. This was not necessary for the studied scenarios. Numerical damping is applied in commonly used Levenberg-Marquadt methods in order to compensate the incidence of ill-conditioned configurations. If a given iteration k reaches a pose in which J(x k ) T J(x k ) is close to singular, one may solve J(x k ) T J(x k ) + λ I ∆x = J(x k ) T l m -l(x k ) instead of (2.24), for a positive λ ∈ R. This procedure improves numerical stability.

Experimental results

The proposed FK algorithm has been implemented in the HRPCable prototype shown in Figure 1.1. As introduced in Section 1.2, this prototype has a 6-DoF platform fully constrained by 8 cables. The control control software has been developed in C++ and runs in an industrial PC Beckhoff C6920 equipped with 2.4GHz i7 core processor. The platform pose control loop runs at 125 Hz. An inner feedback loop running at 2 kHz control the cable tensions according to the control scheme proposed in Section 2.3. Cable tensions were measured by means of force sensors (load pins) placed in the drawing pulleys. Figure 2.3 depicts the performed trajectory. The threshold was defined as = 1 × 10 -6 .

In addition to the components shown in Figure 1.1, a Metris K600 camera system is used to measure the pose (6-DoF) of the platform in real time with a precision of 70 µm. The measurements obtained with the K600 camera system were used as reference to compare with the estimations obtained with the proposed FK algorithm. These results are shown in Figure 2.4. This figure also shows the number of iterations necessary to the FK algorithm to converge for each sampling time during the trajectory. Figure 2.4 indicates the presence of significant errors obtained with the proposed FK algorithm. Such errors are probably a consequence of the assumption that the cables are inelastic. Accordingly, an identification of the elastic behavior of the cables combined with the cable tension measurements may lead to an improved kinematic model. This is, indeed, a subject of possible future works. For a given pair of pose and length vectors x * and l * consistent with the constraints presented in section 2.1.1, rigorously evaluating the capability (denoted here as convergence capability) of an FK iterative algorithm to quickly and reliably find x * for given l * is not a trivial task and is out of the scope of this thesis. The FK problem of a Stewart-Gough platform (topologically very similar to CDPRs) may have up to 40 solutions [START_REF] Wampler | Forward displacement analysis of general six-in-parallel SPS (Stewart) platform manipulators using soma coordinates[END_REF] and redundant configurations do not necessarily present a reduced number of possible solutions compared to the non-redundant ones [START_REF] Husty | A redundant spatial Stewart-Gough platform with a maximal forward kinematics solution set[END_REF]. Therefore, it is necessary to prove that the algorithm is able to converge to x * , and not just to one of the (potentially) many x that are consistent with l * . Nevertheless, it is still interesting to address this capability as follows.

The FK algorithm has been applied to IK solutions corresponding to more than 30 × 10 3 poses equally spaced across the workspace keeping the initial guess constant equal to x g = 0 0 .8 0 0 0 T for all the poses. The workspace of the prototype is a 4 × 2 × 1.4 m 3 (x, y, z) cuboid, considering -10 • γ 10 • , α = β = 0. The algorithm obtained a pose with errors smaller than the proposed tolerance within 7 iterations for all these poses. Clearly, this test was not performed in real-time. The cable length vectors l corresponding to the IK solutions were artificially fed to the FK algorithm.

A Torque-Based Cable Tension Control

The vast majority of CDPR position tracking control strategies for CDPRs define, at some point, a set of desired cable tensions. The present and next sections are dedicated to the control strategies responsible for generating these cable tensions by controlling variables that can be directly managed in the motor drives, namely, motor currents (or torques) and velocities.

As discussed in Section 1.6, the cable tension control in CDPRs is often performed by means of an open-loop scheme based on a model relating the motor torques to the cable tensions. Nevertheless, the precision obtained with such models may be deteriorated for winches presenting significant friction. Clearly, the implementation of a closed-loop strategy measuring the actual cable tensions may reduce this error.

Moreover, fully-constrained CDPRs typically have cable tension sensors for safety reasons. Hence, it is appropriate to consider the inexpensive information of measured cable tensions within these strategies. Therefore, Sections 2.2 and 2.3 focus on closedloop strategies.

Based on the intuitive relation between cable tensions and motor torques, a torquebased scheme is firstly proposed. Several previous works presented satisfying results with torque-based cable tension control schemes operating in open-loop [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF][START_REF] Chellal | Model identification and visionbased H-infinity position control of 6-DoF cable-driven parallel robots[END_REF]. Nevertheless, the conducted experimental tests on the HRPCable prototype indicate that a torque-based closed-loop scheme would have its implementation hindered or substantially complicated by the presence of significant friction in the gear trains. This section aims to provide theoretical and experimental data allowing to analyze the applicability of a torque-based scheme.

Before addressing the control schemes themselves, the experimental set-up used to evaluate the proposed schemes is introduced.

Problem formulation and experimental setup

Sections 2.2 and 2.3 address cable tension control with different approaches. Accordingly, the results presented in both sections are based on the same experimental set-up. This set-up is illustrated in Figure 2.5. Two kinematic chains of the HRPCable prototype (cf. Figure 1.1) are used. The mobile platform is not used. The corresponding cable ends are attached together.

The resulting system may be seen as a 1-DoF mechanism driven by two cables. The position of the attachment point is given by the variable x. The torques applied by the motors in winches 1 and 2 are denoted by γ 1 and γ 2 , respectively. Similarly, motor positions are denoted by q 1 and q 2 . Consider given desired positions x d (t) and desired cable tensions τ d (t) in function of time. The present section proposes a control scheme meant to generate measured cable tensions τ 1 as close as possible to the desired ones controlling motor torques γ 1 . The motor drive of winch 2 is controlled in position mode. The position q 2 (t) is controlled considering the desired position x d (t) and inelastic cables. More precisely, the motor position q 2 is controlled as

q 2 (t) = q 2 (t 0 ) + k qx (x d (t) -x d (t 0 )), (2.25)
with t 0 the initial time and a constant k qx linearly relating cable length variations to the angular displacements of the motors.

Note that, whereas the final goal is to generate the desired cable tension applied on the attachment point P , the considered goal is to track τ d (t) based on τ 1 (t). Therefore, the influence of distributed mass of the cable and the friction on the pulleys are neglected. These are the main causes of divergence between measured tension in the load-pin 1 and the force applied in point P . Even though, this subject is addressed in Section 2.3.4.

As shown in Table 1.1, the diameter of the cables is 4 mm and the diameter of the pulleys is 80 mm. The winches are composed of Beckhoff AM8061 motors coupled to AG2210 gear trains. These gear trains have a reduction ratio equal to 25. Load pins Sensy 5300-1T with integrated amplifier are positioned on the pulley axes. It is worth noting that a sufficiently small influence of noise could be obtained with this configuration of force sensing. Therefore, filtering was not necessary.

Dynamic Model

As stated in the previous section, the cable distributed mass is neglected so that the cable tension is considered constant all over the cable length. Under this assumption, the model of the system depicted in Figure 2.5 is reduced to the model of the winches.

A detailed model of the actuator may take into account each body of the actuator. Take as example the illustrative gear train depicted in Figure 2.6. This gear train consists of three bodies. These bodies are identified by letters a, b and c. Body b presents two gears, identified by numbers 1 and 2. Each body presents a corresponding inertia (I a , I b and I c ). The angular positions are denoted by q a , q b and q c . The external torques γ a and γ c are applied on axes a and c, respectively. The number of teeth of the gears are N a , N b1 , N b2 and N c . By neglecting backlash, the rigid-body dynamic model of the system can be written as

I c qc = (γ a -I a qa ) N b1 N a -I b qb N c N b2 -γ c -γ F ( qa , qb , qc , γ a , γ c ). (2.26)
The torque γ F is a function of ( qa , qb , qc , γ a , γ c ) that suitably models the friction of the system. Denoting

c 1 = Nc N b2 N b1
Na and c 2 = Nc N b2 , all velocities and accelerations may be written in function of the derivatives of

q c          qa = c 1 qc qa = c 1 qc qb = c 2 qc qb = c 2 qc . (2.27)
Therefore, a new function may be written as

γ f ( qc , c 1 τ 1 , τ 3 ) = γ F (c 1 qc , c 2 qc , qc , τ 1 , τ 3 )
. Furthermore, the model (2.26) may be simplified as

(I a c 2 1 + I b c 2 2 + I c ) qc = -γ c + c 1 γ a -γ f ( qc , c 1 γ a , γ c ).
(2.28)

Denoting I = (I a c 2 1 + I b c 2 2 + I c ), γ = c 1 γ a , γ out = γ c
and q = q c , an equivalent single body system is obtained with

I q = γ -γ out -γ f ( q, γ in , γ out ).
(2.29)

Even if model (2.29) seems rather simple, a precise model usually involves an elaborate function γ f [START_REF] Shiqi | Nonlinear dynamic model and governing equations of low speed and high load planetary gear train with respect to friction[END_REF][START_REF] Yang | Hertzian Damping, Tooth Friction and Bending Elasticity in Gear Impact Dynamics[END_REF]. Highly nonlinear phenomena such as stick-slip may present non-negligible influence in the studied system.

For the system illustrated in Figure 2.5, γ is the motor torque multiplied by the gear ratio. The output torque γ out is the torque applied by the cables γ t . Note that γ t = τ r w , where τ is the cable tension and r w is the radius of the winch drum.

Let y be the state vector consisting of the angular position q and velocity q of the winch drum so that

y(t) = y 1 (t) y 2 (t) = q(t) q(t) . (2.30)
In order to take into account the stick-slip phenomenon, the following state-space model is proposed

ẏ(t) =            0 0
, if |γ| γ sd and y 2 = 0;

y 2 1 I (γ -γ kf (y 2 ) -τ r w ) , othewise (2.31)
where γ sd is the static dry friction torque. The scalar γ kf is the kinetic friction torque, which is function of the angular velocity y 2 = q. Different definitions of the function γ kf (y 2 ) were proposed in the literature. It is worth noting that the model typically considered in the control of CDPRs (such as in [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][START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF][START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF]) is described by (2.31) with γ sd = 0. The model considered in [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF] is often used in this field and is given by

γ kf 1 (y 2 ) = f kd tanh(µ y 2 ) + f v1 y 2 , ( 2.32)
where f kd is the kinetic dry friction coefficient and f v1 is a viscous friction coefficient. The function tanh(µ y 2 ) is a continuous approximation of the function sign(y 2 ) and the constant µ is used to tune this function. However, Section 2.2.3 shows that a better consistency with the experimental data was obtained with the following definition of the kinetic friction torque

γ kf (y 2 ) = f kd tanh(µ y 2 ) + f v1 y 2 + f v2 sign(y 2 ) |y 2 |, (2.33)
with a second viscous friction coefficient f v2 . The model presented in (2.33) is used in the remainder of this section. The obtained results presented in Section 2.2.3 show that this model is able to fit experimental and numerical data thanks to the addition of the term f v2 sign(y 2 ) |y 2 |.

Next, before presenting the proposed control scheme and the corresponding results, the parameters of (2.31) and (2.33) need to be identified. The next section addresses this subject.

Parameters Identification

Since the dynamic parameters (inertia and friction constants) of the actuator are independent of the cable tensions, all experimental data described in this section are obtained with the winch running freely (no cable is attached to the actuator). The first identified parameter was the static dry friction torque γ sd . Consider as initial state γ = q = 0. Increasing gradually γ, velocity q will still equal to zero until γ = γ sd . The value of γ sd may be estimated using this procedure. The presented torques are given by the motor driver based on the measured motor current. The experiment depicted in Figure 2.7 was repeated several times in order to obtain a good estimation of γ sd . According to model (2.31), the torque values highlighted in Figure 2.7 are equal to the static dry friction torque. The conducted experiments show that the behavior of the actuator regarding the friction is significantly different considering the direction of the torque. The analysis of the friction torque has thus been done separately for each direction ( q positive and negative). The obtained results are summarized in Table 2.1.

q > 0 q < 0 Mean Value 3.575 3.207 Standard Deviation 0.149 0.257 The remaining friction coefficients have been estimated repeating the experiment depicted in Figure 2.8. Piecewise constant motor torque setpoints are sent to the drives. A constant value is set until the transient behavior ends. For each value of desired motor torque, the related angular velocity has been recorded after the transient interval. Thereby, the relationship between the recorded velocities and torques are independent of the inertia. This procedure have been repeated several times, leading to the chart depicted in Figure 2.9. The coefficient µ was set empirically equal to 10 3 . A gradient based optimization (interior-point method) yielded to the coefficients of model (2.31) minimizing the difference between the predicted torque and the experimental one. The obtained results are also shown in Figure 2.9. Similarly to the estimation of γ sd , different values of friction coefficients were obtained for q > 0 and q < 0. This figure also shows the relevance of the nonlinear term f v2 sign(y 2 ) |y 2 | in (2.33).

Finally, the inertia of the actuator has been identified by applying sinusoidal torques with different frequencies and amplitudes. The value of I = 507.5 kg.cm 2 minimizes the difference between the measured profile of q and the one given by the model, taking measured torques as inputs.

Proposed control scheme

The torque applied on the winch is closely related to the cable tensions. For this reason, a unique control loop where the motor torque is the control input may seem reasonable. Anyhow, the lowest level control input applied on the motor is the torque/current. Accordingly, the present section proposes a control scheme that is based on the studies presented in [START_REF] Chellal | Model identification and visionbased H-infinity position control of 6-DoF cable-driven parallel robots[END_REF][START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF].

The proposed scheme is depicted in the block diagram of Figure 2.10 and it may be expressed as

γ = k P (τ d -τ ) + k I ta t 0 (τ d -τ ) dt -k D q + γ kf ( q) + ẍ I r w , ( 2.34)
with γ kf given by (2.33) and gains k P , k I , k D . Instants t 0 and t a denote the initial and actual time, respectively.

The friction and inertia effects are compensated by the terms γ kf ( q) and ẍ I/r w , respectively. The proportional and integral terms reduce the error between measured tension τ and desired tension τ d . Stability is enhanced by the term -k D q. Note that the addition of a term proportional to the derivative of cable tensions τ is not recommended due to the presence of noise.

Numerical simulation results

This section presents simulations considering the control scheme (2.34), the actuator model (2.31) and the parameters identified in Section 2.2.3. In addition to the actuator model, the cable elasticity needs to be modeled in order to obtain numerical results consistent with the experimental ones.

In contrast to the experimental set-up, the position x(t) may be simulated without the presence of pulleys. Accordingly, the simulated system may be simplified as presented in Figure 2.11. It consists of an actuator (motor, gears and drum) connected to a cable. The other end of the cable is attached to a moving point P . This point moves in one direction and its displacement is given by x. A force sensor measures the cable tensions. The real stretched length of the portion of the cable that is not in contact with the winch drum is l 1 . The unstretched length of this portion of the cable is equal to l 2 . The angular position of the drum is q. The model (2.31) is applied including the linear elastic cable model

t = max EA (l 1 -l 2 ) l 2 , 0 , (2.35)
where the term EA is the specific elastic constant of the cable. It is equal to the Young Modulus of the material multiplied by the area of the cross section of the cable. This model is based on the one proposed in [START_REF] Mohammad | Dynamic modeling and control of parallel robots with elastic cables: Singular perturbation approach[END_REF]. Consider E = 110 GPa [181] the Young Modulus of stainless steel, and the fill rate equal to 0.7 (70% of the area defined by the nominal diameter of the cable is actually filled with steel). The resulting value of EA is equal to 9.68 × 10 5 N.

Simulations take as initial conditions l 1 = l 2 = l 0 and x = q = 0. This way, l 2 = q r w + l 0 and l 1 = l 0 -x. Moreover, simulations considers a cycle time of 1 ms. Figure 2.12 presents the obtained results. The controller is not able to continuously track the desired cable tensions due to the so-called stick-slip phenomenon. This phenomenon results in a discontinuous relation between the input torque and the cable tension. This discontinuity cannot be compensated with a torque-based continuous control policy. More precisely, although γ kf may be compensated in (2.34), the discontinuous system (2.31) introduces the stick-slip torque represented by γ sd that is not considered in the control law. Clearly, a discontinuous torque-based control scheme compensating the stick-slip phenomenon in (2.31) is not viable.

The aforementioned issue is present for velocities q close to zero. This is the condition in which the system is susceptible to the discontinuity represented in (2.31). For sufficiently large velocities, the proposed scheme (2.34) may lead to satisfying results.

Issues related to a torque-based control approach

Section 2.2.5 demonstrated that the application of the proposed torque-based control approach may lead to problems resulting from the presence of friction in the winches. Indeed, preliminary experimental results confirmed this issue. More specifically, the main concern related to the control of the cable tensions is the presence of static dry friction in the actuators. In some cases, this may be a minor problem. For instance, for the CDPR CoGiRo (Figure 1.2 (b)), motor torques of each winch are transmitted to the corresponding drum by a system of pulleys and belt with a reduction ratio equals to 3. Whereas each motor of HRPCable is connected to a 2 stage gear train Beckhoff AM8061 with a reduction ratio equals to 25. This leads to increased friction torques. Figure 2.13 illustrates this problem. In this experiment, one end of the cable is attached to a fixed point and the values of the measured cable tensions are obtained for different motor torques. One can see that a continuous torque variation results in abrupt changes of cable tensions.

As a matter of fact, every experiment performed with variations of the torque-based closed-loop scheme (2.34) led to technical issues. The resulting closed-loop system is not able to track cable tensions when the velocity q is small.

In addition, the estimations presented in Section 2.2.3 are time-consuming and lead to parameters identified with a bad level of precision. Considering that such estimations should be repeated for each one of the winches, a bad compromise regarding the complexity of implementation and the level of precision is obtained.

Within this context, a non-model based scheme able to minimize the influence of friction may be an appropriate alternative to the proposed torque-based scheme. This forms the rationale of the method presented in the next section.

A Velocity-Based Cable Tension Control

Section 2.2.6 discussed the main limitations of torque-based cable tension control schemes when gear trains with large reduction ratio are used. Since this is the case of both HRPCable and Hephaestus prototypes, an alternative approach is proposed in this section.

As the motor torques, one may note that there is an intuitive relation between cable tensions and motor positions. For a given x, variations on q lead to variations on τ according to the cable elasticity. Figure 2.14 illustrates this idea. The presented test is very similar to the one used to obtain the results of Figure 2.13: one cable end is attached to a fixed point and the cable tension is measured while the winch is controlled. The difference between the tests performed in Figure 2.13 and Figure 2.14 is that the motor is controlled in position mode instead of torque mode.

The motor position is changed continuously while the cable tensions are measured. As a matter of fact, the results presented in Figure 2.14 are rather obvious considering the elastic behavior of the cable. Nevertheless, the comparison between Figures 2.13 and 2.14 demonstrates that a more promising scenario is obtained with the later approach.

In general, a trivial position control is able to perform a sufficiently precise position tracking of q independently of the external torques applied on the winches. This means that satisfying positioning precision of q may be obtained independently of the cable tension τ . That said, the motor position q may be controlled in order to generate the desired cable tension τ d . Noting that the relation between the cable tension τ and the position q is continuous and well behaved, one may conclude that better results than the ones presented in Section 2.2 may be obtained controlling the motor position (instead of torque). This rationale would lead to a cascade control scheme in which the outer loop takes as control input the motor position and the inner loop controls the motor torque. Thereby, the influence of the friction within the gear trains would be managed by the position (inner) control loop, which presents good robustness capabilities against model uncertainties. This is indeed the approach proposed in [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF][START_REF] Kraus | Force control of cable-driven parallel robots[END_REF]. Nevertheless, some secondary complications come up with this method. Namely, a precise identification of the system dynamic may be necessary and the suppression of steady-state errors demands the application of an anti-windup feedback strategy (refer to [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF] for further details). For this reason, the present section proposes a scheme in which the motor velocity is taken as control input. This method allies the aforementioned advantage without leading to secondary issues.

Control scheme

Figure 2.15 depicts the proposed velocity-based cable tension control. An outer feedback loop takes the cable tension τ as measured state and delivers the desired winch velocity qd as intermediary control input. This outer loop applies a linear PI correction and a feedforward term to compensate the velocity corresponding to ẋ. An inner loop defines the motor torque taking the desired velocity as reference. Similarly, it applies a linear PI term and applies a feedforward term compensating the inertia and friction. In the considered experimental setup, the servo drives are responsible for this inner loop. Depending on the features available on the motor drive, an alternative feedforward term may be used. The control setpoint delivered to the motor drive is the desired motor velocity. The servo drives may present lower cycle time (between 62.5 µs and 250 µs for the servo drives used in this chapter) than the PLC cycle time. The reduced sampling time used in the torque control loop may lead to better performances.

For a given actual instant t a , the mathematical expressions for qd and γ in Figure 2.15 are

qd = k V P (τ d -τ ) + k V I ta t 0 (τ d -τ ) dt + ẋ r w , (2.36a) γ = k T P ( qd -q) + k T I ta t 0 ( qd -q) dt + γ kf ( q) + ẍ I r w . (2.36b)
The values of ẋ and ẍ are inputs of the block diagram in Figure 2.15. These are the time derivatives of x, which represents the desired cable length. Therefore, ẋ and ẍ can be determined using the inverse kinematics proposed in Section 2.1.

Since the inner loop is performed by the servo drives in the experimental set-up, the control (2.36b) is used uniquely for the simulations. Note that the control implemented in the experimental set-up is the outer loop (2.36a) which consists of a PI feedback control with feedforward ẋ/r w . Therefore, the implementation and tuning are rather simple.

Numerical simulation results

As discussed at the beginning of this section, the influence of friction should be reduced by controlling the motor velocity instead of the motor torques directly. Accordingly, in contrast to the data presented in Section 2.2.5, the numerical simulations using the proposed velocity-based control yielded to satisfying results. As depicted in Figure 2.16, this strategy is able to keep a small error. Maximum error for this simulation is equal to 0.539 N with desired cable tensions varying with an amplitude of 400 N.

The spikes in the cable tension error and torque are consequences of the stick-andslip phenomenon. The discontinuity in (2.31) is faced when the time derivative of the desired cable tensions is close to zero. Nevertheless, in contrast to the scheme proposed in Section 2.2, the position control is sufficiently robust to reduce the influence of the friction.

Real-time experimental results

In accordance with the description given in Section 2.2.1, two winches of the CDPR HRPCable were used as shown in Figure 2.5. Motor 1 is controlled with the velocity control described in section 2.3.1 and cable tension obtained with load pin 1 is taken as feedback signal. Motor 2 is controlled in position mode according to (2.25). This scenario represents the conditions in which the cable tensions should be controlled while the attachment point position is variable. This application is equivalent to the cable tension control necessary to actuate a CDPR. As for the simulations presented in Section 2.3.2, the controller gains used in the present section were tuned by trial and error. N is 0.4 s. The position x was kept constant. Figure 2.18 shows the results for τ d varying with low frequency. This figure shows that the obtained error is negligible so that desired and measured cable tensions are virtually equivalent. Moreover, the small error obtained with a slow variation of τ d shows that the proposed strategy successfully avoids the stick-slip problem. control scheme without the feedforward term would lead to large errors when variable x is used. The small errors obtained when the feedforward term is applied indicates that the compensation of ẋ is satisfying.

Finally, Figure 2.20 presents the results with sinusoidal x and τ d . A small delay may be observed. Nevertheless, the obtained results with the proposed control strategy are considered suitable. Accordingly, the velocity-based cable tension control (2.36) is experimentally validated and may be used in the position tracking discussed in the next chapters.

Further results

Some additional results were obtained during the experiments described above. Figure 2.21 shows several repetitions of the experiment depicted in Figure 2. [START_REF] Martin | Hybrid control with multicontact interactions for 6 DOF haptic foot platform on a cable-driven locomotion interface[END_REF]. One can see that the relation between the cable tension and the displacement of the winch (q r w ) is repeatable. The plot displaying the elongation vs. cable tension indicates that the relation between stress and strain for increasing cable tensions can be reasonably approximated with an affine function. Nevertheless, the values obtained for decreasing tensions shows that significant hysteresis is present.

As shown in Figure 2.5, there is a second load pin measuring the cable tensions in pulley 2. Figure 2.22 presents the values obtained with both sensors. One may see that the cable tensions measured by sensor 2 shift when the direction of movement changes (velocity equals to zero). This shift is probably caused by the friction in pulleys 3 and 4 shown in Figure 2.5. Since τ 1 is the value used in the control loop, it is continuous and close to the desired cable tension.

In order to clarify the influence of the friction in the pulleys, consider a simplified model τ2 = τ 1 -2 τ f sign( ẋ) + τ k , (2.37) where τ f is the the dry friction force due to pulleys 3 and 4, the friction forces in these pulleys being considered equal. The scalar τ k is a constant that depends on the initialization of the system. Parameters τ f and τ k were identified using data from Figure 2.22. The fitted data is shown in Figure 2.23. The pulley dry friction force obtained is τ f = 3.9 N.

This correction may be easily integrated in the cable tension control. However, considering the small influence of τ f = 3.9 N in the interval of feasible cable tensions, this correction is not included in the remainder of this thesis.

Linear Model Predictive Control

A Linear Model Predictive Control (LMPC) for the position tracking of CDPRs is proposed and its performance is evaluated through numerical simulations and real-time experiments. Section 3.1 introduces the control algorithm and the dynamic model on which the control scheme is based. In order to simplify the formulation of the corresponding optimal control problem and its numerical solution, a linear discrete-time approximation is proposed based on the original continuous nonlinear dynamic model. Simulation results presented in Section 3.2 indicate that such approximation is pertinent for the Hephaestus prototype. Finally, Section 3.3 details the experimental validation of the proposed linear MPC. These results prove that the MPC scheme copes with the main concerns discussed in Chapter 1, namely, the cable tension limits are explicitly handled as an integral part of the main controller.

Control Design

The main motivations to use an MPC strategy for the position tracking control of CDPRs were discussed in Section 1.7. In brief, the feedback correction of the stateof-the-art schemes discussed in Section 1.4 are subjected to feasibility issues since the cable tension limits are not considered in the main computation of the feedback loop. A known advantage related to MPC schemes lies in its ability to explicitly handle the system constraints. Considering that, the present section proposes a linear MPC scheme able to handle the cable tension limits as an integral part of the main feedback control. This MPC is based on the dynamic model discussed in Section 3.1.1, while the control strategy is introduced in Section 3.1.2.

Dynamic model

Before addressing the model itself, some notations are introduced. Sequences of vectors are denoted with bold lowercase letters followed by (•) (as, for instance, s(•)). The set of sequences of vectors with dimension n s and infinite elements is denoted as S ns . Similarly, S n s2 n s1 denotes the set of sequences containing n s1 vectors with dimension n s2 . Accordingly, the k th vector of a sequence s ∈ S n s2 n s1 is denoted as s(k) = s k ∈ R n s2 . The set of strictly positive integers is referred to as N. The set of non-negative integers is

N 0 = N ∪ {0}. Moreover, for a, b ∈ R, N a,b = {i ∈ N 0 | a i b}.
This subsection presents the dynamic model of a CDPR consisting of an n-DoF mobile platform driven by m cables, where n ≤ m. The platform pose in time is given by the vectorial function x : R → R n , so that x(t) = p(t) T ψ(t) T T is the pose vector containing both the platform position p(t) and orientation ψ(t) at instant t. The dependence on time is dropped leading to the short notation x whenever there is no risk of confusion. The same practice is used for other vectorial and matricial functions in the remainder of this thesis.

Typically, ψ consists of Euler angles. The common dot notation is used to refer to the time derivatives, so that dx/dt = ẋ and d 2 x/dt 2 = ẍ. In general, the angular velocity ω is not equal to the time derivative orientation vector ψ. For this reason, the following relations need to be introduced

ṗT ω T T = S(x) ẋ pT ωT T = S(x) ẍ + Ṡ(x, ẋ) ẋ, (3.1)
with S(x) a square matrix responsible for transforming ψ into ω. On the other hand, Ṡ(x, ẋ) is the component-wise time derivative of S(x).

Considering that the platform pose and cable forces are known, the magnitude τ i and the direction u i of the force applied by the cable i ∈ N 1,m on the platform are also known. Vector u i is shown in Figure 2.1 and can be computed with (2.8). Accordingly, in this thesis, the vectors u i are uniquely computed in function of the platform pose. In contrast, if cable sagging is considered, the direction u i also varies in function of the cable tension itself. This variation is neglected here. Therefore, we denote these vectors as vectorial functions u i (x).

The vectorial function τ c : R → R m returns the vector of cable tensions in function of time τ c (t) = τ 1 (t) . . . τ m (t) T . Each cable force is applied on the corresponding point B i of the mobile platform (see Figure 2.1). Summing up the contributions of forces and moments of every cable force, the resultant wrench f (x, τ ) applied on the platform is obtained with the linear mapping

f (x, τ c ) = W(x) τ c , ( 3.2)
with the wrench matrix W(x) representing this linear mapping. In a general 6-DoF CDPR, as for vectors u i , the vectors b i in Figure 2.1 may be written as functions of the platform pose with

b i (x) = R(ψ) b pi , ( 3.3)
where R and b pi were defined in Section 2.1.1. The wrench matrix W is then by

W(x) = u 1 (x) . . . u m (x) b 1 (x)u 1 (x) . . . b m (x)u m (x) , ( 3.4)
with v1 denoting a skew-symmetric matrix so that v1

v 2 = v 1 × v 2 is the cross product between vectors v 1 , v 2 ∈ R 3 .
Knowing the wrench applied by the cables on the platform, its motion is modeled using Newton-Euler formalism, which leads to the following expression:

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

M (x) = m p I -m p ĉ(x) m p ĉ(x) H(x) and (3.6) C (x, ẋ) ẋ = m p ωx ωx c(x) ωx H(x) ω x . (3.7)
The scalar m p is the platform mass and I is the identity matrix with dimensions suitably chosen. Denoting as c p the vector going from the platform geometric center to its center of mass, expressed in the coordinate system O p attached to the platform, the corresponding vector in the global coordinate system is computed as

c(x) = R(ψ)c p = c x c y c z T .
(3.8)

The angular velocity ω x is obtained from x according to (3.1). The matrix H is defined as H(x) = R(ψ) I G R(ψ) T + m p ĉ(x) ĉ(x) T where I G is the platform inertia matrix in respect to the coordinate system O p . The vector of gravitational forces is The dynamic model (3.10) is considered as the nominal system in the remainder of this thesis. This is a continuous-time nonlinear system. Nevertheless, as discussed in Section 1.7, the MPC strategy proposed in this chapter needs a linear discrete-time system. To this end, the sequences of vectors x(•), ẋ(•), ẍ(•) ∈ S n and y(•) ∈ S 2n are defined as

g(x) = m p g 0 0 -1 -c y c x 0 T ,
x(k) = x k = x(t 0 + k ∆t) ẋ(k) = ẋk = ẋ(t 0 + k ∆t) ẍ(k) = ẍk = ẍ(t 0 + k ∆t) y(k) = y k = x T k ẋT k T (3.11)
for all k ∈ N 0 , an initial time t 0 and a sampling period ∆t. Without loss of generality, since the continuous system (3.10) is time-invariant (as defined in [62, Chapter 1]) the initial time is considered to be t 0 = 0.

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

(3.12)

Using the Euler method of numerical integration, the state space representation of (3.10) may be approximated by

y k+1 ≈ A I ∆t I 0 I y k + B(y k ) 0 ∆t M(x y k ) -1 W(x y k ) τ k + + 0 ∆t M(x y k ) -1 g(x y k ) -C(x y k , ẋy k ) ẋy k v(y k ) . (3.13)
with x y , ẋy ∈ R n denoting the pose and velocity vectors extracted from a vector y ∈ R 2n . Additionally, 0 denotes a matrix or vector with dimensions suitably chosen and all its elements equal to zero.

In accordance with the notation introduced in the equation above, the considered approximated discrete model is

y k+1 ≈ A y k + B(y k ) τ k + v(y k ).
(3.14)

Control strategy

As introduced in Section 1.7, a key element of an MPC algorithm is the formulation of an OCP based on a bounded time span. This time span is determined by h p ∈ N, called prediction horizon. For given instant k, state y k and finite sequence of cable tensions τ ∈ S m hp , an estimation of the sequence of future states should be defined based on the discrete-model (3.14).

State Predictions

In order to simplify the computation of the estimated future states, the variations of the matrix B and the vector v in (3.14) within the prediction horizon is assumed to be negligible. Therefore, the following approximations are considered

B(y k+i ) ≈ B(y k ) v(y k+i ) ≈ v(y k ) ∀ i ∈ N 1,hp . (3.15)
In other words, for a given time instant k, these matrices are considered to be constant within the prediction horizon and are denoted shortly as B = B(y k ) and v = v(y k ).

Moreover, as Maciejowski discusses in [START_REF] Marian | Predictive control: with constraints[END_REF], considering every possible combination of cable tensions in the prediction horizon may lead to a bad compromise between performance improvement and computational burden. Therefore, it is interesting to consider every combination within a reduced horizon h c < h p , denoted as control horizon, and constraint the last i th components with h c < i h p , as follows

τ k+i = τ k+hc ∀ i ∈ N hc,hp . (3.16)
Using the proposed notations and assumptions, the future states along the prediction horizon are estimated applying recursively (3.14) with the identities (3.15) yielding to

γ(y k ,u)               ỹk+1 ỹk+2 ỹk+3 . . . ỹk+hc+1 . . . ỹk+hp               = D               A A 2 A 3 . . . A hc+1 . . . A hp               y k + E(y k )                 B 0 . . . 0 AB B 0 . . . A 2 B AB . . . 0 . . . . . . . . . A hc B A hc-1 B . . . B . . . . . . . . . . . . A hp-1 B ... A hp-hc+1 B hp-hc i=0 A i B                 u       τ k τ k+1 . . . τ k+hc       + ν(y k )               v A v A 2 v . . . A hc v . . . A hp-1 v              
or, with a compact notation, γ(y, u) = Dy + E(y) u + ν(y).

(3.17)

Note that, since B and v are functions of y, E and ν are also functions of y.

Optimal Control Problem Formulation

Our goal is to define a functional cost of the MPC scheme including (i) the predicted errors, (ii) the magnitude of the cable tensions and (iii) the variation of the cable tensions. The minimization of the cable tensions variation (iii) is important mainly because large variations of cable tensions may be not physically feasible. This issue is also relevant during the initialization of the robot. If the initial cable tensions are not close to an optimal cable tension distribution, the MPC will automatically find a smooth transition between these two sets of cable tensions.

Consider sequences x d (•), ẋd (•) ∈ S n nt describing a desired trajectory, with t f = n t ∆t the trajectory final time. The vector of desired states within a prediction horizon h p for an instant k is γ d (k) = x d (k + 1) T ẋd (k + 1) T . . . x d (k + h p ) T ẋd (k + h p ) T . The weighted expressions of (i), (ii) and (iii) are then respectively given by

J γ (y, u, k) = γ d (k) -γ(y, u) 2 Kγ (3.18) J u (u) = u 2 Ku (3.19) J D (u, τ c (k ∆t)) = k ∆u hc-1 i=0 ∆τ i+k 2 (3.20)
where K γ and K u are positive definite diagonal weighting matrices and k ∆u is a positive scalar. In addition, w 2 K = w T K w denotes the weighted norm of w ∈ R nw with a weighting symmetric positive semi-definite matrix K ∈ R nw×nw . The vector of the variation of cable tensions can be expressed as

∆τ i = τ u (i) -τ u (i -1), if i > 1 τ u (1) -τ c (k ∆t), if i = 1 , ( 3.21)
where τ u (i) ∈ R m is obtained by extracting the i th vector of cable tensions of u. Alternatively, (3.21) may be written as

∆u = ∆τ T 1 . . . ∆τ T hc T = Q u -z, (3.22)
with

Q =          I m 0 . . . -I m I m 0 . . . 0 -I m I m . . . . . . . . . . . . . . . 0 0 . . . . . . -I m I m          , z =       τ c (k ∆t) 0 . . . 0      
.

Note that z depends on the measured cable tensions. Alternatively, the desired cable tension obtained at t-∆t may be used. The dimensions of Q and z are m h c ×m h c and m h c × 1, respectively.

Defining K D = k ∆u I, J D can be rewritten as Using (3.17), the overall cost function J(y, u) is defined as the sum of J γ , J u and J D neglecting the terms independent of u according to

J D (u, τ c (k ∆t)) = u T Q T K D Q u -2 z T K D Q u + z T K D z. ( 3
J(y, u, k, τ c (k ∆t)) = u T Hc (E T K γ E + K u + Q T K D Q) u+ + 2 (D y + F -γ d (k)) T K γ E -z T K D Q d T u.
(3.24)

The proposed MPC scheme consists in computing the optimal sequence of cable tensions by solving in real-time the following QP problem u * = arg min The desired set of cable tensions is used as setpoint for the cable tension control introduced in Section 2.3.1. The cable tension control is responsible for generating actual cable tensions τ c (t) for (k -1)∆t t < k ∆t as close as possible to τ * k .

Figure 3.1 depicts the block diagram of the implemented control scheme, showing that the kinematic model and cable tension control introduced in Chapter 2 are necessary complements to the MPC strategy. The Cartesian velocity ẋMPC predicted by the MPC is used in the cable tension control. This vector is obtained with

y k+1 = A y k + B τ k + v.
Alternatively to the optimization problem (3.25), the constraints may be changed using

u * = arg min U 1 2 u T H c u + d T u s.t. u min u u max |∆u| ∆u max (3.27)
with |∆u| representing the absolute values of ∆u given by (3.22) and a scalar ∆u max limiting the maximal cable tension variations. This way, the variation of cable tensions is explicitly constrained. In this case, the scalar k ∆u may be set to 0.

Numerical Simulations

The present section compares through numerical simulations the performance obtained with the MPC scheme proposed in Section 3.1 with those obtained with stateof-the-art strategies. To this end, the control schemes with which the proposed MPC is compared should be revisited. The results presented in this section were published in [START_REF] Cavalcanti | Model Predictive Control Applied to Large-Dimension Cable-Driven Parallel Robots[END_REF].

Background on the state-of-the-art controllers

The performance obtained with the proposed MPC is compared to those led by strategies commonly used for the control of CDPRs. Namely, MPC is compared to: (i) a linear PID+ controller [START_REF] Paden | Globally asymptotically stable 'PD+' controller for robot manipulators[END_REF] and (ii) a Sliding Mode Control (SMC) [START_REF] Slotine | Sliding controller design for non-linear systems[END_REF]. A brief description of these two control methods is presented in the sequel.

For a given reference trajectory in time t 0 t t f , the desired poses, velocities and accelerations are denoted by x d (t), ẋd (t) and ẍd (t), respectively. At a given instant t, the error in the Cartesian space is expressed as e x (t) = x d (t) -x(t). Similarly, the error in joint space is denoted by e j (t) = l d (t) -l(t), where l d (t) is the vector of desired cable lengths obtained from x d (t) with the inverse kinematics and l(t) is the estimated current cable length vector computed according to motor positions.

PID + computed torque

This control strategy, hereinafter referred to as PID+, applies the following wrench f = M(x) ẍd + C(x, ẋ) ẋdg(x) + W(x) K p e j + K i t t i e j (τ ) dτ + K d ėj , (3.28) where K p , K i and K d are diagonal matrices containing the linear feedback PID gains.

Sliding Mode Control

The SMC strategy defines a sliding surface s = e x + C e ėx , with C e a positive definite diagonal matrix. The wrench control law to be applied on the platform is

f = M ẍd + C d ( ẋd -ẋ) + K sat(s) + Q s + C(x, ẋ) ẋd -g(x), (3.29)
where K, Q and C d are diagonal feedback gain matrices. The function sat(s) is a continuous approximation of the sign function. Each component of this vector-valued function is calculated as follows

sat(s i ) =      1, if s i > ∆ s i ∆ , if |s i | ∆ -1, if s i < ∆ . (3.30)
The resulting function presents the same output than the sign(s) function except for the interval -∆ s ∆ in which a linear interpolation eliminates the discontinuity.

Redundancy Resolution

The two above control strategies define the wrench f to be applied on the platform. The final control output is the vector of cable tensions τ (or motor torques). For fullyconstrained CDPRs, m > n and some TDA should be implemented, as discussed in Section 1.5. In the simulations described in this section, the 2-norm of the vector of cable tensions is minimized according to In some cases, the wrench demanded by the controller may be not feasible. More precisely, in the space of cable tensions, the intersection of the subspaces defined by the constraints (3.32) and (3.33) is empty. Another strategy should then be defined and the following optimization problem is proposed

min τ W τ -f P s.t. τ min τ τ max (3.34)
where the subscript P indicates that the 2-norm is calculated with a weighting positive definite diagonal matrix P, which is necessary since f has components with inconsistent units (forces and moments).

Simulations results

This section presents simulation results comparing the performances obtained with the control strategies presented in subsections 3.1 and 3.2.1. The OCP formulation constraining cable tension variations according to (3.27) is used. The context of these simulations is the project Hephaestus described in Section 1.2. A large-dimension CDPR is intended to automatize several tasks in the construction and maintenance of building facades. The main task of the CDPR is the installation of curtain wall modules. The robot workspace is a rectangular region in front of the building facade. Thereby, the CDPR mobile platform can pick up curtain wall modules on the ground and position them where needed on the building facade. Since the CDPR will operate in an outdoor environment, it will be subjected to external disturbances. One of the main concerns is the incidence of wind gusts. For this reason, the simulations presented in this section are focused on external disturbance rejection performances of the studied control strategies. An impulsive disturbance is applied and the response of the CDPR is analyzed. Note that the simulated trajectory is relatively short. However, the simulation of a longer trajectory would not affect the results which highlight the disturbance rejection capabilities of each control strategy.

The initial and final positions are depicted in Figure 3.2. The path between these two positions is a straight line segment. The trajectory is the fastest possible respecting upper bounds on linear velocities, accelerations and jerks. These bounds are 0.3 m/s, 0.3 m/s 2 and 1.0 m/s 3 , respectively. The resulting trajectory has continuous derivatives up to the acceleration level. The desired orientation of the platform is constant along the trajectory. An impulsive disturbance f d is applied at the instant t = 2 s. This impulsive wrench is f d = 55 55 550 0 0 0 T (N and Nm) and is applied at the reference point of the platform.

The CDPR configuration (cable drawing points, cable-platform attachments, and cable arrangement) can also be seen in Figure 3.2. Moreover, the parameters of the CDPR dynamic model are the following: τ min = 100 N, τ max = 14 kN, m p = 1000 kg, c p = [0 5 0] T m, I G = diag([400 100 400]) kg.m 2 .

In the following, the results obtained with the three studied motion control strategies are presented and discussed. The control parameters used in the simulations are the following:

K p = 71400 I, K i = 71400 I, K d = 71400 I, C e = 2 × 10 -3 I, C d = 36 I, Q = 40 I, K = 0.2 I, h p = 20, ∆ t = 6×10 -4 , ∆ r = 5×10 -1 , K γ,p = 6×10 9 , K γ,v = 1×10 -2 , K u =
2 × 10 -5 I and k ∆u = 0. These gains were tuned with trial-and-error method.

The components of K γ related to the pose errors are set equal to the scalar K γ,p , whereas K γ,v is used for velocity errors. Value of ∆ t is used in (3.30) as ∆ for translational inputs of sat(s), and ∆ r is used for rotational inputs. trajectory.

The histograms of Figure 3.4 present a performance comparison on different aspects. Let e t be the 2-norm of the translational error. Histogram (a) quantifies the maximum value e t over the trajectory. Histogram (b) compares the RMS value of e t along the trajectory. Taking these two performance measures, MPC leads to the smallest error.

Regarding cable tension values, as shown in Figure 3.4-(c), MPC demands the maximal allowed value τ max = 14 kN. This is also visible in Figure 3.5, which depicts the cable tensions near the instant of application of the impulsive disturbance. Indeed, as discussed earlier, the main advantage of MPC is that the controller takes into account the constraints of the system and optimizes the control actions in order to reduce the tracking errors. Here, the maximum allowed cable tension is an active constraint just after the impulsive disturbance is applied. In the case of SMC, the maximum tension value is 13.7 kN along the trajectory, which indicates that this controller response is close to the largest admissible value τ max = 14 kN. If higher gains were used, the wrench f would then be unfeasible. In order to not exceed the limit τ max , the strategy described in Eq. (3.34) would be necessary and the resulting wrench would not be equal to f . However, the MPC strategy does not lead to this risk. The PID+ controller has |τ max | = 9.1 kN. This relatively low tension is a consequence of the use of small gains K p . Indeed, small gains were used because larger gains lead to high frequency oscillations of cable tensions. For instance, increasing the gains of less than 1% with this strategy leads to RMS( τ ) equal to 86 kN.

The time derivative of the cable tensions is a measure of the degree of aggressivity of the control action. Large values of this variable may excite high frequency dynamics which are difficult to control. Figure 3.4-(d) presents the maximum derivative of cable tensions over the trajectory considering all cables. Mathematically, the values in Figure 3.4-(d) are equal to max t max i | τi (t)| . The smallest maximum cable tension derivative is obtained for the PID+, which is an advantage of this strategy. SMC is the most aggressive controller considering this performance measure. The proposed MPC strategy constraints this variable (according to (3.27)). Therefore, any value can be imposed independently of the rest of the controller parameters. In simulations, the value used is 800 kN/s and shows that this value is reached. Figure 3.4-(e) presents the values of the maximum RMS value of τi considering all cables i ∈ N 1,m . More precisely, Figure 3.4 (e) presents values of max i RMS( τi ) . The MPC strategy presents the largest RMS of cable tension derivatives. Note that variables depending on the system states beyond the prediction horizon cannot be taken as constraints in the MPC optimization. The prediction horizon covers only a small part of the whole trajectory over which the RMS values are calculated.

The consumed energy over the trajectory is roughly the same for all the control schemes, cf. Figure3.4-(f). The consumed energy is calculated as

t f t i | l(t) T τ (t)|dt.

Real-time Experiments

As frequently concluded for MPC applications, the numerical solution of OCPs in real-time is a crucial challenge in the implementation on an experimental prototype. Moreover, since industrial software is used, this implementation could not take advantage of commonly used numerical libraries. Finally, the experiments unveil issues that may be neglected in simulations, e.g. the cable tension control and forward kinematics discussed in Chapter 2.

Experimental platform and implementation aspects

The control scheme introduced in Section 3.1 is experimented on the prototype HRPCable, a 6-DoF CDPR driven by 8 cables, installed in LIRMM facilities. Figure 1.1 shows this experimental setup and the corresponding CAD model. Details on the components of this CDPR are given in Table 1.1.

The MPC control scheme has been implemented in an industrial PC and developed in TwinCAT environment, using C++ language. As discussed before, TwinCAT is a software from Beckhoff commonly used in the industry, which turns PC-based systems into real-time controllers using Microsoft Windows kernel. In order to guarantee real-time performances, TwinCAT is not compatible with standard C++ libraries, including basic libraries such as math.h. This issue generally mitigates the application of MPC with constraints in industrial robotics. The control methods described in the previous sections were thus programmed from scratch. Some notes on the numerical implementation are presented in Appendix A. The OCP formulation (3.25) is used.

In the experiments reported below, the control parameters are the following: These gains were tuned manually.

h p = 6, h c = 3, K u = 3.3 × 10 -3 I, k ∆u = 0.3. Matrix K γ is diagonal

Experimental results

The present section experimentally evaluates the efficiency of the proposed MPC scheme with respect to the main goals considered in this thesis. Summarizing some of the key concerns presented in Section 1.2, the proposed method should combine satisfying positioning precision, disturbance rejection capabilities and should safely operate close to system constraints. Section 1.7 shows that this latter characteristic presents particular relevance in the context of this thesis. Indeed, state-of-the-art control strategies are not able to explicitly handle system constraints within the computation of the main feedback correction and, therefore, may be vulnerable to feasibility issues during the robot operation. The proposed MPC strategy is not prone to this issue since it handles the tension limits within the OCP formulation (3.25). In order to show the effectiveness of the proposed controller, Section 3.3.2.1 presents the results for a typical pick-and-place task performed with two different scenarios: τ max = 400 N and τ max = 260 N. The same cable tension limits were assumed for all the cables. The desired trajectory escapes from the robot workspace when the reduced maximum tension τ max = 260 N is applied. The proposed MPC is able to comply with this reduced maximum cable tension while minimizing the Cartesian error between the actual and the desired trajectories. This is an important result related to the safety of the operation of CDPRs. Moreover, Section 3.3.2.2 evaluates the capabilities of the proposed MPC related to robustness against uncertainties.

Pick-and-Place task

The first scenario considered is a typical pick-and-place task. Figure 3.6 shows the sequence of the six desired platform poses defining the pick-and-place path. The trajectory between each subsequent pair of desired poses is defined with a 5 th degree polynomial function. In the following, the results obtained with two different scenarios of cable tension constraints are compared. possible to desired poses complying with these cable tension limits. The translation errors reach values greater than 250 mm and orientation errors are close to 8 • . Figure 3.8 depicts overlaid views comparing the scenarios with nominal and reduced maximum tension. This comparison is done for two different instants. The view in Figure 3.8-(a) was taken when the desired pose was feasible for both scenarios. In contrast, Figure 3.8-(b) depicts an instant in which the desired was feasible for the nominal case and unfeasible for the reduced maximum cable tension. The shift between desired and actual poses is evident in this case.

Nominal constraints:

Strictly complying with the cable tension limits while minimizing the Cartesian errors in following the desired trajectory is the main advantage of the proposed controller. The control strategies discussed in Section 1.4 combined with a tension distribution schemes described in Section 1.5 would fail to fulfill this objective since the main controller would demand an unfeasible wrench. Figure 3.9 depicts the results obtained using the control scheme proposed in [START_REF] Lamaury | Control of a Large Redundantly Actuated Cable-Suspended Parallel Robot[END_REF]. The robot operation is suddenly stopped when the desired trajectory reaches unfeasible poses, exemplifying the risks taken with such control strategy when operating close to system constraints.

Robustness against payload uncertainties

In many applications of CDPRs, the platform should pick weights from one position and release them in another position. In order to validate the applicability of the proposed control strategy in such tasks, it is important to evaluate its robustness against uncertainties on the lifted mass. To this end, the experiment illustrated in Figure 3.10 is proposed. The results are presented in Figure 3.11. An additional mass of 11.5 kg was used. This mass represents 50% of the platform mass m p . The controller is able to keep the tracking errors smaller than 6 mm and 0.14 • , despite this significant uncertainty.

Conclusions

The experimental results presented in this section show that the proposed control strategy is able to perform a trajectory keeping reduced errors while complying with cable tension limits. The real-time control was implemented in an industrial software and hardware environment to enable applicability of the proposed scheme in industry.

The proposed MPC scheme is able to address the cable tension limits explicitly, integrating the redundancy resolution within the main controller. As a result, cable tension limits are not violated even for reduced maximum cable tensions. Indeed, when a desired pose cannot be reached with some given tension limits, the proposed MPC is In order to extend the results obtained in this chapter to system with significant nonlinearities, an NMPC scheme should be developed. The next chapter addresses this subject.

Nonlinear Model Predictive Control

The linear MPC presented in the previous chapter is based on a linear approximation of the CDPR dynamic model. As a result, it may present deteriorated performance in presence of significant nonlinearities. Accordingly, a nonlinear MPC is proposed in this chapter. After a preliminary discussion presented in Section 4.1, the overall NMPC algorithm and its optimal control problem are formulated in the Section 4.2. Since this strategy is based on the original nonlinear dynamic system, the stability of the corresponding closed-loop may be analyzed. This section also highlights additional aspects of the linear MPC proposed in the previous chapter that differ from the standard form of tracking MPC schemes. Notably, in order to attain this standard form, an alternative method aiming to the minimization of the cable tensions is proposed. Section 4.3 discusses issues related to the numerical discretization of the original continuous-time system and proposes numerical solutions for the optimization problems proposed in Section 4.2. The NMPC strategy is validated by numerical simulation results presented in Section 4.4.

Nonlinear Model Predictive Control

The present section outlines the main objectives and the rationale of the NMPC proposed in this thesis. Before addressing the problem itself, some definitions and notations are introduced in the following.

Definitions and notations

In addition to the previously used notations, some additional sets should be introduced

R + := {r ∈ R | r > 0}; R 0 + := {r ∈ R | r 0}; R nv + := {v ∈ R nv | v i > 0 ∀ i ∈ N 1,nv }; R a,b := {r ∈ R | a r b}; U := {τ ∈ R m | τ min τ τ max } (set of feasible cable tensions); U nu := {τ (•) ∈ S m nu | τ min τ i τ max ∀ i ∈ N 1,nu }; Y := y ∈ R 2n y = x T ẋT T
with x and ẋ admissible ;

X := x ∈ R n ∃ ẋ ∈ R n , x T ẋT T ∈ Y .
The set Y should restrict the set of poses and velocities within which the control is able to operate. This definition should use the tools related to the computation of the workspace of CDPRs, e.g. [START_REF] Gouttefarde | Analysis of the wrench-closure workspace of planar parallel cable-driven mechanisms[END_REF][START_REF] Pierre | On the workspace of suspended cable-driven parallel robots[END_REF][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 Sections 4.2 and 4.4.

The set X is the projection of Y in the space of platform poses. This set represents the admissible poses (poses x ∈ R n for which there exists a velocity vector ẋ ∈ R n such that x T ẋT T ∈ Y).

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

K := {α ∈ R + 0 → R + 0 | α continuous, strictly increasing and α(0) = 0}; K ∞ := {α ∈ K | α unbounded }; L := {δ ∈ R + 0 → R + 0 | δ continuous, strictly decreasing and lim t→∞ δ(t) = 0}; KL := {β ∈ R + 0 × R + 0 → R + 0 | β(r, •) ∈ L and β(•, t) ∈ K}. (4.1)
As in Chapter 3, any vectorial function representing the cable tensions in time τ c : R → R m is considered piece-wise constant. Accordingly, a sequence of cable tensions

τ (•) ∈ S m is related to τ c by τ (k) = τ k = τ c (t), for t ∈ t (k -1)∆t t < k ∆t and k ∈ N, (4.2)
for a given sampling period ∆t. The initial time is considered t 0 = 0. Note that, as in Chapter 3, a sequence of vector is denoted with bold letters followed by (•), such that the k th vector of a sequence s(•) ∈ S ns is denoted by s(k) = s k .

Moreover, the euclidean distance between two vectors

v 1 , v 2 ∈ R nv , for any n v ∈ N, is denoted as v 1 v 2 = v 2 v 1 = (v 1 -v 2 ) T (v 1 -v 2 ).
The reader is invited to refer to the Glossary at the end of this manuscript, which is a summary of the used notations.

Problem formulation

Since the cost functional (3.24) considers a linear approximation (3.17) of the continuous model (3.10), the satisfying results presented in Section 3.2 and 3.3 are obtained only for dynamic systems with limited nonlinearities. Therefore, a rigorous stability analysis of a general CDPR nonlinear dynamic model and its closed-loop is hindered. In this context, the goal of the present chapter is to obtain an NMPC scheme based on the original nonlinear model leading to satisfying performance and stability independently of the influence of nonlinearities.

Thanks to the simplicity resulting from the linear model, several issues in the design and implementation of a general NMPC could be addressed in Chapter 3 without a clear separation between themselves. Due to the increased complexity involved in the NMPC controller proposed in the present chapter, it is important to decouple these issues and address them separately. Accordingly, the design and implementation of an NMPC scheme may be divided into four aspects, namely:

(i) Numerical discretization: Definition of a transition-mapping representing a discrete-time approximation of the continuous-time state space model;

(ii) Optimal Control Problem (OCP) formulation: Definition of the cost functional to be minimized, as well as the control input and state constraints;

(iii) Formulation of the OCP as a Nonlinear Programming (NLP) problem: In general, the OCP obtained in (ii) is not formulated as a standard NLP. Typically, the transition mapping defined in (i) is not explicitly considered in the OCP. Therefore, (i) and (ii) should be integrated, leading to a standard NLP problem;

(iv) Numerical solution of the NLP: The output of the controller is obtained by solving numerically the NLP formulated in (iii). Interior-point and Sequential Quadratic Programming (SQP) methods are numerical algorithms commonly used to this end.

Numerical discretization is necessary since most of NMPC schemes assume discretetime systems and the CDPR dynamic model (3.10) is continuous in time. More in detail, consider constants t a ∈ R, x a , ẋa ∈ R n , τ a ∈ R m representing the actual time, pose, velocity and cable tensions, respectively. Considering that the cable tensions τ a are kept constant during a controller cycle period ∆t, the solutions x : R → R n and ẋ : R → R n of the Initial Value Problem (IVP)

ẍ = M(x) -1 g(x) + W(x) τ a -C(x, ẋ) ẋ , (4.3a) x(t a ) = x a (4.3b) ẋ(t a ) = ẋa , (4.3c) lead to a transition mapping φ x : R n × R n × R m → R n such that x(t a + ∆t) ẋ(t a + ∆t) = φ x x a , ẋa , τ a . (4.4)
Using a more compact notation, the vector of the next pose and velocity

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

(4.5)

Typically, the OCP formulation of an NMPC scheme and the subsequent stability analysis consider that the transition mapping representing the dynamic system is known. Accordingly, Section 4.2 considers that the expression for φ y is known. The definition of the numerical discretization method is closely related to the definition of the NLP corresponding to an OCP. A deeper discussion on this matter is postponed to the Section 4.3. This section also discusses the items (iii) and (iv). Before that, the reader should be introduced to the proposed NMPC algorithm. Section 4.2 introduces an NMPC taking the same state vector as the one introduced in Chapter 3. This state vector includes the CDPR platform pose and velocity. Since this set constitutes the minimal set of states with which the dynamic system may be modeled with a state space formulation, the corresponding NMPC scheme is called minimal state NMPC.

Minimal State NMPC

The present section introduces a minimal state NMPC scheme and is organized as follows: Section 4.2.1 details the control algorithm, Section 4.2.2 introduces the optimal control problem. The stability of the proposed scheme is discussed in the Section 4.2.3.

NMPC algorithm

The operation of the NMPC scheme proposed in this chapter may be presented with a block diagram equivalent to the used in Chapter 3 (i.e. Figure 3.1). The difference here consists in the MPC formulation which is present in the corresponding block in Figure 3.1. Instead of a linear MPC, the present chapter proposes an NMPC scheme. In accordance with the definition of receding horizon MPC schemes discussed in Section 1.7.2, an optimal control problem is solved every controller cycle such that a cost functional is minimized satisfying system constraints.

As for the control scheme proposed in Chapter 3, no terminal condition is used (in accordance with the reasons discussed in Section 1.7). Accordingly, the cost functional (denoted by J hp : N × Y × S m hp → R) consists of a sum of individual costs associated to each sampling time within the prediction horizon h p ∈ N. The function responsible to compute these individual costs is called stage cost function and is denoted by : N × Y × U → R. Therefore, for an instant k ∈ N, actual state y a ∈ Y and sequence of cable tensions τ (•) ∈ S m hp , the cost functional is given by The implemented NMPC algorithm is in accordance with the standard formulation of NMPC schemes without terminal conditions. The Algorithm 1 recalls the common implementation of NMPC schemes with the notations used in this chapter. This algorithm is also used to define the control policy τ f b : N × Y → U, which computes the vector of desired cable tensions τ f b (k, y a ) for given instant k ∈ N and actual state y a ∈ Y.

J hp (k, y a , τ (•)) = hp j=1 k, y τ (•) (j, y a ), τ j , ( 4
Considering a stage cost : N × Y × U (which will be defined in Section 4.2.2), this vector is computed based on the solution of the OCP (4.8). Algorithm 1 integrates this control policy within the feedback loop such that the OCP (4.8) is solved for each controller cycle with updated estimations of the states. Accordingly, the NMPC scheme defines a set of desired cable tensions τ d (t) to be applied in the time interval t ∈ [k ∆t, (k + 1) ∆t). These outputs are used as setpoints for a cable tension control scheme.

The estimation of the states y k in step 3 is performed with the tools introduced in Chapter 2. The detailed definition of the OCP (4.8) (mainly the considered stage cost) is presented in the Section 4.2.2. Solve the OCP

Algorithm 1 NMPC Algorithm

Main control loop

τ * ← arg min τ (•)∈U hp hp j=1 j + k, y τ (•) (j, y a ), τ j s. t. y τ (•) (j, y a ) ∈ Y ∀ j ∈ N 1,hp (4.8) 3: Set τ v f b ← τ * 1 ; 4:
return τ v f b . 5: end function

OCP formulation

The definition of a pertinent stage cost is a crucial step for the design of a stable and effective NMPC scheme. In order to meet the theoretical requirements necessary to guarantee stability and satisfying performance of an NMPC scheme, one should seek some particular properties of the stage cost. Postponing the formal discussion on this matter to Section 4.2.3, some key properties may be sketched to explain the rationale used in the design of the proposed OCP. The next paragraphs are dedicated to this end.

Typically, NMPC schemes penalize states and control inputs with respect to the desired behavior of the system. In order to better formulate this desired behavior, the following assumption defines the sequence of desired cable tensions and presents the conditions under which they are feasible. represents the desired trajectory. The sequence y d (•) is considered feasible, in the sense that there exists a sequence of desired cable tensions τ d (•) ∈ S m such that, for y 0 = y d [START_REF] Gibney | Nature -Gigantic Chinese telescope opens to astronomers worldwide[END_REF],

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

τ ∈ U | W (x d (k)) τ d (k) = W (x d (k)) τ and τ < τ d (k) . (4.10)
The cost functional of a tracking MPC (or classical stabilizing MPC) scheme should be null for states and control inputs equal to the desired trajectory and desired control inputs, and strictly positive if the states are not equal to the desired ones. More precisely, the stage cost function should satisfy

(k, y, τ ) = 0 if y = y d (k) and τ = τ d (k), (k, y, τ ) > 0 if y = y d (k). (4.11)
NMPC schemes satisfying this property are called hereinafter in this thesis as Standard Tracking NMPC (ST-NMPC). As a matter of fact, the formulation proposed in Section 3.1 is not in accordance with the standard MPC form. Since the stage cost used in Chapter 3 presents a term proportional to the squared norm of the vector of cable tensions (the term (3.19)), which is always strictly positive due to the minimal cable tension bounds, the corresponding OCP does not have null cost for trajectories coinciding with the desired one. Therefore, this section proposes a different approach aiming to minimize the 2-norm of the cable tension vector while satisfying (4.11).

As for the desired states, ST-NMPC schemes commonly penalize a weighted norm of the difference between control inputs taken as argument of the stage cost and the desired control inputs. For the stage cost (k, y, τ ) studied here, a weighted norm of ττ d (k) could be used. Nevertheless, in order to obtain τ d (•), one would need to solve the tension distribution problem of the whole trajectory beforehand, using a TDA coupled to the inverse dynamics. Although this is a reasonable solution, the OCPs proposed in this thesis exempt the control scheme of this previous computation. To this end, the following stage cost is proposed

(k, y, τ ) = y -y d (k) 2 Ky + τ -τr (x y , τ ) 2 Kτ , (4.12)
where

K y = diag(k y ) and K τ = diag(k τ ), with constant vectors k y ∈ R 2n + and k τ ∈ R m + . (•) K is the norm of (•)
weighted with a symmetric positive definite matrix K. The vector x y represents the pose vector obtained by extracting the first n elements of the y. The vectorial function τr : X × U → U denotes the estimation of desired cable tensions, such that, for given x ∈ X and τ in ∈ U,

τr (x, τ in ) = arg min τ τ 2 2 (4.13a) s. t. W(x) τ in = W(x) τ (4.13b) τ min τ τ max (4.13c)
In words, τr (x, τ in ) is the vector of cable tensions with minimal 2-norm that is able to generate the same wrench than τ in at the pose x. Therefore, for given k ∈ N, y ∈ Y and τ in ∈ U, the term ττr (x y , τ in ) 2 Kτ = 0 iff τ presents minimal 2-norm (τ = τr (x y , τ in )). Otherwise, ττr (x y , τ in ) 2

Kτ > 0, indicating that there exists a τr (x y , τ in ) ∈ U generating the same wrench (thanks to (4.13b)) with τr (x y , τ in ) 2 < τ 2 .

The Algorithm 2 presented in Section 4.3 is able to solve the optimization problem (4.13). As discussed in Section 1.5, one crucial part of a TDA is to compute a feasible cable tension vector able to generate the desired wrench. Since τ in ∈ U, this vector is necessarily feasible and τ min τ in τ max . Moreover, the desired wrench in (4.13b) is, in fact, the one generated by τ in . Therefore, the cable vector of tensions τ in may be used as a feasible initial guess in an iterative optimization, which may be much simpler than a general TDA.

For a given actual state y a at time instant k the proposed optimal control problem is summarized as follows

min τ (•)∈U hp hp j=1 j + k, y τ (•) (j, y a ) , τ i . s. t. y τ (•) j, y a ∈ Y ∀ j ∈ N 1,hp (4.14)
Note that, as a consequence of (4.10) and (4.13),

τr (x d (k), τ d (k)) -τ d (k) = k, y d (k), τ d (k) = 0 ∀ k ∈ N (4.15)
and conditions (4.11) are satisfied with the stage cost (4.12). Thereby, the tension distribution is implicitly performed when solving the OCP (4.14).

Stability analysis

This section presents the stability analysis of the closed-loop system obtained with the NMPC proposed in Sections 4.2.1 and 4.2.2. More precisely, the main goal is to study the conditions under which the corresponding closed-loop system is uniformly asymptotically stable according to the following definition.

y j+1 = φ y (y j , τ f b (j, y j )) , ∀ j ∈ N (4.16)
is called uniformly asymptotically stable to y d (•) on Y if, for each y a ∈ Y and k ∈ N, there exists a β ∈ KL such that y f b (•) ∈ S 2n defined according

y f b 1 = φ y (y a , τ f b (k, y a )) y f b j+1 = φ y y f b j , τ f b (j + k, y f b j ) , ∀ j ∈ N (4.17)
satisfies the following relation

y f b j y d (k+j) β( y a y d (k) , j), ∀ j ∈ N. (4.18)
In words, (4.17) constructs the sequence of states y f b (•) that is obtained considering the closed-loop (4.16) "departing" from y a at sampling time k. The closed-loop (4.16) takes as feedback control policy τ f b defined in Algorithm 1. The definition of uniform asymptotic stability presented above is based on the upper bound of the error

y f b j -y d (k + j) = y f b j y d (k+j) .
Note that, since the sequence y f b (•) starts shifted from k sampling periods from the initial time, the error between y f b (•) and y d (•) should be shifted in time as well. In other words, k represents the initial discrete time of the studied trajectory and j represents its evolution in time starting from sampling time k. The same notations are used in the remainder of this section for trajectories others than y f b (•).

The upper bound used to limit the tracking error is established in terms of the function β ∈ KL. An illustrative example of a function belonging to this class is presented in Figure 4.1. In accordance to the definition presented in Section 4.1.1, the illustrative function β ex (r , t) tends to zero for t → ∞ and is strictly decreasing with respect to the second argument (for a constant r ∈ R + ). Since the second argument in (4.18) 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 a y d (k) , the function of time β( y a y d (k) , •) ∈ L may be increased as well. Definition 4.2 is used in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Definition 2.16] for the analysis of MPC schemes and is closely related to [62, 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 used notations are similar to those proposed in this reference.

Since the OCP (4.14) does not use stabilizing terminal conditions, tools designed to the analysis of NMPC schemes without terminal conditions should be used in the stability analysis of the closed-loop (4.16). The minimal time-varying stage cost denoted by * : N × Y → R is fundamental in this context. This function is defined as

* (k, y) = min τ ∈U (k, y, τ ). (4.19)
The asymptotic stability of ST-NMPCs without terminal conditions is addressed in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 6]. Summarizing some of the results presented in this chapter, the following theorem can be stated.

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

k + j, y τ e (•) (j, y a ), τ e j C σ j * (k, y a ) , (4.21)
for all k, j ∈ N. Then, the nominal closed-loop (4. [START_REF] Sato | Development of string-based force display: SPIDAR[END_REF]) is uniformly asymptotically stable on Y provided that h p is sufficiently large. Theorem 4.3 describes sufficient conditions to obtain a stable closed-loop system with the proposed NMPC for h p sufficiently large. Inequalities (4.20) are closely related to (4.11). These conditions will be quite straightforwardly deduced at the end of this section.

Nevertheless, inequality (4.21) represents an important restriction to the properties of the stage cost. This inequality imposes, for each discrete instant k ∈ N and actual state y a ∈ Y, the existence of a feasible sequence τ e (•) that leads to a stage cost converging exponentially to zero in time. The term k + j, y τ e (•) (j, y a ), τ e j represents this stage cost. The exponential convergence is a result of the term C σ j with σ ∈ (0, 1). Note that the upper bound function C σ j * (k, y a ) is proportional to the initial minimal stage cost * (k, y a ).

It is important to highlight that Theorem 4.3 is based on the existence of any feasi-ble sequence τ e (•) that satisfies (4.21). The control policy τ f b is not considered explicitly in this theorem. In contrast to most of the control architectures, constrained MPC schemes generally do not have an analytical expression for the control policy. For this reason, the closed-loop cannot be studied explicitly. This is done by means of suitable properties of the stage cost. In this chapter, this suitable property is described by (4.21).

Since one of the terms in the stage cost (4.12) is the error yy d 2 Ky , this error need to be bounded as well. Based on that, the following assumption establishes important properties of the controlled system in order to attain (4.21). For each y a ∈ Y and k ∈ N, there exist

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

If Assumption 4.4 is true, then, for each state y a ∈ Y and instant k ∈ N, there exists a feasible sequence of cable tensions τ (•) that generates a trajectory in which the error y τ (•) (k + j, y a )y d (k + j) Ky exponentially converges to zero. Moreover, the upper bound function in (4.23a) is proportional to the initial errors e y (k, y a ). Note that, for given k and y a , the error e y (k, y a ) is a constant with respect to the discrete time j.

It is easy to show that this assumption may be satisfied with a typical control scheme. Consider, for instance, a continuous-time OSC scheme that applies the wrench

f = M(x) ẍd + K P (x d -x) + K D ( ẋd -x) + C(x, ẋ) ẋ -g(x) (4.24)
with diagonal gain matrices K P and K D . The nominal closed-loop system dynamics is given by ẍ

= ẍd + K P (x d -x) + K D ( ẋd -x) (4.25)
or, written in terms of the error e(t)

= x d (t) -x(t), ë = -K D ė -K P e. ( 4.26)
Choosing, for instance, K P = K D = I, the solution to the differential equation (4.26) can be written as e(t) = e -t/2 (a 1 sin (ω t) + a 2 cos (ω t)) , (

with ω = √ 3/2 and constant vectors a 1 , a 2 ∈ R n given by

a 1 = e 0 a 2 = 1 ω ė0 - 1 2 e 0 (4.28)
where e 0 = e(0), ė0 = ė(0). The expressions of the error e and its time derivative satisfy

e(t) = e -t/2 a 1 sin(ω t) + a 2 cos(ω t) e -1/2 t ( a 1 + a 2 ) ė(t) = e -t/2 1 2 a 1 -ω a 2 a3 sin(ω t) + 1 2 a 2 + ω a 1 a4 cos(ω t) e -1/2 t ( a 3 + a 4 )
A tedious but straightforward computation using (4.28) leads to the bound e(t) 2 + ė(t) 2 c ey e -t ( e 0 2 + ė0 2 ), with c ey = (16

+ 2 √ 3)/3. Define k x , k ẋ ∈ R n subvectors of the previously defined k y such that k y = k T x k T ẋ T 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 ) , K min y = min i (k y,i ) , ( 4.29)
the following inequalities hold

e(t) 2 Kx + ė(t) 2 K ẋ K max y ( e(t) 2 + ė(t) 2 ) e -1 t K max y c ey ( e 0 2 + ė2 0 ) e -1 t K max y K min y c ey ( e 0 2 Kx + ė0 2 K ẋ )
which is the continuous counterpart of (4.23a) with

C 1 = c ey K max y /K min y < ∞ and σ 1 = e -1 < 1.
As important as (4.23a) imposing the exponential convergence of error e y (k, •), (4.23b) requires that the sequence τ (•) is feasible. Considering an OSC strategy, this means that the wrench defined in the feedback 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 4.4. The study of the feasibility of a given wrench considering different poses and velocities is addressed in several works (e.g. [START_REF] Gouttefarde | Analysis of the wrench-closure workspace of planar parallel cable-driven mechanisms[END_REF][START_REF] Pierre | On the workspace of suspended cable-driven parallel robots[END_REF][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 it is out of the scope of this thesis. Moreover, it is important to highlight that assumption (4.23b) is often implicitly taken. Studies such as [START_REF] Begey | Dynamic Control of Parallel Robots Driven by Flexible Cables and Actuated by Position-Controlled Winches[END_REF][START_REF] Shang | Synchronization Control in the Cable Space for Cable-Driven Parallel Robots[END_REF][START_REF] El-Ghazaly | Adaptive terminal sliding mode control of a redundantly-actuated cable-driven parallel manipulator: CoGiRo[END_REF] analyze the corresponding closed-loop systems considering that the feedback loop does not lead to unfeasible cable tensions.

Since the term e y (k, y) = yy d (k) 2 Ky in the expression (4.12) 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 ( τ -τr (x y , τ ) 2 Kτ ). (4.30)
Lemma 4.7 will show that the minimum min τ ∈U ( ττr (x y , τ ) Kτ ) is known and (4.30) can thus be simplified. Before stating and proving Lemma 4.7, a preparatory lemma is necessary.

Proof. First, the existence of a τr (x, τ in ) for all x ∈ R n , τ in ∈ U should be proved. The objective function τ 2 2 in (4.13a) is quadratic and strictly convex. The substitution of τ = τ in ∈ U in the constraints of (4.13b)-(4.13c) shows that τ in itself is an element of the set defined by these constraints. Hence, this is a convex non-empty set. Therefore, (4.13) is a feasible strictly convex inequality constrained QP problem, which assumes a global minimum.

Denoting f ∈ R n the wrench such that f = W(x) τ a , the definition (4.13) indicates that both τr (x, τ a ) and τr (x, τ b ) are obtained with

τr (x, τ a ) = τr (x, τ b ) = arg min τ τ 2 2 , ( 4.33)
s. t. W(x) τ = f (4.34) τ min τ τ max (4.35) and, therefore, (4.32).
Lemma 4.9 shows that if two vectors of cable tensions generate the same wrench in a given pose, they lead to the same estimation of desired cable tensions. The following corollary is based on this conclusion. Corollary 4.6. For every τ a , τ b ∈ U and x ∈ R, if Proof. For every y ∈ Y and τ in ∈ U, one may set τ out = τr (x y , τ in ). Thanks to the constraints (4.13c), τ out ∈ U. Moreover, τ out -τr (x y , τ in ) Kτ = 0. Since the expression τ out -τr (x, τ in ) Kτ = 0 is non-negative, this proves (4.38). Substituting (4.38) in (4.30), equation (4.39) is obtained.

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

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

y + = φ y (y, τ b ) = φ y (y, τ a ). (4.41)
Limitations resulting from this Assumption 4.8 are better discussed in Section 4.3. This assumption is also commonly used in the state-of-the-art control schemes presented in Section 1.4.

Based on Assumption 4.8, for a given initial state vector, two sequences of cable tensions that generate identical wrenches on the corresponding poses result in identical trajectories. Therefore, any trajectory performed with an arbitrary sequence of feasible cable tensions may be also generated with an alternative sequence consisting of minimal 2-norm vectors. In this case, the stage cost related to each instant along the trajectory is equal to the minimal time-varying cost. This assertion is rigorously formulated in Lemma 4.9. Lemma 4.9. For each y a ∈ Y, y(•) ∈ S 2n and τ (•) ∈ S m such that, for all k ∈ N,

y k = y τ (•) (k, y a ) ∈ Y (4.42a) τ k ∈ U, (4.42b)
there exists a τ * (•) ∈ S m that satisfies

y τ (•) (k, y a ) = y τ * (•) (k, y a ) = y k , (4.43a) τr (x y k , τ * k ) = τ * k ∈ U (4.43b) (k, y k , τ * k ) = e y (k, y k ) = * (k, y k ) (4.43c)
for all k ∈ N.

Proof. For y a , y(•) and τ (•) stated in the lemma, define τ * (•) ∈ S m and y * (•) ∈ S 2n according to

τ * k = τr (x y k , τ k ), y * k = y τ * (•) (k, y a ) (4.44)
for all k ∈ N. Due to the constraint (4.13b),

W(x ya ) τ 1 = W(x ya ) τ * 1 (4.45)
and, according to Assumption 4.8,

y * 1 = φ y (y a , τ * 1 ) = y τ * (•) (1, y a ) = y 1 . ( 4
(k, y k , τ * k ) = e y (k, y k ) + τ * k -τr (x y k , τ * k ) Kτ = e y (k, y k ) = * (k, y k ), (4.48)
as in (4.43c).

Finally, the main results regarding the stability of the ST-NMPC scheme proposed in Section 4.2.1 may be presented as follows.

Theorem 4.10. Consider the ST-NMPC Algorithm 1 with prediction horizon h p ∈ N and feasible desired trajectory y d (•) ∈ S 2n . Consider also that Assumptions 4.4 and 4.8 hold. Then, the nominal closed-loop system (4. [START_REF] Sato | Development of string-based force display: SPIDAR[END_REF]) is uniformly asymptotically stable on Y provided that h p is sufficiently large.

Proof. The proof consists of the analysis of the conditions presented in Theorem 4.3. 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. According to Lemma 4.7, * (k, y) = y -y d (k) 2

Ky

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 Ky K max y y -y d (k) 2 . ( 4.49)
where K min y and K max y were defined in (4.29).

Defining

α 1 , α 2 ∈ K ∞ by α 1 (r) = K min y r 2 and α 2 (r) = K max y r 2 , (4.49) can be written as α 1 ( y y d (k) ) * (k, y) α 2 ( y y d (k) ) ∀ k ∈ N. (4.50)
In addition, (4.21) should be proved. If Assumption 4.4 holds, then for each y a ∈ Y and k ∈ N, there exist τ (•) ∈ S m , real C 1 , < ∞ and σ 1 ∈ (0, 1) such that

e y k + j, y τ (•) (j, y a ) C 1 σ j 1 e y (k, y a ) (4.51a) τ j ∈ U (4.51b)
for all j ∈ N. Applying Lemma 4.9 with τ (•) = τ (•), the obtained sequence τ * (•) leads to the same trajectory (as implied by (4.43a)) and, therefore, (4.51a) remains valid. In addition, (4.43c) implies that

j + k, y τ * (•) (j, y a ), τ * j = e y (j + k, y τ * (•) (j, y a )) (4.52a) C 1 σ j 1 e y (k, y a ) (4.52b) = C 1 σ j 1 * (k, y a ), (4.52c)
where (4.51a) and (4.39) were used in order to obtain (4.52b) and (4.52c), respectively. Inequality (4.52) proves (4.21), finalizing the stability analysis.

Numerical Implementation of the Control Algorithms

The present section introduces the main numerical algorithms necessary to the implementation of the control scheme discussed in Section 4.2. Section 4.3.1 proposes a numerical discretization method, Section 4.3.2 introduces the algorithm to compute τr and the numerical solution of the optimal control problem is discussed in Section 4.3.3.

Numerical discretization

As discussed in Section 4.1, a discrete-time system approximating the continuoustime model (3.10) is necessary in order to implement the proposed NMPC.

As for general nonlinear systems, it is not possible to compute the exact solution of the initial value problem (4.3) (represented by the transition mapping φ y in (4.5)). Therefore, one may define an estimation φy of the exact transition mapping φ y . In practice, a discrete-time system is obtained using some kind of numerical integration method such as the Euler's or Runge-Kutta's methods [START_REF] Edwards | Differential equations and boundary value problems[END_REF]. The control scheme proposed in Section 4.2 may be implemented with any of these integration algorithms.

Recalling the IVP (4.3), explicit numerical integration methods return ỹ+ = φ(y a , τ a ) based on the computation of ẍ according to (4.3a) for suitably chosen combinations of x and ẋ in the neighborhood of x a and x a . Calling n φ the number of considered combinations of x and ẋ, the integration error may be reduced increasing n φ . This is illustrated in Figure 4.2 comparing the approximation given by the Euler's method for n φ,1 = 3 and n φ,2 = 6. The resulting error for n φ,2 is approximately the half of the one obtained for n φ,1 .

The maximal error obtained with the approximate transition mapping is denoted by ∈ R 0 + , such that φ y (y, τ ) -φy (y, τ )

∀ y ∈ Y, τ ∈ U. (4.53)
It should be noted that the stability analysis presented in the previous sections neglects the error . In order to prove the stability of the closed-loop system considering this divergence, one may analyze the robustness of the control scheme and consider the numerical discretization errors as disturbances. This is indeed, a possible follow-up of the studies presented in this thesis.

Besides the number of evaluations n φ , the error is also influenced by the time step ∆t and the implemented numerical discretization method. Generally, shorter time steps lead to smaller discretization errors. Figure 4.2 illustrates this fact considering a reduced ∆t . Clearly, the original accumulated error (obtained for n φ,2 = 6 and ∆t) is obtained for ∆t = ∆t/3 and n φ,2 = 2.

Similarly, the numerical integration method should be chosen in order to obtain a good compromise between precision and computational burden. Figure 4.3 compares the Euler's and Runge-Kutta's methods. One may note that improved results are obtained with Runge-Kutta's method. Indeed, Runge-Kutta with n φ = 4 (known as classical Runge-Kutta's method) is frequently used for the numerical solution of ordinary differential equations. For most of the applications, this represents a good compromise between computation time and numerical precision. The interested reader is referred to the large amount of textbooks addressing the numerical solution of ordinary differential equations, e.g. [START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I, Nonstiff Problems[END_REF].

In spite of the advantages obtained with advanced numerical discretization methods, some specific concerns related to the proposed NMPC scheme should be evaluated. Firstly, the sampling period used for the position tracking of CDPRs is shorter than typical applications of MPC schemes. Section 1.7 recalled the example of ∆t = 1 min used in [START_REF] Hidalgo | Nonlinear model predictive control of styrene polymerization at unstable operating points[END_REF], whereas the control schemes implemented in this thesis use sampling periods close to 10 ms. Moreover, the approximate transition mapping φ should be integrated in the OCP formulation (4.14) in order to obtain a standard NLP problem (this is discussed in detail in Section 4.3.3.1). Accordingly, intricate formulations of φ may hinder the numerical solution of (4.14).

Therefore, a rather simple method is proposed in this thesis. Similarly to Chapter 3, n φ = 1 is chosen. In this case, Euler's and Runge-Kutta's methods are equivalent. For an actual time t a , the derivatives of x are considered constant within the sampling period, i.e. ẋ(t) ≈ ẋ(t a ), ẍ(t) ≈ ẍ(t a ) for t ∈ R ta,ta+∆t and

x(t a + ∆t) ≈ x(t a ) + ∆t ẋ(t a ) ẋ(t a + ∆t) ≈ ẋ(t a ) + ∆t ẍ(t a ). (4.54)
Hence, the transition mapping introduced in Chapter 3 is also used here: Nevertheless, a general transition mapping does not necessarily satisfy this assumption. Similarly to the influence of the error , the presented stability results are valid for systems in which φ(y, τ ) -φ(y, τ ) are smaller than the disturbance magnitude established by the robustness of the NMPC controller. The same aforementioned robustness study should be done in order to quantify this magnitude.

Estimation of the desired cable tensions τr

The estimation of the desired cable tensions defined according to the optimization (4.13) plays an important role in the stability analysis of the proposed control schemes. The present section introduces an algorithm able to solve this optimization problem.

More in detail, the efficiency of the proposed NMPCs depends on the fast computation of τr (x, τ ). For this reason, special attention is dedicated to some details involved in the numerical solution of (4.13).

Denoting A ineq = I -I Taking τ = τ in , (4.13b) (and (4.59b)) are satisfied. Moreover, in accordance with the domain declared in (4.13), τ in ∈ U. As a result, (4.13c) (and (4.59c)) are satisfied substituting τ in = τ . The cable tension vector τ in is, therefore, in the feasible region defined by (4.59b)-(4.59c) and may be used as a feasible initial guess in an iterative active-set method [START_REF] Nocedal | Numerical optimization[END_REF]. Figure 4.4 depicts an illustrative example of the methodology implemented with the active-set. The feasible region is determined by three linear inequalities, represented by A, B and C. Taking as initial guess the vector p 1 , the algorithm is able to compute the global optimum p * within 3 iterations. In each iteration, an equality constrained optimization problem is solved considering some of the inequalities as equalities. In a given iteration, constraints that are considered as equalities are called active constraints. The remaining constraints, which are neglected, are called inactive constraints. Taking the initial guess p 1 in the first iteration, the inequality A is considered as an equality and inequalities B and C are neglected (inactive). The vector p * 1 is obtained solving this subproblem and is taken as input for the next iteration.

One may conclude that the cost function may be decreased by setting inequality constraint A as inactive (this is done analyzing the Lagrangian multipliers related to constraint A, considered as equality in iteration 1). Therefore, iteration 2 solves a subproblem neglecting this constraint. The obtained solution is p * 2 . However, this vector is not within the feasible region. Accordingly, the next iteration considers p 3 = p 2 + α 2 (p * 2p 2 ), with some α 2 ∈ [0, 1]. Since an additional constraint is faced, iteration 3 considers B as an equality constraint. The solution of the optimization problem taking B as equality constraint returns the vector p * 3 . As in iteration 2, the global optimum is obtained as p * = p 3 + α 3 (p * 3p 3 ).

Similarly, the active-set method solving (4.59) computes a τ * k in an iteration k as the solution of an equality-constrained QP problem in which some of the inequalities (4.59c) are set as equalities 1 . Denoting f in = W(x) τ in , the sub-problems considered in an iteration k has the following form

τ * k = arg min τ τ T τ (4.60a) s. t.       W(x) a a 1,k a a 2,k . . .       A eq,k τ =       f in b a 1,k b a 2,k . . .       b eq,k (4.60b)
Vectors a a i,k are row vectors of A ineq , which have one term equal to 1 or -1 and the remaining elements equal to 0. Similarly, b a i,k are elements of b ineq . The indices a 1,k , a 2,k , ... are given by the set of active constraints A k = {a 1,k , a 2,k , . . . , a n a,k } ⊂ N 1,2m , with n a,k the number of active inequalities in iteration k. The overall number of equalities in (4.60) is n e,k = n a,k + n. Clearly, considering that τ min,i < τ max,i ∀ i ∈ N 1,m , if a constraint τ i = τ min,i is active, then τ i = τ max,i (its counterpart for τ max,i is not active) and vice-versa. Moreover, assuming that matrices A eq,k are full rank, constraints (4.60b) may be seen as a system of n e,k linearly independent equations with m variables. Therefore, n e,k m.

The optimization problem (4.60) is strictly convex with an Hessian matrix equals to the identity matrix. A vector τ ∈ R m satisfying the first-order optimality conditions of this equality constrained problem is a global minimum [START_REF] Nocedal | Numerical optimization[END_REF]Theorem 16.2]. The firstorder optimality conditions (KKT conditions) for (4.60) can be written as τ + A T eq,k µ k = 0 (4.61a)

A eq,k τ = b eq,k , (

where µ k ∈ R n e,k is the vector of Lagrange multipliers.

In order to solve (4.61), the matrix A T eq,k may be factorized with QR decomposition such that

A T eq,k = Y k Z k R k 0 , ( 4.62)
with Y k , Z k orthonormal and R k square and upper triangular.

According to the definition of Y k and Z k , column vectors of Y k Z k form an orthonormal basis for R m . Therefore, any vector in R m may be written as a linear combination of the column vectors of this matrix. Denoting the number of columns of Y k and Z k as n y = n e,k and n z = m -n e,k , respectively, this implies that there exists τ y ∈ R ny and τ z ∈ R nz such that the solution of (4.60) may be written as

τ * k = Y k τ y + Z k τ z (4.63)
Left multiplying (4.61a) by Z T k and recalling that Since R k is upper triangular, vector τ y can be computed with back-substitution in (4.66). One may note that, once Y k , R k and Z k are defined, the solution τ * k = Y k τ y + Z k τ z = Y k τ y of the sub-problem (4.60) can be obtained solving (4.66) without great computational effort. Indeed, the same solution may be obtained with

Z T k A T eq,k = Z T k Y k = 0, Z T k Z k τ z = 0 ⇒ τ z = 0. ( 4
τ * k = A † eq,k b eq,k (4.67)
with A † eq,k the pseudo-inverse of A eq,k . Nevertheless, the back-substitution in (4.66) demands significantly less computational effort than the computation of A † eq,k . Details on the computation of Y k , Z k and R k are further discussed below.

Since only some of the inequalities (4.59c) are considered in (4.60), τ * k may violate the cable tension limits. Hence, the definition of a feasible τ k+1 is done as illustrated in Figure 4.5. Denoting q ∈ R m with q = τ * k -τ k , some α k ∈ [0, 1] may be defined so that τ min τ k+1 = τ k + α k q τ max . More in detail, it is interesting to set α k as large as possible within [0, 1] without violating the cable tension limits.

To this end, one may evaluate the intersection of each hyperplane

Π i = {τ ∈ U | a i τ = b i }, (4.68)
for i ∈ N 1,2m , and the half-line

L = {τ ∈ R m | ∃ α ∈ R 0 + : τ = τ k + α q}. (4.69)
Half-line L and hyperplanes Π i are illustrated in Figure 4.5 for the case m = 2. Clearly, for a j ∈ N 1,m , if q j > 0, the half-line L does not intersect Π j and intersects Π j+m (hyperplanes defining the lower bound and upper bound of τ j , respectively). Using the example in Figure 4.5,

q 1 = τ * k,1 -τ k,1 > 0 such that Π 1 ∩ L = ∅ and Π 3 ∩ L = {τ pl 1 }.
Conversely, if q j < 0, the half-line L intersects Π j and does not intersect Π j+m . In Figure 4.5, Π 2 ∩ L = {τ pl 2 } and Π 4 ∩ L = ∅. Moreover, if q j = 0, the intersection of Π j (or Π j+m ) with L is either empty or the half-line L itself. Therefore, Π j and Π j+m are not considered in the computation of α k if q j = 0.

Thereby, for each pair Π j , Π j+m (upper and lower bound for τ j ) with q j = 0, only one intersection point is computed. Accordingly, the elements of the vector τ b ∈ U are defined as

τ b,j = τ max,j if q j > 0 τ min,j otherwise ∀ j ∈ N 1,m , (4.70)
and the intersection of Π j (or Π j+m ) with q j = 0 may be expressed as τ pl,j = τ k + τ b,j -τ k,j q j α j q.

(4.71)

If every α j > 1, then τ * k ∈ U and τ k+1 = τ * k = τ k + q. Otherwise, the minimal value of α j should be taken. In this case, the inequality corresponding to the minimal α j (called blocking constraint) is added to the set of active constraints. The index of the inequality in (4.59c) representing the blocking constraint is denoted by j b .

More precisely, the scalar α k is computed according to

α k = min 1, min j∈I k τ b,j -τ j q j , ( 4.72)
with I k ⊂ N 1,m defined as

I k = {j ∈ N 1,m | q j = 0}. (4.73)
Once α k is computed, the next iteration takes the feasible τ k+1 = τ k + α k q ∈ U and, if a blocking constraint is met, the corresponding inequality is added as equality in the set A k+1 .

In a given iteration k, a possible solution of (4.60) may be τ * k = τ k itself. In this case, q = 0 and the computation of α k is inappropriate (I k = ∅). Then, if the current τ k is not optimal, the active set A k should be updated. Possibly, dropping one equality constraint from the set of active constraints may allow a further reduction of the cost function without violating the cable tension limits. The Lagrange multipliers in (4.61a) contains sufficient information to perform the evaluation of which equality constraint should be dropped.

The vector of Lagrangian multipliers µ k may be decomposed as µ k = µ T eq,k µ T ineq,k T with µ eq,k ∈ R n and µ ineq,k ∈ R n a,k . The vector µ eq,k is related to the constraints

W(x) τ = f in , ( 4.74)
while the vector µ ineq,k is related to the constraints

    a a 1 ,k a a 2 ,k . . .     τ =     b a 1 ,k b a 2 ,k . . .     . (4.75)
According to [START_REF] Nocedal | Numerical optimization[END_REF]Theorem 16.4], since the Hessian of (4.59) is positive definite, a τ satisfying the first order optimality conditions of (4.59) is the global minimum of this QP problem. Moreover, the vector τ k is optimal if every element of µ ineq,k is nonnegative for an iteration returning τ k+1 = τ * k = τ k [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5]. On the other hand, if one of its elements is negative, this indicates that dropping the corresponding inequality from A leads to a lower value of the cost function without violating the cable tension limits. Accordingly, denoting by j min the index of the inequality that leads to the (negative) minimal element of µ ineq,k , the active set considered in the next iteration is taken as

A k+1 = A k \{j min }. (4.76)
However, it is necessary to compute µ ineq,k in order to determine j min . The substitution of (4.63) Vector

µ k = µ T eq,k µ T ineq,k
T can be computed performing back-substitution of (4.78). Moreover, whereas the vector µ ineq,k is necessary to decide the set A k+1 , vector µ eq,k is unimportant in the numerical algorithm. Therefore, one may perform the back-substitution only for the last n a,k lines of the system (4.78) (the lines corresponding to µ ineq,k ).

Summarizing the handling of the active set for two subsequent iterations (k and k + 1) with non-optimal τ k+1 , three actions may be taken. Namely,

• If the solution τ * k (with τ * k = τ k ) is feasible, then A k+1 = A k ; • If τ * k / ∈ U, the blocking constraint j b is added such that A k+1 = A k ∪ {j b }; • If τ * k = τ k , then j min is dropped such that A k+1 = A k \{j min }.
The Algorithm 2 summarizes the proposed procedure to compute the estimation of desired cable tensions. Algorithm 2 is equivalent to the active set method proposed in [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5]. This algorithm identifies the global minimum of a strictly convex problem in a finite number of iterations. Considering the particular form of (4.59), some adaptations with respect to [START_REF] Nocedal | Numerical optimization[END_REF]Algorithm 16.3] have been made in order to improve the computation time.

Indeed, regarding the complexity of implementation, solutions simpler than Algorithm 2 could be obtained. Besides the explicit solution (4.67), alternative simpler methods such as the projected gradient method could lead to good convergence. Nevertheless, we believe that the proposed method presents a reduced computation time. A rigorous evaluation of different numerical solutions of (4.59) may be the subject of future studies.

For each iteration of states 6 to 29 in Algorithm 2, the computation time related to the determination of τ * k , q and µ ineq,k is small. On the other hand, the QR decomposition necessary to the computation of Y k , Z k and R k may lead to significant computational burden. To lower this burden, the fact that at most one inequality is added or suppressed from two successive A k and A k+1 may be used so that Y k , Z k and R k do not need to be computed from scratch. , the equation above may be rewritten as

The computation of Y

A T eq,k+1 = Y k Z k R k u j b 0 v j b . (4.81)
It is possible to define an orthonormal Q such that v r = Q v j b has all its elements equal to zero except for the first element. More precisely,

Q v j b = v j b 0 = v r . ( 4
A T eq,k+1 = Y k Z k I 0 0 QT Q R k u j b 0 v r R = Q R (4.83)
Noting that Q is orthonormal and R is square and upper triangular, expression (4.83) may be rewritten as

A T eq,k+1 = Y k+1 Z k+1 R k+1 0 , ( 4.84)
which determines the updated Y k+1 , Z k+1 and R k+1 .

In case (ii), an equality is removed from the active set. The matrix A T eq,k+1 is obtained removing one column, say column j r , of A T eq,k . As a consequence, the corresponding column of R k is also removed and A T eq,k+1 is given by

A T eq,k+1 = Y k Z k                             r 1,
0 0 0 0 0 0 0                             R , ( 4.85)
where r i,j is the element in the i th row and j th column of R k .

Even though R is not upper triangular, the set of Householder transformations necessary to render R upper rectangular is cheaper than the QR factorization of A T eq,k+1 computed from scratch. Hence, an orthonormal Q can be determined by means of Householder transformations such that

R = Q R k+1 0 , ( 4.86)
with upper triangular R k+1 . Then, (4.85) may be expressed as

A T eq,k+1 = Y k Z k Q Q R k+1 0 R = Q R , ( 4.87)
which defines Y k+1 , Z k+1 and R k+1 .

OCP numerical solution

Section 4.2 presented an NMPC scheme considering the existence of a numerical algorithm able to solve the corresponding OCP. The present section introduces a solution using the numerical integration discussed in subsection 4.3.1.

NLP Formulation

The discretization (4.55) should be integrated within the OCP formulation so that, for given actual state y a ∈ Y and sequence of cable tensions τ (•) ∈ S m hp , the states y τ (•) (i, y a ) can be explicitly considered within the cost function of an NLP. This subject is addressed in detail in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 11]. A straightforward method is to add the sequence ỹ ∈ S 2n hp of states within the prediction horizon as argument of the NLP and consider the transition mapping (4.55) as equality constraints. This approach leads to the following optimization problem

ỹj = φ(ỹ j-1 , τ j ), ∀ j ∈ N 1,hp (4.88c) ỹj ∈ Y, ∀ j ∈ N 1,hp (4.88d) τ j ∈ U, ∀ j ∈ N 1,hp (4.88e)
Even though (4.88) is a standard NLP problem, its solution is hindered by the increased number of arguments and constraints.

In order to overcome this drawback, the expression (4.55) can be applied recursively (as in (4.7)) so that states y τ (•) (i, y a ) can be explicitly written in function of the actual state y a and the sequence τ (•). More precisely, the recursive application of (4.55) leads to ỹ0 = y a , where, E and Ψ are given by

ỹj = y τ (•) (j, y a ) = A j y a + j i=1 A j-i B(ỹ j-1 ) τ i + v(ỹ j-1 ) , ∀ j ∈ N 1,hp . ( 4
Ψ = ψ T 1 ψ T 2 ... ψ T hp T (4.91a) E =          E 1,1 0 0 . . . 0 E 2,1 E 2,2 0 . . . 0 E 3,1 E 3,2 E 3,
         , ( 4.91b)
with

ψ j = A j y a + j i=1 A j-i v(ỹ i-1 )
E j,i = A j-i B(ỹ j-1 ). Considering the definition of stage cost (4.12), the cost functional in (4.14) may be written as

hp i=1 i + k,y τ (•) (i, y a ) , τ i = = γ -γ d,k 2 Kγ + ũd -u 2 Ku = (E u + Ψ -γ d,k ) T K γ (E u + Ψ -γ d,k ) + (u -ũd ) T K u (u -ũd ) = u T H ỹ,k u + 2 h T ỹ,k u + c γ , (4.94)
with

H ỹ,k = E T K γ E + K u , h ỹ,k = (Ψ -γ d,k ) T K γ E -ũT d K u , c γ = (Ψ -γ d,k ) T K γ (Ψ -γ d,k ) + ũT d K u ũd , ( 4.95)
and K γ ∈ R hp2n×hp2n , K u ∈ R hpm×hpm given by

K γ = diag k T y k T y . . . k T y T , K u = diag k T τ k T τ . . . k T τ T . ( 4
ỹj ∈ Y ∀ j ∈ N 1,hp (4.97d) τ u (j) ∈ U ∀ j ∈ N 1,hp (4.97e)
where τ u (j) is obtained extracting the j th vector with dimension m from u, such that u = τ u (1) T . . . τ u (h p ) T T . In spite of the similarity between (4.88) and (4.97), these two NLPs present two crucial differences.

Optimization problem (4.97) has a reduced set of arguments (ỹ(•) is not an argument). More precisely, (4.88) has (m + n) h p arguments while (4.97) has m h p . Moreover, the constraint (4.88c) demands particular attention in the numerical solution of (4.88) because it implies that arguments τ (•) and ỹ(•) are consistent with respect to the transition mapping φ. In contrast, constraint (4.97c) is implicitly satisfied. Since the sequence ỹ(•) is not an argument of (4.97), each vector ỹj ∀ j ∈ N 1,hp is attributed according to (4.97c) and, thereby, satisfy these equality constraints. The sequence ỹ(•) in (4.97) is nothing more than an internal variable (neither output nor argument of the optimization problem) used to compute the matrices H ỹ,k , h ỹ,k and to check the constraint (4.97d).

Optimization problem (4.97) represents a standard NLP, which can be solved with standard numerical optimization libraries, e.g. [START_REF] Ford | The NAG library 'machine[END_REF][START_REF] Bochkanov | ALGLIB-a cross-platform numerical analysis and data processing library[END_REF]. Note that matrices H ỹ,k , h ỹ,k are dependent of ỹ(•) and, therefore, are not constant. An optimization problem taking as argument u ∈ R m.hp can be classified as a QP problem only if its cost function can be written as f (u) = u T H u + 2 h T u, with constant matrix H and vector h. Hence, (4.97) is not a QP problem.

Numerical Solution of the NLP problem

Section 4.3.3.1 introduced a standard nonlinear optimization problem representing the OCP proposed in Section 4.2. Since the CDPR prototypes used in this thesis are not compatible with third party numerical libraries, the present section briefly describes the implemented numerical solution of (4.97). Results presented in Section 4.5 were obtained using this implementation.

The proposed algorithm is inspired by the SQP method. Nocedal and Wright [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 18] discuss this class of numerical algorithm in detail. Grüne and Pannek [START_REF] Grüne | Nonlinear model predictive control[END_REF]Chapter 12] study the applicability of such algorithms to NMPC solutions. These are iterative algorithms in which, at each iteration, the original nonlinear problem is approximated as a QP problem. Each of these sub-problems may be solved with a standard QP algorithm. An SQP approach is considered appropriate for the solution of (4.97) since the proposed OCP possesses a quadratic cost function.

As discussed in [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 6], the computation of the second derivatives of a nonlinear cost function is often computationally expensive. For this reason, the quadratic approximation taken for each iteration of an SQP method commonly considers an estimation of the exact hessian matrix. Typically, this approximation is done with quasi-Newton algorithms such as the BFGS or SR1 methods [99, Sections 6.1 and 6.2]. Indeed, the implemented algorithm uses the BFGS method.

Quasi-Newton algorithms determine approximations of the Hessian matrix of a given cost function considering its first derivatives. The direct application of such methods for (4.97) may be possible but would lead to a troublesome implementation. The first derivatives of the cost function in (4.97a) should take into account the variations of ũd with respect to u according to (4.93). The variations of τr with respect to x ∈ X and τ in ∈ U in (4.13) may be not continuous. For this reason, (4.97) is solved iteratively considering, for each iteration, a constant vector ûd = ũd , as described in the following.

Figure 4.6 summarizes the implemented numerical algorithm. Broadly speaking, the algorithm solves iteratively the original NLP problem applying, for each iteration, a standard SQP method for a simpler version of (4.97), which takes a constant ûd instead of ũd .

More in detail, for given actual state y a and initial guess u 0 , the vectors γ 0 and ũd,0 can be computed using (4.90) lations. Based on these derivatives, a standard SQP method can be used, leading to a sequence of cable tensions represented by u 1 . This new vector is used to obtain updated γ 1 and ũd,1 . The procedure is repeated until the convergence of ûd is satisfied.

As described above, for each iteration of the overall algorithm, an SQP algorithm takes place. As indicated by its name, a sequential quadratic programming method solves a sequence of subproblems consisting of QP problems. The Hessian matrix of each one of these QP problems is estimated using the BFGS method. The quadratic programming problems themselves should be solved with a given QP algorithm.

The implemented QP algorithm consists of a combination of the gradient projection and conjugate gradient methods. The gradient projection method (in contrast to the active-set method used for the LMPC in Chapter 3) is known to be particularly suitable for QP problems with box constraints, which is the case of the constraint set U. Indeed, if the state constraint set is given by As discussed in [99, Section 16.7], the convergence of the projected gradient method can be significantly improved if it is coupled with the conjugate gradient method ( [99, Chapter 5]). For this reason, these two methods are combined in the proposed algorithm.

Y = R 2n , ( 4

Numerical Simulations

This section presents numerical simulations comparing the linear MPC proposed in Chapter 3 with the strategy introduced in Section 4.2. Results show that the linear approximation considered in Chapter 3 may lead to poor performances when the system dynamics presents significant nonlinearities.

More precisely, for a given instant k, the influence of changes on y k+i ∀ i ∈ N 1,hp in the approximation (3.15) may deteriorate the performance obtained with the corresponding linear MPC. The error resulting from approximation (3.15) is increased for substantial variations of B(y k+i ) and v(y k+i ). In general, this issue is present for high velocities and accelerations. A crucial imprecision that may be obtained with the approximation (3.15) is that the wrench matrix W is considered constant within the prediction horizon. This matrix may change significantly for high velocities, leading to non-negligible variations of B(y k+i ). In contrast, this section shows that the proposed NMPC overcomes this drawback since it considers the original nonlinear model.

In order to stress the differences between these different MPC approaches, results presented in this section simulate the operation of the CDPR HRPCable (Figure 1.1), which can safely operate with relatively high velocities and accelerations in comparison to the Hephaestus prototype.

As for the previous chapter, a typical pick-and-place trajectory is used in order to represent the operation of the robot. Figure 4.8 illustrates the poses visited in this trajectory. The desired trajectory is generated using fifth degree polynomials in time relating two consecutive poses through a straight line path. Maximal linear velocity and acceleration are 0.7 m/s and 0.26 m/s 2 .

It is interesting to note that the workspace dimensions should be considered when comparing the magnitude of developed velocities between different prototypes. Typically, velocities with similar order of magnitude may lead to higher variations of W (and B, consequently) for CDPRs with smaller workspace such as illustrated in Both linear MPC (proposed Chapter 3) and NMPC (introduced in Section 4.2) were simulated using the trajectory illustrated in Figure 4.8 and cable tension limits τ min = 100 N, τ max = 2000 N. For the linear MPC, the following controller parameters have been used: Approximate values of the gains were tuned manually and used in order to define the lower and upper bounds for each gain being optimized. Sampling period and prediction horizon were defined as 10 ms and 12, respectively. Results are shown in Figure 4.9.

k x = 10 6 1.
Due to the reasons discussed above, the linear MPC strategy leads to deteriorated results. The maximal tracking errors reaches almost 36 mm and 0.17 • . This indicates that a nonlinear approach is necessary.

Similarly to the procedure used to define the gains of the linear MPC, a particle swarm algorithm was used in order to tune the NMPC gains. In this case, the considered objective function computes a weighed sum of the maximal translation and The same sampling period is used (10 ms) and the prediction horizon is taken equal to 6. Figure 4.10 depicts the results obtained with the proposed minimal state NMPC.

The errors depicted in Figure 4.10 indicate that the proposed NMPC effectively tracks the desired trajectory even with the incidence of significant nonlinearities due to relatively high velocities and accelerations. As a matter of fact, the obtained errors are virtually null (less than 0.05 mm and 0.9×10 -4 degrees), as expected for a scenario without the occurrence of disturbances or model uncertainties. It is also interesting to highlight that the prediction horizon used in the NMPC scheme is the half of its counterpart used in the linear MPC (6 versus 12). Moreover, increasing further the prediction horizon of the linear MPC controller does not lead to improved performance since the linear approximation is more susceptible to imprecision for longer prediction horizons. The considered value of 12 was the best compromise found.

As in Section 3.3, it is also pertinent to evaluate the ability of the proposed control scheme to perform the tracking of unfeasible desired trajectories. To this end, a shortened trajectory is studied. The first, second and third poses in Figure 4.8 were As expected, the tracking of the trajectory can be performed without violating the cable tension limits. For time in 7.3 t 12.1 (s) the tension of cable 3 is equal to 500 = τ max,3 and the desired trajectory cannot be tracked without errors. Accordingly, tracking errors ranging up to 35 mm and 1.3 • are obtained. Therefore, one may expect to extend the stability analysis results of Section 4.2.3 for cases in which the feasibility considered in Assumption 4.1 is not satisfied. This is, indeed, a subject currently under study.

It is worth noting that the results presented in this section consider the state constraint set equal to Y = R 2n . As a result, constraint (4.97d) does not affect the solution of (4.97). This indicates that, in practice, the robot workspace may be indirectly delimited by the constraint (4.97e) (cable tension limits). This exempts the control designer from the potentially arduous task of defining the set Y. Moreover, the solution of (4.97) is also facilitated by the reduced number of constraints.

Real-Time Experiments

The NMPC scheme proposed in Section 4.2 was implemented on the HRPCable prototype leading to the experimental results detailed in the present section. The performance of the LMPC and NMPC schemes (proposed respectively in Section 3.1 and 4.2) are compared using typical pick-and-place trajectories. These control schemes are also compared with respect to their robustness against payload uncertainties.

The control parameters were tuned by trial and error and are given as follows:

k x =

Pick-and-place trajectories

Typical pick-and-place trajectories were performed taking the visited poses shown in Figure 4.12. The desired trajectories are defined using 5 th degree polynomials based In accordance with the goals defined in the present chapter, the comparative results presented in Figure 4.14 indicates that the proposed NMPC scheme leads to significantly better tracking errors than the LMPC approach. Indeed, Table 4.1 shows that significant improvements were obtained considering both the RMS of the tracking errors and the maximal errors.

Reduced maximal cable tension

Chapter 3 highlighted that the main contribution of the LMPC scheme proposed in Chapter 3 lies in its capability to handle cable tension limits explicitly. Thanks to this feature, this LMPC is not prone to the feasibility issues discussed in Section 1.7 and exemplified in Figure 3.9. In order to verify this property, an unfeasible desired trajectory is used in Section 3.3. The same behavior should be validated for the NMPC scheme proposed in the present chapter.

Therefore, as for Section 3.3.2.1, a second scenario is considered taking a maximal cable tension reduced to τ max = 250 N. The obtained results are shown in Figure 4.15. The whole trajectory could be performed without interruption and the restricted cable tension limits were satisfied. Therefore, this experiment demonstrates that the proposed NMPC scheme is not prone to the feasibility issues depicted in Figure 3.9.

Robustness against payload uncertainties

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 3.10. An additional payload of 15.1 kg was used, representing 65.6% of the platform mass. The obtained results are depicted on Figure 4.16. The same procedure was performed using the LMPC proposed in Chapter 3. Comparative results are presented in Figure 4.17. Table 4.2 shows that, as for the results presented in Section 4.5.1.1, significantly better tracking errors were obtained with the proposed NMPC scheme.

Conclusions

This chapter discussed the main drawbacks of the LMPC scheme proposed in Chapter 3. Namely, this approach hindered the stability analysis of the resulting closed-loop system and it may lead to deteriorated tracking errors. To overcome these drawbacks, an NMPC is introduced based on a nonlinear discrete time model and a new cost func-tional formulation. The stability of the corresponding closed-loop system can be analyzed. Sufficient conditions leading to uniform asymptotic stability were deduced. Numerical simulations and real-time experiments validated the applicability of the proposed NMPC scheme. Significantly better tracking errors were obtained for a typical pick-and-place trajectory. The proposed NMPC scheme also led to improved precision regarding its robustness against payload uncertainties. A video summarizing the main experiments performed with the proposed NMPC scheme can be watched here. posed velocity-based scheme is able to tackle these issues with a non-model based approach. Accordingly, the proposed strategy may be implemented without the identification of the dynamic parameters of the winches. Details on the numerical and experimental validation are presented. This approach proves to be appropriate for CD-PRs presenting winches with large reduction ratio, which leads to significant friction and model uncertainties. Using the aforementioned kinematic model and the velocitybased cable tension control, the implementation of different MPC position tracking control schemes is enabled.

Assuming limited incidence of nonlinearities on the CDPR dynamic system, Chapter 3 introduces a linear model able to estimate the mobile platform motion within the prediction horizon. Based on this approximation, a linear MPC scheme is proposed. Numerical simulations and experimental tests proved that the linear MPC may safely operate close to system constraints. This capability is validated applying a desired trajectory that cannot be performed without violating the cable tension limits. The proposed linear MPC scheme is able to perform a trajectory as close as possible to the desired one while satisfying cable tension bounds. Conversely, state-of-the-art control schemes are not able to suitably respond under such conditions. Comparing the behavior obtained with the proposed linear MPC and a state-of-the-art control scheme, one may conclude that the capability to operate close to system constraints represents an important result related to the safety of operation of CDPRs.

Additionally, this achievement may enable improvements on the design of CD-PRs. Typically, acknowledging the risks taken when operating close to system constraints, the design of a CDPR should lead to oversized dimensions. The proposed MPC scheme, being able to safely operate in these conditions, enables the designer to reduce such oversized dimensions.

Nevertheless, it was noted that the proposed linear MPC may be sensitive to increased nonlinearities. The positioning precision may be deteriorated for trajectories presenting relatively high velocities. Accordingly, Chapter 4 introduces an NMPC able to consider system nonlinearities. In contrast to its linear counterpart, the proposed NMPC can be classified as a standard tracking NMPC, and the stability of the resulting closed-loop system could be analyzed. In order to fit with the required characteristics of standard tracking NMPCs, an alternative method for the minimization of the cable tensions is proposed. Details on its numerical implementation are presented. The improved performance is validated experimentally and with numerical simulations.

After the construction of the Hephaestus CDPR prototype, the very first tests using this prototype were performed on November 2019. Acceptable results were obtained using a kinematic control (discussed in Section 1.4.1). For this reason, the consortium dedicated low priority to the implementation of advanced control strategies with respect to other primordial tasks. Accordingly, the time available (one week) for the implementation of the proposed MPC schemes on the Hephaestus prototype was not sufficient to obtain significant results. Such implementation remains a focus of possible future works.

It is worth recalling that all the experimental results presented in this thesis were obtained using software and hardware commonly used in the industry. This demon-strates the applicability of the proposed methods in commercial applications.

Perspectives

Satisfying results were obtained with the numerical algorithms proposed in this thesis. Nevertheless, as highlighted in the course of this manuscript, many of them should be subjected to a rigorous analysis in order to be thoroughly validated. This is the case, for instance, of the numerical solution of the NLP (4.97) proposed in Section 4.3.3.2. Even though the algorithm described in Section 4.3.3.2 led to satisfying results in the studied CDPRs and trajectories, a thorough evaluation of its convergence capabilities may validate its applicability in a general case. Similarly, the Algorithm 2 used to compute τr and the numerical solution of QP problems proposed in Section A.3.2 should be compared to other numerical methods in order to evaluate if the computation time can be reduced.

In addition to the studies presented in this manuscript, some other works were initiated and may be pursued in future works. An MPC scheme based on path following was proposed. Whereas a trajectory defines desired poses for each time instant, a path defines a set of poses that are used as reference without any connection with time. Therefore, in contrast to position tracking schemes, this path following MPC does not consider a specific desired state for each sampling time. The main goal is to attain a final desired pose as fast as possible following a predefined path. The path is defined as a curve in an n dimensional space (with n the number of DoF of the CDPR) representing the platform poses.

Additionally, preliminary results were obtained using the proposed cable tension control in order to operate a CDPR with an admittance control scheme. Thereby, the robot pose may be passively controlled by external efforts applied on its mobile platform. Moreover, hybrid position-force control was also successfully tested. Some degrees of freedom may be tracked with high stiffness whereas the remaining degrees of freedom are driven by the influence of external wrenches. This approach may be helpful to perform CDPR aided assembly tasks. Tests using the HRPCable prototype indicates that the platform may be positioned by hand while gravitational effects are compensated.

The reader may note that an important result obtained with the control schemes proposed in Chapters 3 and 4 is the capability to track unfeasible trajectories. Nevertheless, Assumptions 4.1 and 4.4 are not satisfied in this case. Accordingly, improved versions of the NMPC proposed in Chapter 4 could be examined in order to guarantee stability in cases where Assumptions 4.1 and 4.4 are violated. The corresponding stability analysis would use tools designed to economic MPC schemes. In accordance with the brief overview on EMPC presented in Section 1.7.2, the stage cost of such control schemes does not satisfy (4.11). As also discussed in Section 1.7.2, while results for constant and periodic reference states are well-known [START_REF] Müller | Economic model predictive control without terminal constraints for optimal periodic behavior[END_REF][START_REF] Grüne | Economic receding horizon control without terminal constraints[END_REF][START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF][START_REF] Kohler | On Periodic Dissipativity Notions in Economic Model Predictive Control[END_REF][START_REF] Kohler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF], there are very few studies addressing the case of a general desired trajectory [START_REF] Dong | Analysis of economic model predictive control with terminal penalty functions on generalized optimal regimes of operation[END_REF][START_REF] Martin | Dissipativity and economic model predictive control for optimal set operation[END_REF]. Therefore, a first stability analysis using the theory addressing EMPC would probably consider a time-invariant system (with a constant desired pose).

A.1 Motivation

Limitations of memory allocation in real-time are difficult to manage and to predict [START_REF] Puaut | Real-time performance of dynamic memory allocation algorithms[END_REF]. For this reason, it is preferable to allocate the necessary memory before the real-time operation. Accordingly, this is also recommended for TwinCAT real-time applications [190]. Hence, variables are often declared using static allocation in robotic applications. Nevertheless, the implementation of the proposed MPC schemes would most probably be troublesome with this approach. This is illustrated in the following.

Generally, the use of interpreted programming languages (such as MATLAB and Python) or libraries of matrix manipulation for compiled programming languages (such as Eigen [191] for C++) enables the user to numerically manipulate matrices in a way similar to the notation commonly used for symbolic matrix operations. For instance, using one of the aforementioned programming solutions, one may compute C as the sum of two known matrices A and B with one line of code similar to the following script: where the notation defined in MATLAB was used for illustrative purpose.

One may note that this approach may significantly facilitate the implementation of numerical algorithms using complex matrix operations. Nevertheless, dynamic memory allocation is a crucial step for each one of the operators. Using the illustrative operation (A.II), an equivalent C++ code with static memory allocation can be represented as the following: • Matrix Unit(Matrix A): Calculates the unit vector corresponding to each column of the input matrix;

• Basic IIR and FIR filters;

• Rotations matrices and their derivatives in the group SO(3).

As a matter of fact, the operation (A.II) can be performed with the implemented library in TwinCAT with the line of code presented below: In contrast to the matrices used by the operators presented above, a Matrix may present a fixed allocated memory. This is the case for most of the matrices used in more elaborate implemented functionnalities. For instance, matrices used in the dynamic model of the CDPR (such as M, C, and W) are preallocated before the real-time operation.

Besides operations with matrices containing floating point numbers, other functionalities dealing with matrices containing integers as elements (class IntMatrix) were also implemented. A short list of implemented features is presented below:

• Columns and row operations: access, add or delete non-consecutive columns or rows;

• Intersection and Union of two sets represented by two vectors;

• GetLines(B, A, I1): Assign to matrix B the lines I1 of the matrix A;

• iMinVect(vin) and iMaxVect(vin): Returns the index of the minimal and maximal member of a vector;

• FindZerosVector: Finds elements that has absolute value smaller than a threshold.

Thanks to the tools obtained with MatLib, an additional module able to simulate the control scheme was implemented in C++ (compatible with TwinCAT). This module simulates the closed-loop system representing the CDPR behavior. This significantly facilitates the debugging during the implementation.

A.3 Advanced Matrix Computation

In order to implement the proposed MPC and other control strategies, basic features introduced in Section A.2 were used as a basis to develop methods of matrix algebra. The theory and the implementation in C++ (compatible with TwinCAT) of these methods are briefly discussed in this section. The theory used here is detailed in [START_REF] Gene | Matrix computations[END_REF][START_REF] Nocedal | Numerical optimization[END_REF][START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF].

A.3.1 Matrix factorization

Several features of matrix algebra are frequently based on matrix factorization. For instance, the solution of linear systems (determined, undetermined or overdetermined), matrix inversion, calculation of the null space, the Moore-Penrose inverse and the Singular Value Decomposition (SVD) are typically based on matrix factorization. This section discusses three important methods: LU, Cholesky and QR. All of them present advantages and drawbacks relevant in the context of this thesis. Accordingly, they were implemented based on [START_REF] Gene | Matrix computations[END_REF][START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF].

For a given nonsingular square matrix A, the aim of the LU decomposition is to find an upper triangular square matrix U and a lower traingular matrix L so that Π A = L U, with Π a permutation matrix resulted from pivoting. This decomposition is particulaly useful for the solution of determined linear systems. Using LU factorization, a linear set may be decomposed as

A x = b ⇒ L U x = Π -1 b (A.1)
which can be solved by

L y = Π -1 b (A.2a) U x = y.
(A.2b)

Since both matrices L and U are triangular, equations (A.2) may be solved with backsubstitution (firstly (A.2a) and then (A.2b) using the obtained y). LU decomposition was implemented using pivoting and Crout's algorithm. This method is known as the preferrefed way to solve a linear set of equations [START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF].

If the matrix of interest is symmetric and positive definite, a more convenient factorization may be used, namely, the Cholesky decomposition. In this case, matrix A may be decomposed as A = L L T and the solution of the linear system may be computed similarly to (A.1)-(A.2). Generally, this method is approximatelly two times faster than a conventional LU decomposition and presents good numerical stability without the need of pivoting. Symmetric positive definite matrices are often present in convex optimization problems. For this reason, Cholesky decomposition was used as described in the Section A.3.2.

Additionally, the QR factorization also presents advantages that were important in several steps in this thesis. Here, a rectangular matrix B is decomposed as B = Q R with Q square orthonormal and R upper triangular. This decomposition was implemented using Householder reductions as described in [START_REF] William H Press | Numerical Recipes in C ++ The Art of Scientific Computing Second Edition[END_REF]Section 11. Note that, since no permutation matrix is present in (A.3), pivoting is not necessary. However, other applications of the QR factorization motivated its application in this of iterations in order to obtain the solution of (A.6) in comparison to the conventional method discussed in Section 4.3.2 (and illustrated in Figure 4.5). The probable cause for this difference is illustrated in Figure A.1. Whereas the conventional method described in Section 4.3.2 would demand two iterations in order to obtain u, this vector is obtained in the first iteration using (A.8). This difference related to the necessary number of iterations was noticed experimentally and in numerical simulations. However, the finite termination guaranteed by a conventional method [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5] is not valid in this case. The vector u is used as input for the next iteration of the active-set method changing the set of active constraints, similarly to the algorithm described in Section 4.3.2.

B.1 Tension Distribution Algorithm

Considering the tension distribution problem of a CDPR with n degrees of freedom driven by m cables (m > n), this section presents a closed-form solution for the minimization of the infinity norm of the cable tensions vector. More precisely, the proposed non-iterative algorithm is able to solve (B.1):

min τ ∈R m τ ∞ s. t. W τ = f τ min τ (B.1)
for given wrench matrix W ∈ R n×m , desired wrench f ∈ R n and minimal cable tension τ min (it is assumed that τ min,i = τ min , τ max,i = τ max ∀ i ∈ N 1,m in order to simplify the presentation of the algorithm). Note that the definition of the infinity norm leads to

τ ∞ = max i (τ i ).
The proposed algorithm is based on the theory of wrench feasibility discussed in [START_REF] Bouchard | On the Ability of a Cable-Driven Robot to Generate a Prescribed Set of Wrenches[END_REF][START_REF] Gouttefarde | Characterization of Parallel Manipulator Available Wrench Set Facets[END_REF]. The tension distribution constraints imposed for cable tensions τ are

W τ = f τ min τ τ max (B.2)
Using the hyperplane shifting method [START_REF] Bouchard | On the Ability of a Cable-Driven Robot to Generate a Prescribed Set of Wrenches[END_REF][START_REF] Gouttefarde | Characterization of Parallel Manipulator Available Wrench Set Facets[END_REF], (B.2) can be written as (B. The τ * max obtained with (B.7) is the minimum value for τ max that may be used in (B.2) for which there is a tension distribution τ * satisfying these constraints. Using the theory of the hyperplane shifting method, it is possible to calculate the tension distribution τ * τ * max which satisfies (B.2). In other words, the definition presented in (B.7) calculates the optimal τ * max and the completion of the TDA demands the calculation of each component of the cable tension vector. The optimality of the algorithm described below is proved in Section B.

Recalling that (B.8) defines τ * + and τ * -, the remaining elements may be calculated using (B.10). In summary:

         τ * + = τ * max τ * -= τ min τ * 0 = W † 0 f -W + τ * + -W -τ * - , (B.11)
where τ * max is defined according to (B.7). For a full-rank W, Section B.2 proves that solution of (B.1) and TDA proposed in (B.11) are both unique and equivalent.

As an example, for a robot with 6 DoF, n = 6. Therefore, once each row of C is composed of a vector spanning the nullspace of five columns of the wrench matrix,

#I 0 h = 5,
with #(•) denoting the cardinality of the set (•). Considering a full rank W, a cable robot with m = 8 cables has #I 0 h = 5 #(I + h ∪ I - h ) = 3.

(B.12)

B.2 Optimality Proof

The present section proves the equivalence between the tension distribution obtained with the general optimization (B.1) and the closed-form solution (B.11).

Consider given W, f and τ min . Take also the matrix C computed as in [START_REF] Gouttefarde | Characterization of Parallel Manipulator Available Wrench Set Facets[END_REF]. Let τ * max be defined according to This proves the two first equations of (B.11). The remaining unknown components τ * 0 should satisfy (B.17):

W 0 τ * 0 = f -W + τ * + -W -τ * -. (B.25)
As stated in the beginning of this section, W is a full-rank matrix. Consequently, W 0 also has full rank. Besides, W 0 has dimensions of (n) × (n -1). Thereby, there is a unique τ * 0 that satisfies (B.25). This vector can be computed as follows

τ * 0 = W † 0 f -W + τ * + -W -

B.3 Experimental Results

The TDA proposed in Section B.1 was implemented in the HRPCable prototype using the programming tools described in Appendix A. This section presents the main experimental results obtained with this algorithm. The trajectory depicted in method (minimization of the infinity norm) and the cable tensions obtained with the minimization of the 2-norm of the cable tensions vector (the method described in [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n + 2 Cables[END_REF] was used). As discussed earlier, (B.11) represents a closed-form expression able to solve the tension distribution problem in the control of CDPRs. This is an important advantage of the proposed method. Nevertheless, the main drawback resulted from the minimization of the infinity norm is the presence of discontinuous desired cable tensions. This drawback is visible in Figure B.3. At t ≈ 25 s, one of the desired cable tensions jumps from 290 N to 230 N. However, this same figure also shows that the control of cable tensions can respond sufficiently fast. At any rate, discontinuities on the desired cable tensions should be avoided. Methods able to overcome this issue are still under investigation.

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 and real constants C 1, σ ∈ (0, 1) such that

1. 1

 1 HRPCable prototype and its CAD view. . . . . . . . . . . . . . . . . . . . 1.2 (a) The five-hundred-meter aperture spherical radio telescope FAST [1, 2] and (b) the prototype built in the framework of the research project CoGiRo (Control of Giant Robots) -TECNALIA/LIRMM [3]. . . . . . . . 1.3 The commercial solution Inca, a cable-driven haptic device from Haption [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 A motion simulator prototype developed in the Max Planck Institute for Biological Cybernetics [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Hephaestus prototype and an illustration of its main application, the installation of curtain wall modules. . . . . . . . . . . . . . . . . . . . . . 1.6 Illustration of the CWM assembly. . . . . . . . . . . . . . . . . . . . . . . 1.7 A typical mechanism on which the drawing pulleys are mounted, illustrating that the drawing point is not fixed. . . . . . . . . . . . . . . . . . .

3. 3

 3 Norm of the (a) translational and (b) rotational errors. . . . . . . . . . . . 3.4 Comparative results: (a) maximal errors, (b) RMS errors, (c) maximal cable tensions, (d) maximal cable tension derivatives, (e) RMS of cable tension derivatives, and (f) consumed energy. . . . . . . . . . . . . . . . . 3.5 Simulated cable tensions for (a) PID+ Controller, (b) Sliding Mode Controller and (c) Model Predictive Control. . . . . . . . . . . . . . . . . . . . 3.6 Illustration of a pick-and-place path in the CAD view of HRPCable. . . . 3.7 Pick-and-place task results with (a) τ max = 400 N and (b) τ max = 260 N. The actual pose is shown in continuous lines and the desired one in dashed lines. Note that the axes limits for TE and OE are different for (a) and (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .List of Tables1.1 Hephaestus and HRPCable main features. . . . . . . . . . . . . . . . . . . 2.1 Estimation of the static dry friction torque (measured values in Nm). . . 4.1 Summary of comparative errors between LMPC and NMPC. . . . . . . . 4.2 Summary of comparative errors between LMPC and NMPC regarding the robustness against payload uncertainties. . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -HRPCable prototype and its CAD view.

Figure 1 . 2 -

 12 Figure 1.2 -(a) The five-hundred-meter aperture spherical radio telescope FAST [1, 2] and (b) the prototype built in the framework of the research project CoGiRo (Control of Giant Robots) -TECNALIA/LIRMM [3].

Figure 1 . 3 -

 13 Figure 1.3 -The commercial solution Inca, a cable-driven haptic device from Haption [4].

Figure 1 . 4 -

 14 Figure 1.4 -A motion simulator prototype developed in the Max Planck Institute for Biological Cybernetics [5].

Figure 1 . 5 -

 15 Figure 1.5 -Hephaestus prototype and an illustration of its main application, the installation of curtain wall modules.

Figure 1 . 6 -

 16 Figure 1.6 -Illustration of the CWM assembly.

Figure 1 . 7 -

 17 Figure 1.7 -A typical mechanism on which the drawing pulleys are mounted, illustrating that the drawing point is not fixed.

Figure 1 . 8 -

 18 Figure 1.8 -The NIST ROBOCRANE, controlled with a Joint Space control [6].

Figure 1 . 11 -

 111 Figure 1.11 -Summary of main studies addressing position tracking of CDPRs.

LociFigure 1 . 13 -

 113 Figure 1.13 -Comparison of Loci with constant ττ p for different values of p.

Figure 1 . 15 -

 115 Figure 1.15 -An illustrative system used to introduce the main elements of an OCP.

Figure 1 . 16 -

 116 Figure 1.16 -A simple example of an application of the Dynamic Programming principle introduced in [10].

Figure 1 . 18 -

 118 Figure 1.18 -Illustration of feasibility issues resulted from the application of terminal constraints under the incidence of disturbances.

Figure 2 . 1 -Figure 2 . 2 -

 2122 Figure 2.1 -Illustration of the main geometric notations.

Figure 2 . 3 -

 23 Figure 2.3 -CAD view of the HRPCable prototype presenting the performed trajectory.

Figure 2 . 4 -

 24 Figure 2.4 -Experimental results: Solid lines represent the poses estimated with the proposed FK algorithm and dashed lines represent the poses measured by the Metris K600 camera system. On the right, the number of iterations necessary for the FK algorithm to converge for each sampling time during the trajectory is presented.

Figure 2 . 5 -

 25 Figure 2.5 -An illustrative schematic view of the experimental set-up.

Figure 2 . 6 -

 26 Figure 2.6 -Drawing of a simplified gear train.

Figure 2 . 7 -

 27 Figure 2.7 -Identification of the static dry friction torque.

Figure 2 . 8 -

 28 Figure 2.8 -Estimation of the kinetic friction coefficients.

Figure 2 . 9 -

 29 Figure 2.9 -Estimation of the kinetic friction coefficients.

Figure 2 . 10 -

 210 Figure 2.10 -Block diagram of the proposed torque control scheme.

Figure 2 . 11 -

 211 Figure 2.11 -Illustration of the simulated simplified system.

Figure 2 . 12 -

 212 Figure 2.12 -Simulation results with the torque control: evolution of the cable tension and the torque versus time.

Figure 2 . 13 -

 213 Figure 2.13 -Simultaneous plot of the applied torque and measured cable tension (with the cable attached to a fixed support).

Figure 2 . 14 -

 214 Figure 2.14 -Simultaneous plot of variations in the unstretched length and the measured cable tension (with the cable attached to a fixed support).

Figure 2 . 15 -

 215 Figure 2.15 -Block diagram of the proposed velocity-based cable tension control scheme.

Figure 2 .Figure 2 . 16 -

 2216 Figure 2.17 shows the results obtained with the velocity-based control with a step input. The lower level is 190 N and the upper is 440 N. The rise time to reach 440

Figure 2 . 17 -

 217 Figure 2.17 -Measured cable tension with a step input in the desired tension.

Figure 2 .Figure 2 . 18 -Figure 2 . 19 -

 2218219 Figure 2.19 depicts the results for a sinusoidal x and constant τ d . This figure compares strategies with and without the feedforward term ẋ/r w . This figure shows that a

Figure 2 . 20 -

 220 Figure 2.20 -Measured and desired cable tensions (on the left) and displacements x (on the right).

Figure 2 . 21 -

 221 Figure 2.21 -Experimental results highlighting the hysteresis present on the behavior of the cable elasticity.

Figure 2 . 22 -

 222 Figure 2.22 -Difference between cable tensions measured by load pins 1 and 2 (τ 1 and τ 2 ).

Figure 2 . 23 -

 223 Figure 2.23 -Identification of the pulley friction: measured cable tension τ 2 (experimental data) and predicted value τ2 given by (2.37) (fitted data).

Contents 3 . 1

 31 Control Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.1.1 Dynamic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.1.2 Control strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.2 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.2.1 Background on the state-of-the-art controllers . . . . . . . . . . 74 3.2.2 Simulations results . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.3 Real-time Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.3.1 Experimental platform and implementation aspects . . . . . . . 79 3.3.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.3.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

  with g the gravitational acceleration. Denoting M(x) = M (x) S(x) and C(x, ẋ) = C (x, ẋ) S(x) + M (x) Ṡ(x, ẋ), (3.9) the dynamic system (3.5) can be rewritten as M(x) ẍ + C(x, ẋ) ẋ = g(x) + W(x) τ c . (3.10)

Figure 3 . 1 -

 31 Figure 3.1 -Block diagram of the proposed control scheme

u 1 2 u

 2 T H c u + d T u s.t. u min u u max (3.25) where u min , u max ∈ R n hc are given by u min = τ T min . . . τ T min T and u max = τ T max . . . control inputs over the control horizon. The vector of desired cable tension is set as τ d (t) = τ * k for (k -1)∆t t < k ∆t. The solution of (3.25) is repeated at each sampling time, updating matrices H c and d in function of y k and τ c (k ∆t).

  . W τ = f (3.32) τ min τ τ max(3.33) Satisfying the constraints above, the tension distribution generates the desired wrench (3.32) with cable tensions in an admissible interval(3.33).

Figure 3 .

 3 Figure 3.3 shows the evolution versus time of the norms of the translational and rotational errors. All the control strategies are able to compensate the tracking error caused by the applied external disturbance, but PID+ presents an oscillatory behavior and an increased settling time. SMC responds faster and without oscillations. MPC presents the fastest response, resulting in the smallest tracking errors along the whole

Figure 3 . 2 -

 32 Figure 3.2 -Illustration of the initial and final positions of the simulated trajectory.

Figure 3 . 3 -

 33 Figure 3.3 -Norm of the (a) translational and (b) rotational errors.

Figure 3 . 4 -

 34 Figure 3.4 -Comparative results: (a) maximal errors, (b) RMS errors, (c) maximal cable tensions, (d) maximal cable tension derivatives, (e) RMS of cable tension derivatives, and (f) consumed energy.

Figure 3 . 5 -

 35 Figure 3.5 -Simulated cable tensions for (a) PID+ Controller, (b) Sliding Mode Controller and (c) Model Predictive Control.

Figure 3 . 6 -

 36 Figure 3.6 -Illustration of a pick-and-place path in the CAD view of HRPCable.

  and positive definite. Diagonal elements are repeated for each set of 12 elements (size of y). These 12 elements are equal to k y = [ k T x k T ẋ ] T = = [ 2 × 10 7 2 × 10 7 3.3 × 10 7 2 × 10 7 . . . 2 × 10 7 2 × 10 7 1.2 × 10 -3 0.6 × 10 -3 . . . 0.6 × 10 -3 0.6 × 10 -3 1.3 × 10 -3 0.6 × 10 -3 ] T . (3.35)

Figure 3 . 7 -

 37 Figure 3.7 -Pick-and-place task results with (a) τ max = 400 N and (b) τ max = 260 N. The actual pose is shown in continuous lines and the desired one in dashed lines. Note that the axes limits for TE and OE are different for (a) and (b).

Figure 3 . 7 -

 37 (a) shows the experimental results obtained with τ max = 400 N and τ min = 100 N. The proposed controller is able to keep the translation errors (TE) smaller than 4.5 mm and the orientation errors (OE) smaller than 0.3 • .Reduced maximum tension:Figure 3.7-(b) shows the experimental results obtained with τ max = 260 N and τ min = 100 N. The proposed MPC finds poses as close as

Figure 3 . 8 -

 38 Figure 3.8 -Overlaid views of both scenarios (nominal constraints and reduced maximum tension). The upper view (a) depicts an instant in which the desired pose is feasible for both cases. Conversely, the lower view (b) depicts an instant in which the desired pose is feasible for τ max = 400 N and unfeasible for τ max = 260 N.

Figure 3 . 9 -Figure 3 . 10 -

 39310 Figure 3.9 -Experimental results obtained with the control scheme proposed in [8] for unfeasible desired an trajectory. Pick the payload

. 6 )

 6 where the state y τ (•) (j, y a ) represents the j th term of the sequence of states y pa (•) ∈ S 2n hp resulting from the application of the cable tensions τ (•) taking as initial state y a . This sequence of states is obtained applying recursively the transition mapping (4.5) j times, such that y pa (1) = y τ (•) (1, y a ) = φ y (y a , τ 1 ) y pa (2) = y τ (•) (2, y a ) = φ y (y pa (1), τ 2 ) y pa (3) = y τ (•) (3, y a ) = φ y (y pa (2), τ 3 )

1 :loop 3 : 4 : 5 :

 1345 Set k ← 0; 2: Estimate the actual states y k ; Set desired cable tensions τ d (t) ← τ f b (k, y k ) for t ∈ [k ∆t, (k + 1) ∆t); Set k ← k + 1; 6: end loop Control policy Inputs: y a ∈ Y and k ∈ N; Output: τ v f b ∈ U; 1: function τ f b (k, y a ) 2:

Assumption 4 . 1 .

 41 Consider the sequences x d (•), ẋd (•) ∈ S n of desired poses and velocities, respectively, and y d (•) ∈ S 2n such that y d (k) = x d (k) T ẋd (k) T T

Definition 4 . 2 .

 42 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

Theorem 4 . 3 .

 43 Consider the NMPC Algorithm 1 with prediction horizon h p ∈ N and minimal stage cost satisfying

Figure 4 . 1 -

 41 Figure 4.1 -An illustrative example of a function β ex ∈ KL.

Assumption 4 . 4 .

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

Lemma 4 . 5 .

 45 Consider τ a , τ b ∈ U and x ∈ X. If W(x) τ a = W(x) τ b , (4.31) then, τr (x, τ a ) = τr (x, τ b ).

  ) then τ b = τr (x, τ b ). (4.37) Proof. Taking τ a , τ b and x from the corollary statement such that τ b = τr (x, τ a ), constraints (4.13b)-(4.13c) in the definition of τr implies that W(x) τ a = W(x) τ b and τ b ∈ U. Therefore, applying Lemma 4.5, τr (x, τ a ) = τr (x, τ b ). Thanks to (4.36), this leads to τ b = τr (x, τ b ).

Lemma 4 . 7 .

 47 For every y ∈ Y, τ in ∈ U and k ∈ N, min τ ∈U ( ττr (x y , τ in ) Kτ ) = 0, (4.38) and the minimal time-varying stage cost is given by * (k, y) = yy d (k) 2 Ky = e y (k, y). (4.39)

Assumption 4 . 8 .

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

Figure 4 . 2 -

 42 Figure 4.2 -Illustration of how the numerical integration precision is influenced by ∆t and n φ for the Euler's method.

Figure 4 . 3 -

 43 Figure 4.3 -Comparison between Euler's and Runge-Kutta's methods for the same n φ .

  x y ) -1 g(x y ) -C(x y , ẋy ) ẋy v(y) . (4.55) It is worth noting that estimation (4.55) satisfies Assumption 4.8, as proved in the following. For the state y ∈ Y and cable tension vectors τ a , τ b ∈ U stated in Assumption 4.8, denote by f ∈ R n the wrench vector given by f = W(x y ) τ a = W(x y ) τ b , (4.56) where (4.40b) was used. The application of (4.55) leads to φy (y, τ b ) = φy (y, τ a ) = A y + 0 ∆t M(x y ) -1 f + v(y), (4.57) which is equivalent to (4.41).This means that the presented stability results are valid for a closed-loop system defined asy k+1 = φy (y k , τ f b (k, y k )) , (4.58)with the NMPC feedback policy τ f b defined in Algorithm 1.

Figure 4 . 4 -

 44 Figure 4.4 -Active-set method: illustrative example of steps taken for the solution of a constrained optimization problem.

. 64 )

 64 Similarly, since A eq,k Z k = 0, (4.63) substituted in (4.61b) leads toA eq,k Y k τ y = b eq,k (4.65) and using (4.62) with the orthogonality of Y k (Y T k Y k = I), R T k τ y = b eq,k . (4.66)

Figure 4 . 5 -

 45 Figure 4.5 -Notations used in order to compute α k .

1 :Set A 1 1 .

 111 function τr (x, τ in ) as the lines of A ineq τ 1 = b ineq that are satisfied; 5: Set I 1 according to (4.73); 6:and the remaining equal 0.Since the matrix Yk Z k is orthonormal, Y k Y T k + Z k Z T k = I and (4.79) may be written asA T eq,k+1 = Y k Z k R k Y T k e m,j b 0 Z T e m,j b . (4.80) Denoting the row vectors of Y k and Z k with Y k = u T

. 82 )

 82 Denoting v b = v j b , (4.82) may be rewritten as v j b = QT v b 0 T and (4.81) may be decomposed as follows

j

  + k, ỹj , τ j (4.88a) s. t. ỹ0 = y a (4.88b)

. 89 )Denoting u = τ T 1 .

 891

( 4 .

 4 91c) Similarly, for given instant k ∈ N, cable tension sequence τ (•) ∈ S m hp and states ỹ ∈ S 2n hp , vectors γ d,k ∈ R 2n hp and ũd ∈ R m hp are defined as γ d,k = y d (k + 1) T . . . y d (k + h p ) T T (4.92) ũd = τr (x ỹ1 , τ 1 ) T . . . τr (x ỹhp , τ hp ) T T . (4.93)

. 96 )

 96 Therefore, for given y a ∈ Y and k ∈ N, the OCP (4.14) may be expressed asu * = arg min u∈R m hp u T H ỹ,k u + 2 h T ỹ,k u (4.97a) s. t. ỹ0 = y a (4.97b) ỹj = φ ỹj-1 , τ u (j) ∀ j ∈ N 1,hp(4.97c)

Figure 4 . 7 -

 47 Figure 4.7 -Comparison of the variation on the cable orientation for illustrative 2-DoF CDPRs with the same platform displacement.

K

  u = k u I = 1.0 × 10 -6 I (4.99c) k ∆u = 0.026. (4.99d)Gains (4.99) were obtained using a particle swarm optimization algorithm (built-in MATLAB function). More in detail, for a given set of controller gains, the considered objective function computes a weighted sum of 1. The maximal translation error; 2. The maximal orientation error; 3. The RMS value of the variation of cable tensions (computed according to Section 3.2).

Figure 4 . 8 - 9 T

 489 Figure 4.8 -Visited poses of the simulated trajectory

Figure 4 . 9 -

 49 Figure 4.9 -Simulation results obtained with the linear MPC proposed in Chapter 3.

Figure 4 . 10 -

 410 Figure 4.10 -Simulation results obtained with the NMPC proposed in Section 4.2.

K

  u = k u I = 6 × 10 -7 I (4.101c) h p = 4.(4.101d)

Figure 4 . 11 -

 411 Figure 4.11 -Application of the proposed NMPC with an unfeasible desired trajectory.

Figure 4 . 12 -

 412 Figure 4.12 -Visited poses for the real-time experiments.

Figure 4 . 13 -

 413 Figure 4.13 -Experimental results for τ max = 400 N.

Figure 4 .

 4 Figure 4.13 depicts the experimental results obtained using the proposed NMPC scheme taking cable tension limits of τ min = 100 N and τ max = 400 N. The same trajectory was performed with the LMPC proposed in Section 3.1.Figure 4.14 compares the corresponding results.

Figure 4 .

 4 Figure 4.13 depicts the experimental results obtained using the proposed NMPC scheme taking cable tension limits of τ min = 100 N and τ max = 400 N. The same trajectory was performed with the LMPC proposed in Section 3.1.Figure 4.14 compares the corresponding results.

Figure 4 . 14 -

 414 Figure 4.14 -Tracking errors for LMPC and NMPC.

Figure 4 . 15 -

 415 Figure 4.15 -Experiments taking τ max = 250 N.

Figure 4 . 16 -

 416 Figure 4.16 -Experimental results obtained with the NMPC scheme addressing the robustness against payload uncertainties.

Figure 4 . 17 -

 417 Figure 4.17 -Tracking errors for the LMPC and NMPC schemes in the experiment testing the robustness against payload uncertainties.

  C = A + B;(A.I)More in detail, in order to compute (A.I), the operator + dynamically allocates the memory equivalent to a matrix with the same size of A (and B), and computes the sum of each elementA[i][j] + B[i][j]. Subsequently, the operator = attributes to C the obtained values and frees the dynamically allocated memory. Using operators similarly defined, more elaborate matrix operations may be performed. For instance, several operations may concatenated such as in the script below: C = (A' -B * 8) * (B + A * B');(A.II)

  transpose(At * ,A * ,na,na); multiplication_by_scalar(B1 * ,B * ,8,na,na); difference(AtB * ,At * ,B1 * ,na,na); transpose(Bt * ,B * ,na,na); multiply(ABt * ,A * ,Bt * ,na,na); sum(BABt * ,B * ,ABt * ,na,na); multiply(C * ,AtB * ,BABt * ,na,na); (A.III) with additional statically allocated matrices At, B1, AtB, Bt, ABt and BABt. These matrices present the same number na of rows and columns than A, B and C. In addition, column;

C

  = (A.trp() -B * 8) * (B + A * B.trp()); (A.IV)

  2] and [88, Section 5.1]. Clearly, the solution of the linear system B x = b can be computed applying backsubstitution for R x = Q T b. (A.3)

3 )= 0 j= i ∈ I | c j w i = 0 (B. 6 )

 3006 of the matrix C is perpendicular to n -1 columns taken among the m columns of W. More specifically, there are n c = m n-1 possible combinations of vectors for a W ∈ R n×m . Denote W k 0 the submatrix formed by each one of these combinations with k ∈ N 1,nc . Each W k 0 results in two lines (say c l and c p ) of C: c T l = null (W k 0 ) T c p = -c l . (B.4) Each line c i corresponds to the orientation of a hyperplane which restricts the feasible wrenches. Vector d translates these hyperplanes. Their components d j are the translation normal to each hyperplane in relation to the origin. These components are given by (B.5) d j = τ max i∈I + j c j w i + τ min i∈I - j c j w i . (B.5) Sets I + j , I - j and I 0 j are subsets of I = N 1,m and are defined according to i ∈ I | c j w i > 0 I - j = i ∈ I | c j w i < 0 I For given τ min , W and f , Hussein et al. [31] propose the minimization of τ max as follows τ * max = max j c j f -i∈I - j τ min c j w i i∈I + j c j w i . (B.7)

2 .

 2 Naming h the argument j that maximizes the expression present in (B.7). The τ * = τ * 1 τ * 2 . . . τ * m T satisfying (B.1) should also satisfy (B.8):if i ∈ I + h τ * i = τ min if i ∈ I - h τ min < τ * i < τ * max if i ∈ I 0 h . (B.8) Therefore, elements {τ * i , i ∈ (I + h ∪ I - h )}are trivially defined using (B.8). The remaining {τ * i , i ∈ I 0 h } should be computed considering the equality constraint in (B.2). Define submatrices and subvectors of W and τ * relating to the sets I + h , I - h and I 0 [w i ], i ∈ I 0 h W + = [w i ], i ∈ I + h W -= [w i ], i ∈ I - .2) impose that W 0 τ * 0 + W + τ * + + W -τ * -= f . (B.10)

c h f = i∈I - h τ min c h w i + i∈I + h τ * max c h w i (B. 14 ) 1 W

 141 τ * max = max j c j f -i∈I - j τ min c j w i i∈I + j c j w i = c h f -i∈I - h τ min c h w i i∈I + h c h w i (B.13)where h is equal to the argument which maximizes (B.13). The equation (B.13) can be written as Suppose that optimization (B.1) obtains a τ ∞ = t b < τ * max . Let d b be the vector calculated using τ max = t b in (B.5). See thatc h f = , h th line of C f < d b is violated. This proves that τ *max is equal to the solution of (B.1). After that, it is necessary to prove that TD proposed in (B.11) satisfies (B.2) with τ max = τ * max .Indeed, substituting τ max = τ * max in (B.2), f is feasible. As a consequence, there is at least one τ * so thatW τ * = f τ min τ * τ * max (B.16)Therefore, using the same notation introduced in Section B.0 τ * 0 + W + τ * + + W -τ * -= f . (B.17)Multiplying (B.17) by c h and using (B.14)c h W + τ * + + c h W -τ * -= i∈I - h τ min c h w i + i∈I + h τ * max c h w i . (B.18)Notice that c h W 0 = 0. Equation (B.18) can then be written asi∈I - h c h w i (τ min -τ * -,i ) + i∈I + h c h w i (τ * max -τ * +,i ) = 0. (B.19)According to (B.6), the definition of I + h and I - h imposes thatc h w i < 0, i ∈ I - h c h w i > 0, i ∈ I +The combination of (B.20) and (B.21) results in i∈I - h c h w i (τ min -τ * -,i ) 0 (B.22) i∈I + h c h w i (τ * max -τ * +,i ) 0. (B.23) According to (B.19), the sum of l.h.s. of (B.22) and (B.23) should be equal to zero. Once c h w i = 0 for i ∈ (I + h ∪ I - h ), τ * -,i = τ min τ * +,i = τ * max ⇒ τ * -= τ min τ * + = τ * max . (B.24)

  τ * -. (B.26) Therefore, the optimization (B.1) has a unique solution which is equal to the TDA (B.11).

Figure B. 1 -

 1 Figure B.1 -CAD view of HRPCable and visited poses.

Figure B. 2 -

 2 Figure B.2 -Performed trajectory.

  Figs. B.1 and B.2 was used in this test.

Figure B. 3

 3 Figure B.3 compares cable tensions obtained with the proposed tension distribution

Figure B. 3 -

 3 Figure B.3 -Desired and measured cable tensions using the minimization of (a) infinity norm and (b) 2-norm.

cM

  y τ (•) (n, y), τ n C σ n * c (y) ∀ n ∈ N 1,N -1 , (C.8)then condition (C.7b) is satisfied withγ K = C 1 -σ K 1 -σ (C.9)andsup K γ K < ∞ ∀ K ∈ N.Condition (C.8) 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 (C.7) are equivalently given as (C.7a) and (C.8). Considering (C.4) (and (C.6), consequently), the conditions of Theorem 4.3 are simplified to (a) There exist α 1 , α 2 ∈ KL ∞ such that α 1 ( yy d (k) ) * c (y) α 2 ( yy d (k) ) ∀ y ∈ Y. (C.10) (b) For all y ∈ Y, there exist feasible τ e (•) ∈ S m , real C < ∞ and σ ∈ (0, 1) satisfying c y τ e (•) (j, y), τ e j C σ j * c (y) , ∀ j ∈ N. (C.11) These conditions are equivalent to (C.7a) and (C.8).The time-varying case addressed in Theorem 4.3 is obtained considering [13, Assumptions 6.29 and 6.30], which indicate that the results above can be straightforwardly extended to time-varying systems. Mass matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m p Mass of the mobile platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . n Number of degrees of freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . p Mobile platform position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ψ Vector of Euler angles representing the orientation of the platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . q Motor position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . r p Pulley radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . τ max Maximal admissible cable tensions . . . . . . . . . . . . . . . . . . . . . . . . τ min Minimal admissible cable tensions . . . . . . . . . . . . . . . . . . . . . . . . . τ c Cable tensions function of time . . . . . . . . . . . . . . . . . . . . . . . . . . . . τ d Desired cable tensions function of time . . . . . . . . . . . . . . . . . . . . τ d Sequence of desired cable tensions . . . . . . . . . . . . . . . . . . . . . . . . . τ (•) Sequence of cable tension vectors . . . . . . . . . . . . . . . . . . . . . . . . . . u Vector of cable tensions within the prediction horizon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . W Wrench matrix .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Mobile platform pose vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x d Desired pose vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x(•) Sequence of platform pose vectors . . . . . . . . . . . . . . . . . . . . . . . . . x Mobile platform pose function of time . . . . . . . . . . . . . . . . . . . . . y Minimal state vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . y + Next state (for given actual state and cable tension) . . . . . . . . 0 Matrix or vector with dimensions suitably chosen and all elements equal to zero . . . . . . . . . . . . . . . . . . . . . . . . . . Functions J : N × Y × U hp → R, cost functional . . . . . . . . . . . . . . . . . . . . . . . . . : N × Y × U, stage cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * : N × Y, minimal time varying stage cost . . . . . . . . . . . . . . . . . . . φ y : Y × U → R 2n transition mapping . . . . . . . . . . . . . . . . . . . . . . . . . R : R 3 → R 3×3 , rotation matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . τ f b : N × Y, NMPC control policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . τr : X × U, estimation of desired cable tensions . . . . . . . . . . . . . . . y τ (•) (•, y a ) : N → R 2n , trajectory for initial state y a and sequence of cable tensions τ (•) . . . . . . . . . . . . . . . . . . . . . . . . . v p , p-norm of a vector v . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v , 2-norm of a vector v . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v K , weighted norm of a vector v given by √ v T K v . . . . . . . . . . . . . v 1 v 2 , Euclidean distance between two vectors v 1 and v 2 . . . . . . . . .

  .12 Illustration that, in general, a cost function ττ 2 (i.e. p = 2) leads to a smaller margin between cable tensions and unfeasible tensions than its counterpart with p = ∞[START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Comparison of Loci with constant ττ p for different values of p. . . 1.14 A simple open-loop cable tension control scheme [8]. . . . . . . . . . . . 1.15 An illustrative system used to introduce the main elements of an OCP. . 1.16 A simple example of an application of the Dynamic Programming principle introduced in [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.17 MPC strategies solve a finite horizon OCP at each control cycle taking updated measured states. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.18 Illustration of feasibility issues resulted from the application of terminal constraints under the incidence of disturbances. . . . . . . . . . . . . . . Identification of the static dry friction torque. . . . . . . . . . . . . . . . . 2.8 Estimation of the kinetic friction coefficients. . . . . . . . . . . . . . . . . 2.9 Estimation of the kinetic friction coefficients. . . . . . . . . . . . . . . . . 2.10 Block diagram of the proposed torque control scheme. . . . . . . . . . . . 2.11 Illustration of the simulated simplified system. . . . . . . . . . . . . . . . Simultaneous plot of variations in the unstretched length and the measured cable tension (with the cable attached to a fixed support). . . . . . 2.15 Block diagram of the proposed velocity-based cable tension control scheme. 60 2.16 Simulation results with the velocity-based control: evolution of the cable tension and the torque versus time. . . . . . . . . . . . . . . . . . . . . . . 2.17 Measured cable tension with a step input in the desired tension. . . . . . 2.18 Measured and desired cable tensions. Displacement x is constant and desired cable tension is determined according τ d = 320+220 cos(2π t/62.5) [N]. One may note that the obtained error is rather small. . . . . . . . . . 2.19 Experimental results for a constant desired cable tension and sinusoidal displacement x. The chart on the left presets the desired and measured tensions comparing schemes with and without the feedforward term ẋ rw . The chart on the right presents the sinusoidal displacement x . . . . . . .

1.11 Summary of main studies addressing position tracking of CDPRs. . . . . 12.1 Illustration of the main geometric notations. . . . . . . . . . . . . . . . . 2.2 Pulley geometry and coordinate systems. . . . . . . . . . . . . . . . . . .

2.3 CAD view of the

HRPCable prototype presenting the performed trajectory. 2.4 Experimental results: Solid lines represent the poses estimated with the proposed FK algorithm and dashed lines represent the poses measured by the Metris K600 camera system. On the right, the number of iterations necessary for the FK algorithm to converge for each sampling time during the trajectory is presented. . . . . . . . . . . . . . . . . . . . . . . . 2.5 An illustrative schematic view of the experimental set-up. . . . . . . . . 2.6 Drawing of a simplified gear train. . . . . . . . . . . . . . . . . . . . . . . 2.7 2.12 Simulation results with the torque control: evolution of the cable tension and the torque versus time. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Simultaneous plot of the applied torque and measured cable tension (with the cable attached to a fixed support). . . . . . . . . . . . . . . . . . ix 2.14 2.20 Measured and desired cable tensions (on the left) and displacements x (on the right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.21 Experimental results highlighting the hysteresis present on the behavior of the cable elasticity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.22 Difference between cable tensions measured by load pins 1 and 2 (τ 1 and τ 2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.23 Identification of the pulley friction: measured cable tension τ 2 (experimental data) and predicted value τ2 given by (2.37) (fitted data). . . . . . 3.1 Block diagram of the proposed control scheme . . . . . . . . . . . . . . . 3.2 Illustration of the initial and final positions of the simulated trajectory. .

Table 1 . 1 -

 11 Hephaestus and HRPCable main features.

		Hephaestus	HRPCable
	Payload	1300 kg	100 kg
	Footprint	14 × 4 × 12 m	8 × 4 × 6 m (x × y × z)
	Number of cables	8	8
	Degrees of freedom	6	6
	Maximum cable tension	16 kN	2 kN
	Diameter of the cables (steel) 12 mm	4 mm
	Diameter of the pulleys	300 mm	80 mm
	Platform mass	23 kg	400 kg
	Motors and Hardware	Beckhoff	Beckhoff
	PLC	C6650 3.6 GHz i7	C6920 2.4 GHz i7
	Motors	AM8071	AM8061
	Gear trains (reduction) Planetary in line (76.2:1) AG2210 (25:1)

Table 2 . 1 -

 21 Estimation of the static dry friction torque (measured values in Nm).

  1 , Z 1 and R 1 is performed with Householder transformations, as proposed in [176, Section 2.10] and detailed in [176, Section 11.2]. Appendix A Algorithm 2 Computation of τr (x, τ ) Inputs: x, τ in . Internal variables: τ min , τ max . Output: τr (x, τ in ).

  1 r 1,2 . . . r 1,jr-1 r 1,jr+1 . . . r 1,n e,k 0 r 2,2 . . . r 2,jr-1 r 2,jr+1 . . . r 2,n e,k

	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	0	0 . . . r jr,jr-1 r jr,jr+1 . . .	. . .
	0	0 . . .	0	r jr+1,jr+1 . . .	. . .
	0	0 . . .	0	r jr+2,jr+1 . . .	. . .
	0	0 . . .	0	0	. . .	. . .
	0	0 . . .	0	0	0 r n e,k ,n e,k
	0	0 . . .	0	0	0	0
	. . .	. . .	. . .	. . .	. . .	. . .	. . .

  and (4.93), respectively. An approximation of the NLP (4.97) is constructed replacing ũd by a constant ûd in (4.94) and (4.95). The resulting cost function possesses continuous first derivatives that can be computed numerically or analytically. Results presented in Section 4.5 were obtained with first derivatives computed analytically, which were defined with tedious (but straightforward) calcu-

					START
				Update values of
	Overall			and	using active-set
	algorithm
					Convergence of	YES
					satis ed?
					NO
	Standard SQP algorithm	conditions satis ed? SQP convergence YES	NO	Estimate the Hessian for constant
					using BFGS
	Standard QP	YES
	algorithm	QP convergence conditions satis ed?	NO	Gradient projection iteration
				Conjugate gradient
					iteration

Optimal

Figure 4.6 -Summary of the implemented numerical solution of (4.97).

  .98) constraints (4.97d) can be neglected and the optimization problem (4.97) is constrained uniquely by the box constraints (4.97e). Section 4.4 discusses that assumption (4.98) may be typically considered. Therefore, Sections 4.4 and 4.5 consider Y = R 2n .

Table 4 . 1 -

 41 Summary of comparative errors between LMPC and NMPC.

		RMS		Maximal error
		TE	OE	TE	OE
	LMPC	4.28 mm 0.15 • 7.47 mm 0.27 •
	NMPC	1.52 mm 0.08 • 2.34 mm 0.14 •
	Improvement	64.5%	44.6%	68.6%	48.1%
	4.5.1.

1 Nominal cable tension limits

Table 4 . 2 -

 42 Summary of comparative errors between LMPC and NMPC regarding the robustness against payload uncertainties.

		RMS		Maximal error
		TE	OE	TE	OE
	LMPC	4.66 mm 0.52 • 12.37 mm 1.05 •
	NMPC	3.92 mm 0.31 •	7.83 mm	0.63 •
	Improvement	16.0%	40.7%	36.7%	40.0%

The used notations are in accordance with those in Figures 1.9 and 1.10.

Even though the proposed numerical solution generates a sequence of cable tensions {τ 1 , τ 2 , . . .}, this sequence does not present any relation with the evolution of tension in time (as it was the case of sequences τ (•) in Section 4.2).

Acknowledgments

Firstly, I'd like to thank my thesis director, Marc Gouttefarde, for his supervision. By the end of this doctorate, I find myself a better professional and a better person. Much of these professional and personal growths were motivated and guided by him. I also thank my co-supervisor, Ahmed Chemori, for his support.

General Conclusion

The present thesis contributes to the state-of-the-art on trajectory position control of CDPRs. In the following, the main obtained results are summarized and an overview of possible future works are presented.

Glossary

Acronyms CDPR Cable-Driven Parallel Robot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CWM Curtain Wall Module . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . DDP Differential Dynamic Programming . . . . . . . . . . . . . . . . . . . . . . . . DoF Degrees of Freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . DPP Dynamic Programming Principle . . . . . . . . . . . . . . . . . . . . . . . . . . EMPC Economic MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . FK Forward Kinematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . GPC Generalized Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . IK Inverse Kinematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . JSC Joint Space Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . KKT Karush-Kuhn-Tucker conditions . . . . . . . . . . . . . . . . . . . . . . . . . . LP Linear Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LQG Linear Quadratic Gaussian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LQR Linear Quadratic Regulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LMPC Linear Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . MPC Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . NLP Nonlinear Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . NMPC Nonlinear MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . OCP Optimal Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . OE Orientation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . OSC Operations Space Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QP Quadratic Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . RHC Receding Horizon Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SMC Sliding Mode Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SQP Sequential Quadratic Programming . . . . . . . . . . . . . . . . . . . . . . . TDA Tension Distribution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . TE Translation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . WCW Wrench Closure Workspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Kinematic Model and Cable Tension Control

Several variables considered in the position tracking control strategies proposed in Chapters 3 and 4 do not have a direct relation with variables that can be directly measured and controlled in the CDPRs studied in this thesis. The variables that can be directly measured and controlled in real-time are the motor currents, positions and velocities. Interestingly, the model on which the control is based presents numerous parameters and variables that are not directly related to the motors operation. The present chapter is devoted to the development of strategies and algorithms relating these variables, which are a priori abstract to the control system, with the variables that can be directly accessed through the motor drives. Section 2.1 introduces an algorithm able to estimate the platform poses and velocities based on the motor positions and velocities. Sections 2.2 and 2.3 discuss two different control strategies responsible to generate the desired cable tensions. In the first one, the control inputs are the motor currents whereas the second one is based on motor velocities.

Pick the payload Pick the payload Release the payload Release the payload able to find a trajectory as close as possible to the desired one while strictly respecting the tension limits. The obtained experimental results also show its robustness against payload uncertainties.

Chapter 1 showed that state-of-the-art control schemes lead to safety issues when operating close to system constraints. The results presented in this section demonstrate that the proposed linear MPC scheme overcome this issue. Therefore, the implemented control algorithm fulfills the main purpose for which it was conceived.

Whereas the necessity to safely operate close to system constraints is a specific requirement in the Hephaestus project, this characteristic is also favorable for the design of CDPRs. Considering the limitations of the existing state-of-the-art schemes, the designer must oversize the margin dedicated to model uncertainties in order to guarantee the safety of the robot operation. In contrast, a more efficient design may be used for a CDPR operating with the proposed MPC.

Nevertheless, the performance of this linear MPC may be deteriorated with the incidence of significant nonlinearities. More precisely, the vector γ obtained with the affine system (3.17) may be a bad approximation of the future states if (3.15)

if q = 0 then Determine index j min of (4.59c) leading to minimal element of µ ineq,k ;

A k+1 ← A k \{j min }; if there is a blocking inequality then

22:

Determine the index j b of (4.59c) that leads to the blocking constraint;

23:

end if

29:

end for 30:

return τ k ; 31: end function briefly discusses this subject. Subsequently, for k > 1, Y k-1 , Z k-1 and R k-1 are used to obtain their updated counterparts (Y k , Z k and R k ). These matrices need to be updated in two possible scenarios: (i) one equality is added to the active set (step 23 in Algorithm 2) or (ii) one equality is removed (step 16 in Algorithm 2).

Consider firstly case (i). In the iteration k, an inequality j b m (an analogous procedure is used for cases with m < j b 2 m) was added such that A T eq,k+1 = A T eq,k e m,j b , (

where e m,j b ∈ R m given by e m,j b = 0 . . . 1 . . . 0 T with the j th b element equal to 1

Summary of the work

The present thesis was developed within the framework of the european union H2020 project named Hephaestus. Chapter 1 introduced the main requirements related to the position tracking of the CDPR prototype used in this project. Besides positioning accuracy, one important requirement is that the CDPR should safely operate close to the system constraints. More in detail, the robot operation includes trajectories demanding nominal cable tensions close to the cable tension limits. Moreover, a significant incidence of disturbances and modeling errors is expected.

In this context, particular attention should be devoted to the cable tension limits. A conventional control scheme may return feedback corrections that are not consistent with the interval of admissible tensions. Indeed, the overview of the state-of-the-art on position tracking of CDPRs presented in Chapter 1 indicates that the existing control schemes are prone to this issue. More precisely, these control schemes operating in the aforementioned conditions expected for the Hephaestus project may either return unfeasible desired cable tensions, or simply crash (attain a condition without a welldefined output). Since MPC is one of the few control strategies able to explicitly handle system constraints, this control design methodology is chosen.

Before addressing the MPC control strategies, some preliminary issues had to be solved. Some of these issues are addressed in Chapter 2. In order to meet high positioning precision, an FK algorithm considering the pulley kinematics is proposed. An explicit expression for the differential kinematics enabled the implementation of a numerical solution of the nonlinear least squares system representing the FK problem. Its convergence capabilities are evaluated experimentally and numerically.

Moreover, the applicability of different closed-loop cable tension control strategies is studied. Such control schemes aim to generate the desired cable tensions using the measurements obtained with force sensors. Firstly, a torque-based control scheme is analyzed. Simulations and experimental results indicate that this approach leads to crucial issues related to the friction present in the gear trains. In contrast, the pro-

Appendices

A

Real Time Implementation -TwinCAT

The application of the proposed MPC strategies demands non-trivial numerical operations and algorithms. Most of these operations consist of matrix manipulations. There are several open source libraries that may be used for this purpose. Nevertheless, the software TwinCAT is not compatible in real-time with third party mathematical libraries. More precisely, TwinCAT uses an alternative mathematical library instead of the basic header math.h. Therefore, C++ codes were developed based on the alternative TwinCAT mathematical library. The present appendix discusses some details of this implementation. Section A.1 compares typical programming approaches that can be used for matrix manipulation, highlighting the constraints imposed by the real-time applications in TwinCAT. Section A.2 briefly introduces the main rationale involved in the implementation of the C++ library. Some details of the implemented matrix factorization methods and optimization algorithms are presented in Section A.3. A * represents the pointer to the first element of matrix A (and similarly to the remaining matrices).

Contents

Clearly, although (A.II) and (A.III) lead to equivalent outputs, the comparison between these scenarios indicates that the implementation of elaborate matrix operations should be facilitated with an approach similar to the one used in (A.II). However, as recommended in [190], memory allocation before the real-time operation should be preferred.

A.2 Basics

The main goal sought in the implementation of the proposed library (called MatLib) is to attain similar programming advantages to the interpreted programming languages without the need of dynamic memory allocation. To this end, a large memory is allocated in the initialization of the control scheme (out of real-time) and MatLib is responsible to explicitly manage this memory. By this means, MatLib enables the user to declare matrices with different sizes using the memory which has been allocated before the real-time operation.

Every operator and function available in this library (and its dependent libraries) are compatible with any 2 dimensional matrix. With the same philosophy, operations with multidimensional matrices are possible. Nevertheless, since the proposed control can be implemented with 2 dimensional matrices, the multidimensional case was not addressed.

Each matrix is declared as an object of the class Matrix. This class presents, among others, the members rows and columns of type int determining the number of rows and columns of the matrix. Additionally, the member p of type double * points to the first element of the corresponding matrix. Once a Matrix is declared, the members p, rows and columns may change.

Besides these members, this class presents several overloaded operators so that the matrix operations may be programmed in an intuitive manner (such as in (A.II)). These operators use an object belonging to the class ProvisoryMatrix. This object has access to a large memory which is allocated before the real-time operation. In order to illustrate the main used rationale, consider the following definition of the operators + and = which are applied between two objects belonging to the class Matrix:

Matrix Matrix::operator+(Matrix M) { MpAllocate(rows, columns, Mp, MpG); for(i = 1; i <= rows; i++) { \\ for (j = 1; j <= columns; j++) { Mp(i, j) = ( * this)(i, j) + M(i, j); } } return Mp; } void Matrix::operator=(Matrix M) { for (i = 1; i < M.rows + 1; i++) { for (j = 1; j < M.columns + 1; j++) { ( * this)(i,j) = M(i, j); } } MpG.free = 0; } Objects Mp and MpG belong, respectively, to the classes Matrix and ProvisoryMatrix. The function MpAllocate(rows, columns, Mp, MpG) reserves the necessary memory in MpG and sets-up Mp to access this memory. This is done setting the pointer Mp.p to the address of the first element within this memory slot and attributing Mp.rows = rows and Mp.columns = columns. The necessary memory reserved is equal to rows * columns numbers of type double. The memory is managed by the variable MpG.free, which indicates how many numbers of type double are being currently used within the memory allocated by MpG. Accordingly, the MpAllocate(rows, columns, Mp, MpG) sets MpG.free = MpG.free + rows * columns. Conversely, this variable is set to 0 whenever the operator = is called, indicating that the memory of managed by MpG may be overwritten. The operator () taking two integers (say i and j) as arguments returns the reference of the element in row i and column j of the corresponding matrix. This way, the returned variable from the operator () may be used to read and write the corresponding element of a matrix.

Using a similar methodology, the common matrix operations were programmed using intuitive operators without the need of dynamic memory allocation. The implemented notation is inspired by the one used in MATLAB. Following, a non-extensive list of implemented trivial matrix operations and declarations is presented:

• Matrix operator() (int i1, int i2, int j1, int j2): Returns a submatrix of rows between i1 and i2, columns between j1 and j2. A call M(i1,i2,j1,j2) is equivalent to M(i1:i2,j1:j2) in MATLAB; • Basic matrix operations +, -and * (sum, difference and multiplication);

• Method Matrix trp(): returns the transposed of the matrix;

• Method Matrix Norm(): returns the matrix with the norm of each column vector;

• Matrix eiFun(int iIn, int mIn): returns a vector with mIn lines with every element equal to zero, except for the iIn th , which is equal to 1;

• Matrix Eye(int iIn, int jIn): returns a iIn×jIn matrix with every element equal to zero, except for those in the main diagonal, which are equal to 1;

• Matrix Zeros(int iIn, int jIn): iIn×jIn zero matrix;

• Matrix Ones(int iIn, int jIn): iIn×jIn matrix with every element equal to 1;

• Matrix Sum(Matrix Min, int Dir): Calculates the sum of each row or thesis. As discussed in [START_REF] Gene | Matrix computations[END_REF], this factorization is particularly efficient for the solution of least squares problems. Accordingly, this method was used in Section 2.1.3 for the solution of the FK problem.

In addition, for B ∈ R n 1 ×n 2 , n 1 > n 2 , the matrices Q and R can be written with submatrices as

with columns(Q 1 ) = rows(R 1 ) (number of columns of Q 1 is equal to the number of rows of R 1 ), it is easy to verify that

This indicates that Q 2 is in the null-space of B T . Indeed, QR factorization may be used in order to compute the null-space of a rectangular matrix. This utility was used in Section 4.3.2. This section also details that matrices Q and R may be updated with reduced computational cost if one row is added or removed from B, which is useful in active-set optimization algorithms. Two methods were considered: active-set and projected gradient [START_REF] Nocedal | Numerical optimization[END_REF]. Projected gradient is simpler to implement and may be pertinent for box constraints such as (A.6b). Nevertheless, active-set finds the optimal solution in a finite number of iterations and can easily consider additional constraints (this later advantage may be used in future works). Therefore, the active-set technique was chosen. This method solves a simplified version of (A.6) taking a subset of the constraints as equality constraints at each iteration (a detilaed discussion addressing this technique is presented in [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5]). The resulting subproblem in a given iteration can be written as of elements of F, n m the number of elements of M and n n the number of elements of N . H, U and d in (A.7) can be written as follows

A.3.2 Quadratic programming

where

and u n ∈ R nn . Note that each element of u m is equal to the corresponding element of u min and each element of u n is equal to the corresponding element of u max . The solution of (A.7) is (u

It is worth noting that H f is symmetric and positive definite. As a consequence, H f may be factorized with Cholesky decomposition such that H f = L L T , with L lower triangular, as described in the last section. Advantages of using this technique were also presented. The optimal u * f is obtained from two successive back substitution procedures:

, and then L T u * f = s.

Before starting the next iteration of the active-set method, components of u * f should comply with cable tension limits. Good results were obtained using a simple saturation function. More precisely, each element u i of the vector u used in the next iteration is obtained as a saturation of the corresponding element u * i of u *

It was noted that, generally, the application of (A.8) reduces the necessary number

C

On the Proof of Theorem 4.3

This appendix discusses how the Theorem 4.3 is obtained based on the results described in [START_REF] Grüne | Nonlinear model predictive control[END_REF].

The general formulation presented in [13, Chapter 6] is applied using the notations introduced in Chapter 4. Theorem 4.3 is based mainly on [START_REF] Grüne | Nonlinear model predictive control[END_REF]Theorem 6.24]. The later addresses timeinvariant closed-loop systems (with a constant desired state vector). Accordingly, first, a simplified case of Theorem 4.3 is considered taking a constant desired state. Sufficient conditions presented in [START_REF] Grüne | Nonlinear model predictive control[END_REF]Theorem 6.24] are reviewed and shown to be equivalent to the conditions of the time-invariant case of Theorem 4.3. The time-varying case is obtained in accordance with the results of [13, Section 6.5].

Theorem 4.3 is based mainly on [START_REF] Grüne | Nonlinear model predictive control[END_REF]Theorem 6.24]. As discussed in Chapter 4, 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 computed based on an instant k and state y according to

with the cost functional J K defined in (4.6). The scalar V K (k, y) is also known as the cost-to-go for given instant k and state y. Note that, if the global minimum for the following optimal control problem exists,

then τ * (•) satisfies the equality

In order to deduce Theorem 4.3, let us first consider a simplified case with a constant desired state

with a constant vector y d ∈ Y. In this case, the closed-loop 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 [START_REF] Grüne | Nonlinear model predictive control[END_REF]Theorem 6.24] may be directly applied as it is. Moreover, the value function V K is constant with respect to time and may be simplified to

Similarly, the stage cost and minimal stage cost can be simplified according to

The conditions stated in [13, Theorem 6.24] address time-invariant systems and are summarized as (i) There exist α 1 , α 2 ∈ K ∞ such that for every y ∈ Y the following relation holds

(C.7a) (ii) ( [13, 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