Vincenzo Di
email: vincenzo.dipaola@edu.unige.it

Alexandre Goldsztejn
email: alexandre.goldsztejn@ls2n.fr

Matteo Zoppi
email: matteo.zoppi@unige.it

St Éphane Caro
email: stephane.caro@ls2n.fr

Design of a Sliding Mode-Adaptive PID Control for Aerial Systems with a Suspended Load Exposed to Wind Gusts

Wind gusts are among the well-known disturbances that afflict aerial systems. A tracking task for a load suspended by massless rigid bars and quadrotors is considered in this paper. The complexity of the tracking task lies in the presence of wind gusts. Therefore, a self-tuning gain techniques for a Proportional-Integral-Derivative (PID) control are introduced to accomplish tracking tasks when wind gusts appear. To automatically tune the PID gains, a sliding mode condition is exploited. The stability of this approach, considering bounded disturbances, is guaranteed by introducing a supervisory control designed through the Lyapunov method. The proposed solution and its performances were tested through simulation and compared with both PD and PID controls.

mass of the load x i , ẋi , ẍi ∈ R 3 position, linear velocity, acceleration of the i th drone in F O

INTRODUCTION

Nowadays quadrotors are increasingly attracting attention from researchers due to their agility and potential ability to accomplish tasks autonomously. These peculiarities make them suitable for a wide range of applications such as filming, grasping, and collaborative transportation [START_REF] Tognon | Theory and Applications for Control of Aerial Robots in Physical Interaction Through Tethers Springer Tracts in Advanced Robotics[END_REF]. The deployment of a team of quadcopters enables them to overcome the payload limits that affect them as individuals. However, the use of multiple quadrotors to steer a load leads to an increase in complexities in terms of control. Within this context, among the control techniques developed so far, both geometric control [START_REF] Bullo | Geometric control of mechanical systems: modeling, analysis, and design for simple mechanical control systems Springer[END_REF] and flatness [START_REF] Jean | Analysis and control of nonlinear systems: A flatness-based approach[END_REF] theories recorded quite large success in the research community. Indeed, in [START_REF] Lee | Geometric control of cooperating multiple quadrotor uavs with a suspended payload[END_REF], a geometric control was designed to command a team of agents transporting a point-mass load with rigid links. Furthermore, in [START_REF] Koushil | Geometric control and differential flatness of a quadrotor uav with a cable-suspended load[END_REF], the systems (i.e. cables-quadrotors-platform) were proven to be flat. Subsequently, these works gave rise to several contributions where flatness and/or geometric control were used to perform aerial transportation tasks with a rigid platform [START_REF] Wu | Geometric control of multiple quadrotors transporting a rigid-body load[END_REF][START_REF] Lee | Geometric control of quadrotor uavs transporting a cable-suspended rigid body[END_REF] and also with flexible cables [START_REF] Taeyoung | Stabilization of a rigid body payload with multiple cooperative quadrotors[END_REF][START_REF] Kotaru | Differential-flatness and control of quadrotor(s) with a payload suspended through flexible cable(s)[END_REF]. Recently, in [START_REF] Li | Cooperative transportation of cable suspended payloads with mavs using monocular vision and inertial sensing[END_REF], the geometric control [START_REF] Lee | Geometric control of quadrotor uavs transporting a cable-suspended rigid body[END_REF] was implemented in a real hardware that exploited the abilities of onboard cameras to estimate the system's state and accomplish tracking tasks. As a result, this work allowed disengaging flights from Motion Capture System (MOCAP). Indeed, one of the main objectives among researchers is to empower aerial systems to be as autonomous as possible in order to let them navigate in outdoor environments autonomously. Hence, a preliminary study investigating how to cope with external disturbances by means of varying the system stiffness was conduct in [START_REF] Di Paola | A preliminary study of factors influencing the stiffness of aerial cable towed systems[END_REF]. With this in mind, this work tackles the issues arising from the presence of wind during aerial transportation tasks. By nature, wind includes turbulence and gusting [START_REF] Cole | Spatio-temporal wind modeling for uav simulations[END_REF]. Consequently, the wind is a stochastic and chaotic phenomenon. Thus, by analogy, the control of chaotic systems [START_REF] Yan | Design of robust controllers for uncertain chaotic systems with nonlinear inputs[END_REF][START_REF] Chang | Adaptive robust pid controller design based on a sliding mode for uncertain chaotic systems[END_REF] was inspirational for this work. In particular, a Sliding-Mode Adaptive PID (SM-APID) control that confers to a system with massless rigid links, quadrotors and a point-mass load, the ability to track a given trajectory under the effect of wind is introduced for the first time. The Dryden wind model [START_REF] Waslander | Wind disturbance estimation and rejection for quadrotor position control[END_REF] was hereby considered to test the SM-APID ability.

Therefore, one contribution of this work is related to the introduction of aerodynamic forces in the dynamical model of aerial systems with massless rigid links, quadrotors and a point mass load. Another major contribution regards the design of the SM-APID control, for the mentioned system, to address the tracking problem when gusts appear. Moreover, numerical simulations demonstrate the validity and robustness of the proposed approach with respect to the existing control methods.

The paper is organized as follows: Section 2 introduces the dynamical model of the system. Section 3 describes the control design. Section 4 presents a case study and simulation results. Finally, conclusions and future work are drawn in Section 5.

Q 2 Q n z x y 2 2 2 z x n n Q 1 z x y 1 1 1 O X Y Z q n q 2 y n q 1
x des L Fig. 1: Aerial System: general architecture of the aerial system considered.

DYNAMIC MODEL

In this work, the tracking problem was addressed by considering the system as composed of a team of n quadrotors and n massless rigid links that collaborate to carry the point-mass load along the desired trajectory as illustrated in Figure 1. The equations of motion were obtained by using the energetic method [START_REF] Lee | Geometric control of cooperating multiple quadrotor uavs with a suspended payload[END_REF][START_REF] Lee | Geometric control of quadrotor uavs transporting a cable-suspended rigid body[END_REF][START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF][START_REF] Tang | Aggressive maneuvering of a quadrotor with a cablesuspended payload[END_REF]. Throughout all this paper the used symbols were summarized in the Nomenclature Section.

Since each quadrotors lies on the surface of a sphere of radius l i i.e. length of the i th link, the configuration manifold associated with the mentioned system is R 3 ×(S 2 ×SO(3)) n where S 2 = {q ∈ R 3 | ||q|| = 1} represents the 2-sphere set whereas SO(3) = {R ∈ R 3×3 | R T R = I, det(R) = 1} defines the Special Orthogonal Group. Recalling that the action integral for the non conservative forces is

a = t f t 0 (T -U + W f + W m ) dt, (1)
where T is the kinetic energy, U is the potential energy and W f and W m represent the works of external forces and moment, respectively. Now, Hamilton's principle of least action states that the path a conservative system takes from configurations Ψ(t 0 ) to Ψ(t f ), where Ψ is the set of generalized coordinates belonging to the configuration manifold, is the one that makes the action stationary. In mathematical terms

δa = t f t 0 (δT -δU + δW f + δW m ) dt = 0, (2)
where δ expresses the variation in the integrated quantities. Hence, integrating all these terms by parts 3 returns the equations of motion

Q i z x y i i i v i v i v i
                               d dt x L = ẋL , (3)
( n i=1 m i q i q T i + m L I)(ẍ L + ge 3 ) = n i=1 v i -m i l i ||ω i || 2 q i + (q i • f i,w )q i + f L,w , (4)
qi = ω i × q i , (5) ωi
= 1 l i qi (ẍ L + ge 3 ) - 1 m i l i qi v ⊥ i - 1 m i l i qi f i,w , (6)
Ṙi = R i Ωi , (7)
m i = J i Ω + Ω i × J i Ω i , (8)
where f i,w and f L,w are the aerodynamic forces acting on the quadrotors and the load, respectively, I is the identity matrix of dimension three, the vector v i represents the control force of the i th robot such that v i = f i R i e 3 whereas the vectors v i and v ⊥ i denote the projection of v i along and normal to q i , respectively

v i =(I + q2 i )v i = (q i • v i )q i = q i q T i v i , (9a)
v ⊥ i = -q2 i v i = -q i × (q i × v i ) = (I -q i qT i )v i . (9b)
The symbol • denotes the hat map:• : R 3 → so(3) where so(3) corresponds to the Lie algebra associated to the Lie group SO [START_REF] Jean | Analysis and control of nonlinear systems: A flatness-based approach[END_REF]. The relationship between Eqs. (9a) and (9b) is depicted in Figure 2 and it can be expressed as

v i = v ⊥ i + v i . ( 10
)
Remark 2.1. The need to distinguish between parallel and orthogonal components comes from the structure of the equations of motion. In fact, Equation (4) governs the translational dynamics of the load and then it contains the parallel component of the control force. On the other hand, Equation (6) governs the attitude of the load and therefore includes the orthogonal component of the control force. Those properties play a central role in designing the control.

Remark 2.2. Compared to previous works, Equations (3)-( 8) present the introduction of the wind forces acting on the quadrotors and the payload.

CONTROL

The overall objective of the control is to ensure that the load follows the desired trajectory while pursuing the desired system configuration under the action of external disturbances. The relative formation control [START_REF] Lee | Geometric control of cooperating multiple quadrotor uavs with a suspended payload[END_REF] is used. Thus, the objective is to find a set of control inputs U = {f i , m i ∀ i = 1, . . . , n} such that x L -→ x des L and q i -→ q des i while t -→ ∞. Furthermore, to ensure insensitivity to bounded external disturbances, a SM-APID controller is designed to define a proper action that steers the load. The fusion of the two above-mentioned approaches enables to control ACTSs even in presence of external disturbances.

In the first Subsection 3.1 the structure of the SM-APID controller and aspects concerning its design are recalled. Subsection 3.2 concerns with the contribution of this work since it shows how to exploit the SM-APID control to guide the position of the load subject to external disturbances. In Subsection 3.3, the design of the orthogonal component to reach the desired configuration is done and in Subsection 3.4 the attitude control for the quadrotors is recalled for sake of completeness.

Design of the Sliding Mode-Adaptive PID control

Let us consider the description of a dynamical system in general form

ẋ = f (x) + ∆f (x) + g(x)u + δ, (11)
where x, ẋ ∈ R n form the system' state

x = [x, ẋ] T ∈ R 2n , M x , f (x
) and ∆f (x) represent the system model and the its uncertainties, g(x)u is the affine terms with u as control input and, the last terms δ identifies the external disturbances. Observe that both disturbances and system uncertainties are assumed to be bounded above

δ(•) ≤ δ u (•), f (•) ≤ f u (•) and ∆f (•) ≤ ∆f u (•).
A common and direct approach that can be used to define a control input consists in using the so-called feedback linearization techniques namely

u * = -g(x) -1 f (x) + ∆f (x) + δ + ẍdes + k 0 ẽ + k 1 ė . ( 12
)
where e = [ẽ, ė] T ∈ R 2n and ė = [ ė, ë] T ∈ R 2n are the state error vectors which are the time derivatives of e = yx with y = xdes , ẋdes representing the desired state vector. Despite its simplicity, this control action is not robust technique because of the linerization it requires for guaranteeing the stability. Hence, this aspect constitutes another way to stress the necessity of introducing a robust technique.

Remark 3.1. Note that u * depends on both f (x) and ∆f (x). This will be useful in what follows to understand the condition of insensitivity of the proposed methodology to them.

Control Architecture

Generally, when dealing with adaptive control, the control input is split in two terms [START_REF] Hsueh | Supervisory controller design based on lyapunov stable theory[END_REF][START_REF] Yuri | Sliding Mode Control and Observation[END_REF]

u = u P ID + u s , (13)
where u P ID covers the role of the so-called equivalent control u eq whereas u s is the supervisory control. This latter is an additional control that keeps the system state within some defined boundaries and guarantees the stability of the dynamical system while the u P ID exploits a gradient-based adaptation law for updating its gains providing robustness of the control. The design of both u s and u P ID that is intended to be used here has been defined in [START_REF] Chang | Adaptive robust pid controller design based on a sliding mode for uncertain chaotic systems[END_REF] and therefore only the key concepts are reported in the followings.

Design of the Supervisory Control

The primary role of the supervisory control u s is to keep the system state inside a designed constraint set

C = {x ∈ R 2n | ||x|| ≤ M x }, (14)
where M x is a pre-specified parameter usually chosen such that M x ≥ ||y|| ∞ . Its design relies on the errors dynamics and the Lyapunov criterion to guarantee the stability of the system. Hence, to asymptotically attain the zero-error condition, the following Lyapunov function candidate is considered

V e = 1 2 e T Φe. (15)
Observe that V e represents an output error energy measuring the distance from a representative point in the state space to the desired output. Therefore, a plausible strategy to reach it is to exercise a control action that strictly decreases V e . Before continuing this analysis, it is better to recall the role of the matrix Φ. Generally, it is defined as a positive definite and symmetric matrix that springs out as the solution of the Lyapunov equation [20]

A T Φ + ΦA = -Q, (16)
where matrix Q ∈ R 2n×2n is a given, positive definite symmetric matrix whereas matrix A ∈ R 2n×2n pops out by rewriting the error dynamic as

ė = Ae + B(u * -u P ID -u s ) (17)
where the structure of both A and B is

A = 0 I -k 0 -k 1 and B = 0 g(x, ẋ) , (18)
note that k 0 ∈ R n×n and k 1 ∈ R n×n are diagonal matrices whose coefficients are chosen such that the roots of the characteristic polynomial, associated with the differential equation ë + k 1 ė + k 0 ẽ = 0, belong to the left-half complex plane. This aspect is fundamental to guarantee the stability of the system and the convergence to zero of the error as time goes to infinity. Hence, the design of u s must comply with the decrease of V e , which means that it must meet

Ve < 0. ( 19
)
The derivative of the Lyapunov function can be computed as follows

Ve = 1 2 ( ėT Φe + e T Φ ė) = 1 2 (e T (A T Φ + ΦA)e + 2e T ΦB(u * -u P ID -u s )) = - 1 2 e T Qe + e T ΦB(u * -u P ID -u s ) ≤ - 1 2 e T Qe + |e T ΦB|(|u * | + |u P ID |) -e T ΦBu s . (20)
Therefore, to satisfy, the supervisory controller can be chosen as

u s = sgn(e T ΦB)(|u * | + |u P ID |). (21)
However, with this design of u s , it constantly intervenes in the control process. Further, the presence of sgn function leads to chattering. Therefore, to adhere with its supervisory definition and to reduce the chattering the set of constraints C is used. In particular, the Indicator function I f is introduced in the u s design as follows

ũs = I f u s where I f = 1, V e > V M , 0, V e ≤ V M , (22)
and

V M = 1 2 λ min (Φ)(M x -||y|| ∞ ) 2 , (23)
where λ min (Φ) is the minimum eigenvalue of the Lyapunov matrix Φ [START_REF] Chang | Adaptive robust pid controller design based on a sliding mode for uncertain chaotic systems[END_REF][START_REF] Zabczyk | Mathematical Control Theory[END_REF]. The introduction of the Indicator function completes its design. Indeed, thanks to this modification, ũs would intervene in the control process only when the system's state diverges from the desired over the limits imposed by M x .

Remark 3.2. Is worth observing that the value of u * inside Eq.( 21) is bounded above. Indeed, rearranging Equation (12) it becomes possible finding an estimate for u * of the type

u * ,u = g -1 δ u + f u + ∆f u + | ẍdes | + |k 0 ẽ| + |k 1 ė| .

PID Adaptive Laws

The main control activity is then delegated to an Adaptive PID (APID) which adjusts its control action by automatically tuning its gains. The adaptation laws are derived with the aim to reach the so-called sliding mode S = 0 (i.e. insensitivity to external disturbances) where S is the so-called sliding surface defined as [START_REF] Yuri | Sliding Mode Control and Observation[END_REF]

S = ẋ -x r , (24)
where the reference signal x r is

x r = ẋdes + k 1 ẽ + k 0 ẽ dt. ( 25
)
Remark 3.3. When the sliding mode occurs one has

S = 0 and ẋ = x r , (26)
this, substituted in the derivative of x r , ẋr = ẍdes + k 1 ė + k 0 ẽ, returns ë + k 1 ė + k 0 ẽ = 0 thus revealing that, when the system is in the sliding mode, it becomes insensitive to external disturbances. On the other hand, when S = 0, ë + k 1 ė + k 0 ẽ = g(x, ẋ)(u *u P IDu s ) and therefore the system is not fully insensitive to external disturbances since they appears inside u * .

To guarantee approaching the sliding mode the Lyapunov function approach is exploited and the Lyapunov function candidate is then chosen as

V = 1 2 S 2 . ( 27
)
Hence requiring that S(t) → 0 for t → ∞ coincide with reducing V . Consequently, the gradient method is applied to choose the gains which take directions of maximum slope over V . Now, a common expression for a PID controller that can compensate for the gravity acceleration is

u P ID = K P ẽ + K I ẽ dt + K D ė + ge 3 , (28)
using the gradient method and the chain rule, it is possible to obtain the adaptation laws for the control gain matrices K P , K I and

K D KP,ii = -γ ∂S Ṡ ∂K P = -γ ∂S Ṡ ∂u P ID ∂u P ID ∂K P = -γSẽ , (29a)
KI,ii = -γ ∂S Ṡ ∂K I = -γ ∂S Ṡ ∂u P ID ∂u P ID ∂K I = -γS ẽ dt, (29b)
KD,ii = -γ ∂S Ṡ ∂K D = -γ ∂S Ṡ ∂u P ID ∂u P ID ∂K D = -γS ė , ( 29c
)
where minus is placed opposite to the energy flow V and γ ∈ R + is called learning rate1 .

To conclude, the behaviour of the controller can be resumed as follows: if u s = 0, the PID gains adapt themselves to decrease V to zero (i.e. reach the sliding mode) whereas if u s = 0, the PID gains are not able to decrease V and then also V e under V M .

Design of the Parallel Components

The parallel component is chosen to track the desired trajectory of the load, in other words it is necessary to define v i s.t. x L -→ x des L while t -→ ∞. Let's consider the load Equation (4). Since the translational dynamic is influenced only by the terms parallel to q i then, it is possible to derive v i as

v i = m i l i ||ω i || 2 q i + t i q i + m i m L q i q T i n j=1 t j q j -(q i • f i,w )q i , (30)
where t i represent the tension along the i th link. To find the tensions to be applied along the links, it is enough to substitute (30) into ( 4) and make some arrangements to see that

(ẍ L + ge 3 ) = 1 m L n j=1 t j q j + n i=1 m i q i q T i + m L I -1 f L,w, (31)
and considering

n i=1 m i m L q i q T i + I -1
as a dimensionless distribution mass-matrix M, one has

(ẍ L + ge 3 ) = 1 m L n j=1 t j q j + 1 m L M f L,w , (32)
which leads to

ẍL = -ge 3 + 1 m L w + 1 m L δ(t), (33)
where δ(t) represents the effect of the external disturbance whereas w is the wrench that act on the load. Rearranging this equation in the form of Newton' law, one has

w = m L (ẍ L + ge 3 ) - 1 m L δ(t) , (34)
grouping the terms at RHS inside u, returns

w des = m L u, (35)
where the superscript des in w des emphasizes that the control action coincides, up to the m L constant, with the wrench that should be generated on the load to track the desired trajectory. In other words, it represents the action necessary to compensate for the external forces, disturbances and then correct the tracking error. As a consequence, the tensions necessary to steer the load, for n ≥ 3, can be found as

t = W † w des + Nλ, (36)
where W = [q 1 , q 2 , ..., q n ] ∈ R 3×n is the so-called wrench matrix, W † = W T (WW T ) -1 is its pseudoinverse, N ∈ R (n-3)×n contains the vector that span the kernel of W and λ ∈ R (n-3) .

Design of the Orthogonal Components

The normal component is chosen to reach the desired relative configuration of the system, in other words it becomes necessary to define v ⊥ i s.t. q i -→ q des i while t -→ ∞. This correspond to the tracking problem on S 2 which was studied in [START_REF] Bullo | Geometric control of mechanical systems: modeling, analysis, and design for simple mechanical control systems Springer[END_REF][START_REF] Wu | Variation-based linearization of nonlinear systems evolving on so(3) and S 2[END_REF]. In particular, the direction and angular velocity errors, are defined as follows e q i = q des i × q i and e ω i = ω i + q2

i ω des i ,

whereas, the angular acceleration is

ωi = -K q i e q i -K ω i e ω i -(q i • ω des i ) qi -q2 i ωdes i , (38)
for positive gains K q i and K ω i . Now, substituting (31) into ( 6) yields to

ωi = 1 m L l i qi n j=1 t j q j - 1 m i l i qi v ⊥ i - 1 m i l i qi f i,w , (39)
and v ⊥ i = -q 2 i v i leads us to the expression for the normal component

v ⊥ i [4] v ⊥ i = m L l i qi (-K q i e q i -K ω i e ω i -(q i • ω des i ) qi -q2 i ωdes i ) - m i m L q2 i n j=1, j =i t j q j + q2 i f i,w . (40)

Quadrotor's Attitude Control

To fully control the overall system, the attitude control of each quadrotor must be included. Hence, as done in [START_REF] Lee | Geometric tracking control of a quadrotor uav on se(3)[END_REF], the control moment for each quadrotor is

m i = -K R i e R i -K Ω i e Ω i + Ω i × J i Ω i -J i ( Ωi R T i R c i Ω c i -R T i R c i Ωc i ), (41)
where • c stands for controlled variables, • ∨ denotes the vee map:

• ∨ : so(3) → R 3 , K R i , K Ω i are diagonal
positive-gain matrices and e R i , e Ω i are the errors for the attitude dynamics defined as

e R i = 1 2 (R T c i R i -R T i R c i ) ∨ and e Ω i = Ω i -R T i R c i Ω c i , (42)
Equation ( 41) together with

f i = v i • R i e 3 , (43)
constitute the control actions for each agent of the system. The overall control scheme is depicted in Figure 3. Fig. 3: Control scheme: it illustrates the various blocks that define the control logic and the quantities that are exchanged while the system is working.

SIMULATIONS

To shed some light on the peculiarities of the proposed approach numerical simulations are made. In particular, the tracking performances of different control are compared. An aerial system, with three rigid links and three quadrotors transporting a point-mass load, is considered. The task to be completed involves guiding the load along an elliptical trajectory. Both the load and the quadrotors were subjected to a wind that persisted for the entire simulation. To compute the aerodynamic forces f L,w and f i,w empirical equations are used for the former and some experimental measures for the latter, respectively. In particular,

f L,w = 1 2 C d A L ρ(v w -|| ẋL (t)||) 2 e
∞ where e ∞ represent the direction of the relative velocity (i.e. unit vector) and whose parameters were reported in the Table 1 together with other useful parameters for the simulation. The value of f i,w is difficult to derive empirically since the shape of quadrotors is complex and not unified. Consequently, its value was estimated through experimental tests carried out at the LS2N arena [START_REF] Erskine | Wrench analysis of cable-suspended parallel robots actuated by quadrotor unmanned aerial vehicles[END_REF]. In particular, thanks to the use of a fan (shown in Figure 4a) and a first-order observer on a quadrotor, it was possible to determine a link between wind force and wind speed. This, together with the Dryden's model, allowed us to derive the wind forces acting on the drones as shown in Figure 4b. All this continues to be an approximation of reality but their proper modelling is beyond the scope of this paper. What is of interest instead, is to show that the presence of a persistent and time-varying force of aerodynamic nature, can be managed by the proposed control.

Parameter Value

Parameter Value

A x 0.8 m m 1 = m 2 = m 3 1.15 kg A y 0.5 m l 1 = l 2 = l 3 1 m m L 0.5 kg k 0 diag 0.1 0.1 0.1 D 0.2 m k 1 diag 1 1 1 C d 0.5 J 1 = J 2 = J 3 diag 0.01 0.1 0.14 kg • m 2 ρ 1.2 kg m 3 K P (0) = K I (0) = K D (0) diag 5 5 5 A L πD 2 4 K P = K I = K D diag 5 5 5 γ 10 K q i = K ω i = K R i = K Ω i diag 5 5 5 T 20s Ω i (0) = ẋL (0) = ω i (0) = ω des i 0 0 0 R i (0) I x des L A x cos(ωt) A y sin(ωt) 0 m i des i e 1 x L (0) A x 0 0.5 m q 1 (0) sin 80 • 0 -cos 80 • q des 1 0 -sin 40 • -cos 40 • q 2 (0) 0 sin 80 • -cos 80 • q des 2 sin 40 • 0 -cos 40 • q 3 (0) -sin 80 • 0 -cos 80 • q des 3 0 sin 40 • -cos 40 •
The initial configuration of the system and some phases of the task are shown in the following figure with the purpose of showing the overall behaviour of the system during the execution of the task Figure 5. It is noteworthy that the system properly reconfigures itself to reach the desired output while counteracting the external disturbances. In addition, the sliding variables trends and the gains curves were reported in Figures 7 and6, respectively. Figure 6 shows that k P and k I gains stabilize after a transitory phase that last, approximately, 5 seconds whereas k D gains continue to adjust in order to reach the sliding surface. Indeed, from Figure 7, one can see that the sliding variable values first decrease and then oscillate about zero for the entire duration of the simulation. The attainment of the sliding condition in finite time is not guaranteed with this design. Indeed, in Section 3.1.3, the main requirement states as S → 0 for t → ∞. It should also be noted that, since this condition is not reached in the first 10 seconds, when the wind changes intensity, Figure 4b, the sliding variables locally increase their value and then decrease again. Thus, although the system may never achieve effective insensitivity from external disturbances, it does reduce their effect on the system by making the controller more robust. This can be observed with the following comparison. Comparisons made by the controller developed in [START_REF] Lee | Geometric control of cooperating multiple quadrotor uavs with a suspended payload[END_REF] (a PD control) and a revised version of it [START_REF] Li | Cooperative transportation of cable suspended payloads with mavs using monocular vision and inertial sensing[END_REF] (PID). To qualitatively compare their performance, reference can be made to Figures 8 and9. The first two images show the trajectory followed by the load during the simulation with different control techniques. The third image shows the error trends along with the three directions (x, y and z). Therefore, from these pictures, it can be inferred that both PD and PID controllers fail to accomplish the task. In fact, the PD only manages to keep the error within a certain value. Instead the PID, thanks to the integral part, seems to recover the trajectory. However, during the recovery phase, the quadrotors perform a sudden manoeuvre which could not be achievable in reality since it could lead the system to crash. Indeed, the gains of these controllers (PD and PID) are chosen sufficiently high in order to guarantee the performance shown here and then to make comparisons. In general, with lower gain values both the controllers are not able to conclude the task. On the other hand, this result confirms the robustness of the proposed controller.

Additional Figures 10 and 11 also show the error trend of the attitude of each quadrotor and the error of the configuration (link directions) as well. It should be noted that the attitude and configuration errors are quite similar for both the PD and PID. In contrast, in the case of the SM-APID, it can be seen that the e R SM-APID PD q 1 q 7 q 8 q 1 q 7 q 8 q 1 q 7 q 8 q 1 q 7 q 8 q 1 q 7 q 8 q 1 q 7 q 8 Fig. 11: Configuration errors: the link direction error for each bar in the three study cases is reported. The error equation on S 2 is: e q = 1q i • q des i .

"#$ PD R 1 R % R Q R 1 R % R Q 0 5 10 15
H H H H! 0 e 1 R 1 R % R Q
converges slightly faster at zero while e q error remains lower than PD and PID controllers. However, in all cases e q is not zero because the controller does not compensate for the external disturbances. The reason why it is lower for the proposed controller is due to the presence of the term q2 i f i,w in Equation (40). This is absent in the other two controllers.

On the other hand, a tangible way to compare the performances is provided by Table 2. It makes a statistical analysis of the errors. In this case, Table 2 refers to the position error only since, in this example, a tracking task is considered. To conclude, Table 2 qualitatively confirms that, in general, the errors (MEAN, STD or RMSE) along each direction are smaller in the case of the SM-APID controller.

CONCLUSION

The analogy between chaotic systems and wind considered as such enabled to make current controllers more robust in presence of external disturbances. In particular, the design of a controller to guide an aerial system subjected to wind gusts was tackled. Both an online tuning gains PID and a supervisory control were integrated with existing techniques to steer the system immersed in a wind field. Furthermore, simulations showed the validity of our approach thus confirming the validity of this analogy.

Despite the presence of many existing techniques for the estimation and compensation of aerodynamic forces, this is the first time that a controller of this type has been used in this context. Futures effort aims to develop a model that takes into account the flexibility of the cables and in doing experimental validation of our approach.

Fig. 2 :

 2 Fig. 2: Control force decomposition: representation of the i th projection of the control force along a generic link.

  Wind force estimation: quadrotor windswept during the experimental estimation of the aerodynamic force. Wind forces: wind force components acting on the aerial system.

Fig. 4 :

 4 Fig. 4: Wind force estimation and profiles.

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig. 5: Simulation phases: snapshots of the system during the simulation with SM-APID control. The snapshots show the system in (a) the initial condition t = 0s, (b) after one turn along the ellipses were the wind suddenly disappears t = 11s and (c) in the final phase of the task t = 19s. The red arrows, where present, show the presence and direction of the aerodynamic force acting on the quadrotor and load.

Fig. 8 :

 8 Fig. 8: Trajectories comparison: spatial and planar view of the undertaken load trajectory.

Fig. 9 :

 9 Fig. 9: Position errors: the components of the position errors for the three study cases are hereby depicted.

Fig. 10 :

 10 Fig. 10: Attitude errors: orientation error for each quadrotors in the three study cases is given. The error equation on SO(3) is: e R = 1 2 trace(I -R T c i R i ).

Table 1 :

 1 Simulation parameters.

Table 2 :

 2 Error analysis: the Mean, STD and RMSE values for the trajectory tracking problems is computed for the three study cases.

	Trajectory Control Type Statistics	Load position (m)
				x	y	z
			MEAN	0.137 0.173 0.008
	Ellipses	PD	STD	0.124 0.163 0.094
			RMSE	0.185 0.238 0.094
			MEAN	1e -4 2e -4 1e -4
	Ellipses	PID	STD	0.093 0.107 0.103
			RMSE	0.093 0.107 0.103
			MEAN	0.002 0.001 3e -4
	Ellipses	SM-APID	STD	0.049 0.034 0.090
			RMSE	0.049 0.034 0.090

To avoid cumbersome notation the subscripts ii (emphasizing the diagonal structure of the gain matrices) has been inserted only to the left hand side of the equation intending the derivative operation to be carried out component-wise. As typical in the literature[START_REF] Noordin | Adaptive pid controller using sliding mode control approaches for quadrotor uav attitude and position stabilization[END_REF].


Thomas Duigou

Melchior Joan Hérisson

Du Lac

Mahnaz Kenza Bazi-Kabbaj

Sabeti Azad

Gizem Buldum

Olivier Telle

Yorgo El-Moubayed

Pablo Carbonell

Neil Swainston

Mahnaz Sabeti

Galaxy-SynBioCAD: tools and automated pipelines for Synthetic Biology Design and Metabolic Engineering

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Centre

Île-de-France -Jouy-en-Josas -Antony Domaine de Vilvert 78 352 Jouy-en-Josas Cedex contact : thomas.duigou@inrae.fr www.inrae.fr

Galaxy-SynBioCAD: tools and automated pipelines for Synthetic Biology Design and Metabolic Engineering


Bozhao Wang

Philippe Cardou
email: philippe.cardou@gmc.ulaval.castephane.caro@ls2n.fr

Stéphane Caro

An Approach for Predicting the Calibration Accuracy in Planar Cable-Driven Parallel Robots and Experiment Validation

HAL is

Introduction

Cable-driven parallel robots (CDPRs) are a group of parallel robots that are actuated through flexible cables instead of rigid links. Compared with serial robots, this particular type of robots benefits from its high flexibility [START_REF] An | An All-in-one Cable-driven Parallel Robot with Flexible Workspace and Its Auto-calibration Method[END_REF], high payload to weight ratio [START_REF] Picard | Control Strategies for a Cable-Driven Parallel Robot with Varying Payload Information[END_REF], reconfigurability [START_REF] Gagliardini | Discrete reconfiguration planning for Cable-Driven Parallel Robots[END_REF], and potentially large translational workspace [START_REF] Li | Optimal Force Distribution Based on Slack Rope Model in the Incompletely Constrained Cable-Driven Parallel Mechanism of FAST Telescope[END_REF]. On most CDPRs, the moving-platform (MP) is usually connected to several cables, then through pulleys, winches and then linked to the base structure. There are two configurations according to the number and spatial position of cables used: suspended and fully-constrained CDPRs [START_REF] Picard | Pulleys and Force Sensors Influence on Payload Estimation of Cable-Driven Parallel Robots[END_REF].

Parallel mechanisms, because of their large numbers of links and passive joints [START_REF] Wang | On the accuracy of a Stewart platform. I. The effect of manufacturing tolerances[END_REF], may not necessarily have a high accuracy. Therefore, kinematic calibration is important for such robot architectures. Previous works have implemented plenty of robot calibration methods, with nonlinear least square method (NLLS) being the most common one. However few works focused on the calibration with the combination of several sensors. Besides, few existing studies aimed at identifying the Cartesian coordinates of the exit points, the initial cable lengths and the movingplatform poses. For CDPRs with incremental encoders, the initial cable length, therefore the initial platform pose is unknown, and is typically different each time for applications. Identifying both the initial cable lengths and the Cartesian coordinates of the cable exit points helps have a better knowledge of the initial platform pose. Daney et al. [START_REF] Daney | Interval methods for certification of the kinematic calibration of parallel robots[END_REF] used a method based on interval arithmetics for the kinematic calibration of parallel robots. In [START_REF] Zhang | Calibration of geometric parameters and error compensation of non-geometric parameters for cable-driven parallel robots[END_REF], a high precision and robustness iterative calibration method is proposed to significantly improve the end-effector position errors. Klimchik et al. [START_REF] Klimchik | Design of Calibration Experiments for Identification of Manipulator Elastostatic Parameters[END_REF] used a calibration technique to compensate the elastic deflections of the manipulator components. The method proposed in [START_REF] Verner | Optimal calibration of parallel kinematic machines[END_REF] takes into account of geometric errors and exploits the least error sensitive regions to perform optimal calibration. Wu et al. [START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF] developed an irreducible geometric model to reduce measurement noises during calibration and proposed an approach to find the optimal robot calibration configurations. For CDPRs, most of previous works relied on non-linear least square (NLLS) methods for parameter identification, as it applies conveniently to the minimization of the cable length residuals. In [START_REF] Qian | Kinematic calibration of a cable-driven parallel robot for 3D printing[END_REF], the NLLS method is used on a 6-cable, 3-degree of freedom (DoF) CDPR, with a parallelogram, which is able to reach a larger workspace. The application is proven to be effective by simulations and experiments. The authors also proposed an algorithm to select optimal measurement poses. The authors of [START_REF] Zhang | Kinematic Calibration of Cable-Driven Parallel Robots Considering the Pulley Kinematics[END_REF] used NLLS method for a CDPR called TBot and also considered pulley kinematics. They proposed a measurement pose optimization method, which consists in minimizing the identification matrix condition number. In [START_REF] Sandretto | Calibration of a Fully-Constrained Parallel Cable-Driven Robot[END_REF], several identification methods derived from NLLS are proposed and tested. The other similar approach is orthogonal distance regression (ODR) [START_REF] Boggs | A Stable and Efficient Algorithm for Nonlinear Orthogonal Distance Regression[END_REF].

Sensors used for robot calibration can be divided into two categories: exteroceptive and proprioceptive, which can be used in combination or not. Proprioceptive motor position sensors are used in [START_REF] Miermeister | Auto-calibration method for overconstrained cabledriven parallel robots[END_REF] for auto-calibration. Borgstrom et al. [START_REF] Borgstrom | NIMS-PL : A Cable-Driven Robot With Self-Calibration Capabilities[END_REF] used proprioceptive tension sensors and cable encoders to achieve two novel jitter-based and tension-based self-calibration methods. Zhang et al. [START_REF] Zhang | Calibration of geometric parameters and error compensation of non-geometric parameters for cable-driven parallel robots[END_REF] relied exclusively on external sensors to perform their iterative calibration method. A single theodolite is used in [START_REF] Zhuang | Kinematic calibration of a Stewart platform using pose measurements obtained by a single theodolite[END_REF] for kinematic calibration of a parallel mechanism. Martin et al. [START_REF] Martin | A laser-based direct cable length measurement sensor for CDPRS[END_REF] used laser-based cable length measurement sensor to improve calibration quality. In [START_REF] Zhang | Kinematic Calibration of Cable-Driven Parallel Robots Considering the Pulley Kinematics[END_REF], an external laser tracker is used to obtain accurate robot poses as the input of the calibration problem. A efficient calibration method concerning partial pose measurements, without end-effector orientation is proposed in [START_REF] Klimchik | Geometric and elastostatic calibration of robotic manipulator using partial pose measurements[END_REF]. Because of the high accuracy of laser tracker devices, they are also often used to provide ground truths of unidentified parameters.

Little research has been done on the combination of different types of sensors in CDPRs.

Calibration with exteroceptive sensors has drawbacks such as poor measurement accuracies over volumes, and the difficulty to identify the end-effector pose in certain cases [START_REF] Renaud | Kinematic calibration of parallel mechanisms: A novel approach using legs observation[END_REF]. The use of sensors embedded on the MP can mitigate those disadvantages. The calibration in this study integrates a laser displacement sensor, an inclinometer that are installed on the MP, and the motor encoders. To the best of the authors' knowledge, all those proprioceptive sensors have not been used altogether before for CDPR calibration. This work is expected to refine the sensor selection process for a prescribed CDPR calibration accuracy. Redundant internal joint encoders are used for the self-calibration of a parallel mechanism [START_REF] Zhuang | Autonomous calibration of hexapod machine tools[END_REF]. Renaud et al. [START_REF] Renaud | Kinematic calibration of parallel mechanisms: A novel approach using legs observation[END_REF] proposed to perform the kinematic calibration of a parallel mechanism by observing its legs with a camera. The method combines the advantages of both exteroceptive and proprioceptive sensors. Andreff et al. [START_REF] Andreff | Vision-based kinematic calibration of an H4 parallel mechanism: Practical accuracies[END_REF] proposed a kinematic calibration method with computer vision for H4 mechanism, a four-chain parallel mechanism that offers three translational and one rotational DoFs [START_REF] Pierrot | H4 parallel robot: modeling, design and preliminary experiments[END_REF]. In ref. [START_REF] Daney | Interval methods for certification of the kinematic calibration of parallel robots[END_REF], Daney et al. used internal sensors to provide the leg length differences of a parallel mechanism.

This paper presents a method of simulating and predicting the calibration accuracy of CD-PRs. From the 2-cable, 3-DoF planar CDPR under study, a calibration method is proposed in this work, which concerns the combination of a laser displacement sensor, an inclinometer attached to the moving-platform, as well as the motor encoders. The proposed method aims to identify the Cartesian coordinates of the cable exit points, the initial cable lengths and the moving-platform poses. The calibration task is formulated as an optimization problem, the calibration method uses an NLLS algorithm to minimize the cable length residuals. The calibration process is first simulated for acquiring the expected system parameter error ranges, and then verified by experiments. Other than the specific CDPR under study, this method of simulation and calibration accuracy prediction model can be used on other CDPR configurations. The CDPR calibration accuracy is directly affected by the sensor measurement errors and polymer cable elongations. During the simulation of the calibration process, the sensor measurement errors are modelled and compensated, and the effect of each sensor on the calibration quality is studied. Based on the authors' previous work [START_REF] Wang | An Approach for Predicting the Calibration Accuracy in Planar Cable-Driven Parallel Robots[END_REF], this study implements the experimental validation process, where a laser tracker is used to obtain the system variable ground truths to eliminate the effects of cable elongations and moving-platform (MP) manufacturing errors. In addition, the CDPR geometric model considers the pulley geometry, random MP poses and the floor elevation. Thus, the calibration quality is improved. As a result, the calibrated system parameter errors fall within the simulated ranges.

The rest of the paper is organized as follows: Section 2 focuses on the CDPR geometric modelling and on the accuracy tests that contribute to calibration quality. The robot position control scheme is also mentioned. Section 3 describes the identification methodology of the current study. The simulation results are discussed in Section. 4. Section 5 presents the calibration experiment methodology in detail. Then section 6 presents an analysis of the results in this experiment. Finally, some conclusions are drawn.

CDPR modelling and control

The CDPR under study operates on the CRAFT CDPR prototype located in LS2N, Nantes, France, measuring 3.8 m × 4.3 m × 2.8 m. Figure 1 shows the real 3-DoF moving-platform and the detailed MP structure. The MP uses a Bosch Rexroth 40×40 strut with a effective length of 500 mm. It is suspended by the two cables, which are fixed on the MP with angle brackets and then go through two pulleys in the diagonal direction of the prototype. In such a way, the planar CDPR workspace is formed. The related cable exit points next to the pulleys are noted as A 1 and A 2 , respectively. Pins and supports are designed to hold the MP steadily on the stand. The CDPR base frame F 1 and workspace frame F 2 are accurately defined with a laser tracker, with a 5.2 m by 2.8 m workspace frame size. The origin O 2 of F 2 is the projection of exit point A 1 on π 1 , the horizontal plane of F 1 . The x 2 axis is horizontal and the z 1 axis is vertical. Finally the two cables are led to the actuated winches fixed on the corresponding bottom corners in the base frame. The horizontal movement of the winch exit point because of cable winding is neglected, and the cables are assumed to be massless and straight.

The CDPR schematic is shown in Fig. 2. The MP has one rotational and two translational DoF in the planar workspace. However it is held by only two cables, which makes the robot under-constrained. With given cable lengths, the MP still has one degree of freedom to move, but remains at the pose where its gravitational potential energy is the smallest. For calibration, the MP is equipped with an inclinometer and a laser displacement sensor connected to the the MP bar through a revolute joint, as shown in Fig. 2. Thanks to this revolute joint, the laser sensor always points vertically to the floor, to directly measure the MP height. The concrete floor reflects enough laser light to perform proper measurements. It is not perfectly flat, however, so a laser tracker is used to identify its unevenness. Motors control the cable lengths, and are equipped with encoders that measure their angular positions, and thus the cable length variations. The total MP mass is 2.5 Kg.

The geometric model of the current CDPR, including detailed pulley modelling is presented in Sec. 2.1. The detailed random pose generation process is discussed in Sec. 2.2. The sensor accuracy tests are carried out and the results are summarized in Sec. 2.3. The floor elevation along the workspace direction is accurately measured with the laser tracker and the corresponding results are detailed in Sec. 2.4. Finally the position control scheme of the CDPR is presented in Sec. 2.5.

Geometric model

Figure 3 shows the ith loop of the current CDPR. The vectors pointing from the workspace origin O 2 to the ith pulley exit point and to the MP center are a i and p j , respectively; with i = 1, . . . , m and j = 1, . . . , n. m is the number of cables (2 in our case), and n is the total number of measurement poses. The vector pointing from P to the ith anchor point is b i . Therefore, the loop closure equation corresponding to each cable is expressed as:

l i = l i u i = --→ A i B i = p + b i -a i (1)
where l i is the i-th cable vector and l i is the i-th cable length. Then the i-th unit cable vector is written as:

u i = l i l i (2)
Detailed pulley modelling is considered in the geometric model, which results in the shift of cable exit point from A i to A i [START_REF] Wang | Exit Point, Initial Length and Pose Self-calibration Method for Cable-Driven Parallel Robots[END_REF][START_REF] Picard | Pulleys and Force Sensors Influence on Payload Estimation of Cable-Driven Parallel Robots[END_REF], shown in Fig. 4. The frame attached to the pulley is denoted as F pu = (A i , x pu , y pu , z pu ). The axis x pu is horizontal and goes through the pulley center O pi , the axis z pu is vertical and the axis y pu is parallel to the pulley rotating axis. The vector from the pulley exit point A i to the cable anchor point B i is l ci . The jth MP rotational angle compared with the horizontal plane is θ j . As a result, the cable length from the cable exit point A i to cable anchor point B i is:

l ti = l pi + l ci (3)
where l ci denotes the cable length from exit point A i to B i , l pi is the cable length wrapped on

l pi = r p π -β i -γ i ( 4
)
where r p is the pulley radius, with tan(

β i ) = √ m i m T i -r 2 p rp and sin(γ i ) = a iz -b iz ||m i || 2
. m i is denoted as the vector pointing from the pulley center O pi to the anchor point B i , and is expressed as:

m i = p + b i -a i + r p R b i x pu (5)
where R b i is the rotation matrix from the base frame to the pulley frame. Following the ith loop of the CDPR geometry, the unwind cable length l ti can be expressed as:

l ci = m i m T i -r 2 p (6)
Then through eq. ( 3), the total length l ti can be obtained. The desired motor joint position is calculated given l ti and winch radius r w :

q d = l t r w (7)

Generation of nominal poses

Because the CDPR is under-actuated, equilibrium pose of its MP depends on both the cable lengths and gravity [START_REF] Wang | An Approach for Predicting the Calibration Accuracy in Planar Cable-Driven Parallel Robots[END_REF][START_REF] Abbasnejad | Direct Geometrico-static Problem of Underconstrained Cable-Driven Parallel Robots With n Cables[END_REF]. Each such pose corresponds to a minimum in the gravitational potential energy possible with a certain set of cable lengths.

In order to generate different MP poses, random cable lengths are first generated. According to the workspace size, the cable lengths vary between 1 m and 5 m. The MP pose is obtained by CDPR forward kinetostatics. More specifically, the MP y p axis coordinate p jy is minimized by Matlab fmincon function, while the MP static equilibrium is taken into account as an optimization constraint. Figure 5 shows the forces acting on the MP. The angle between cable i and axis x p of the platform frame is noted as ϕ i , and the z 2 axis coordinate of the corresponding anchor point is h a . Upon assuming that the MP center of mass is P , the sum of moments at P should remain zero because of static equilibrium, therefore:

m P = b p T 1 E T τ 1 cos(ϕ 1 -θ) sin(ϕ 1 -θ) + b p T 2 E T τ 2 cos(ϕ 2 -θ) sin(ϕ 2 -θ) = 0, with E = 0 -1 1 0 . ( 8
)
with the force static equilibrium:

f = τ 1 + τ 2 + G = τ 1 cos(ϕ 1 -θ) sin(ϕ 1 -θ) + τ 2 cos(ϕ 2 -θ) sin(ϕ 2 -θ) + G cos(-θ) sin(-θ) = 0 2 (9)

Sensors used and the accuracy tests

Two sensors are used to realise robot calibration: the WitMotion BWT61CL inclinometer and the SICK DT50-2 laser displacement sensor. The inclinometer measures the rotational angles along three axes in a Cartesian coordinate system, with a resolution of 0.001 • and a ±90 • measurement range. The laser displacement sensor measures the distance to a surface within the range from 200 mm to 5000 mm with a resolution of 1 mm. The calibration method highly relies on the sensor measurement quality. Better knowledge on the measurements will certainly contribute to an improved calibration accuracy. However, the measurement accuracy given by the sensor datasheet does not necessarily correspond to what can be observed in the lab. Therefore, in order to characterise our sensors, their measurements are recorded and compared to the predetermined angles and distances.

WitMotion BWT61CL inclinometer

The inclinometer is tested with a indexing head which is placed on a level surface plate, as shown in Fig. 7(a). The sensor is secured to a thick aluminum bar and is kept at the same position for each of the 10 repetitions at every single spindle angle. The indexing head used has a 40:1 ratio between the crank and the spindle, so that one turn of the crank results in a 9 • rotation of the spindle. Before testing, the inclinometer measurement is reset to 0 • . Then the inclinometer is tilted successively at 0 • , 9 where θ s is the inclinometer measurement, and θ r is the nominal spindle angle set by the indexing head. Figure 9(a) shows that the obtained errors e θ absolute values keep increasing and form a fairly smooth curve. According to the inclinometer measurement principle (Fig. 8), instead of measuring the tilt angle directly, the inclinometer measures the lateral component of the gravitational acceleration a, namely,

e θ = θ s -θ r , (10)
g sinθ = a ( 11
)
where g and a are the Euclidean norms of g and a, respectively. The total differentiation of Eq. 11 yields

g cosθdθ = da (12)
therefore, According to Eq. ( 13), e θ has a relation proportional to secθ. The curve fitting for the inclinometer measurements is conducted by choosing empirically the constant da g . As a result, da g = -0.24 fits e θ the best, especially for the range lower than 60 • . Then in order to obtain more accurate sensor measurements, only the generated poses with θ values below 60 • are chosen for later simulations and experiments. From data fitting, the inclinometer accuracy is determined to be ±0.1 • with a 0.05 • repeatability.

Figures 7(b) and 10 show the DT50 accuracy test setup photo and schematic, respectively. In Fig. 9(b), it can be seen that the e θ values are above 0 in most of measurement range, meaning that the measured distances are larger than the real ones. The curve does not show a clear trend of the error distribution throughout the sensor range. If we assume the error to be roughly uniformly distributed, then the average of all the measurement errors, e d,f it = 2.87 mm can be used over the whole range. If we take the mean e d value at each test position, the standard deviation of each measurement with respect to e d,f it is 1.83 mm. Upon assuming that the maximum measurement error is equal to 3 times the standard deviation, a 6 mm error tolerance value is obtained. According to the maximum differences among the 10 measurements at each fixed position, the repeatability of the DT50 is determined as 1.5 mm.

Floor elevation measurements

The floor surface of the CRAFT prototype is neither perfectly flat nor parallel with the x 2 axis of F 2 . Therefore the DT50 sensor measurements will be affected as they depend on the floor flatness. In order to improve the DT50 measurement accuracy inside F 2 , the floor elevation along the direction of x 2 is measured with a laser tracker. The laser tracker is placed next to the workspace plane and it emits laser to the laser tracker T-probe (Fig. 11), which measures 3D Cartesian coordinates with an accuracy of ±0.1 mm.

The projection of point P on the floor is G 0 . Starting from G 0 , a line along the direction of x 2 is drawn on the floor with distance markers every 100 mm. The floor elevation is measured at each of the distance markers with 5 times of repetition. All measurements present high repeatability of values within ±0.05 mm. The mean measurement values are calculated for each measured position and a spline representing the floor elevation is created from those mean values. The measured and calculated floor elevations are summarized in Fig. 12.

Robot position control

Figure 13 shows the position control scheme of the CDPR. During the robot movement, a sampling time of 0.001 s is used. Knowing the predefined MP poses x d , through the CDPR inverse kinematics, the desired motor joint position q d is calculated. In order to prevent thrusts of the MP movement between two poses, fifth-order polynomial interpolation is used to obtain the desired joint velocity qd and acceleration qd . I m is the motor and gearbox inertia, and Γ is the output torque.

Calibration methodology

The robot calibration is achieved while recording the signals from all the sensors throughout all the MP measurement poses. At the initial position, the cable lengths are denoted as l i0 . When the MP is moved to all the n different poses within the workspace, the cable length variations ∆l ij are recorded by the encoders. Thus, the actual cable lengths l ij are calculated by:

l ij = ∆l ij + l i0 (15)
The CDPR system variables are the outputs of the problem, including the Cartesian coordinates of cable exit points a ix , a iz , the initial cable lengths l i0 , the Cartesian coordinates of the MP center p jx , p jz and the MP rotation angle θ j . Among these variables, the laser displacement sensor measures p jz , and the inclinometer measures θ j . All the sensor measurements are the inputs of the identification problem. The rest of the output variables are considered to be the unknowns, included in a 3m + n dimensional vector x:

x = [a 1x a 1z , a 2x a 2z , l 10 l 20 , p 1x , . . . , p nx ] (

Non-linear least square method, based on solving non-linear equations to reduce the cable length residuals, is commonly used for robot calibration [START_REF] Zhang | Kinematic Calibration of Cable-Driven Parallel Robots Considering the Pulley Kinematics[END_REF][START_REF] Qian | Kinematic calibration of a cable-driven parallel robot for 3D printing[END_REF][START_REF] Sandretto | Calibration of a Fully-Constrained Parallel Cable-Driven Robot[END_REF]. With NLLS method, in this work, the identification problem is formulated while integrating the measurements from laser displacement sensor, inclinometer and motor encoders. Upon comparing the real cable length l ij from Eq.( 15) with the estimated cable length l ti from Eq.(3), a system of mn equations can be obtained [START_REF] Fortin-Côté | An admittance control scheme for haptic interfaces based on cable-driven parallel mechanisms[END_REF][START_REF] Picard | Pulleys and Force Sensors Influence on Payload Estimation of Cable-Driven Parallel Robots[END_REF]:

f ij (x) = (l pi + l ci ) 2 -(∆l ij + l i0 ) 2 = 0, i = 1, 2 , j = 1, . . . , n (17)
To solve the nonlinear system of equations defined by Eq.( 17), the number of inputs must be larger than or equal to the number of unknowns:

mn + 2n ≥ 2m + m + 3n (18)
As the CDPR has m = 2 cables, the least number of measurement poses is n = 6. A similar simulation process is used in [START_REF] Wang | Exit Point, Initial Length and Pose Self-calibration Method for Cable-Driven Parallel Robots[END_REF][START_REF] Fortin-Côté | An admittance control scheme for haptic interfaces based on cable-driven parallel mechanisms[END_REF]. Arbitrary errors are added on x to simulate the approximately known system variables. The identification problem is then formulated as the nonlinear least square problem:

min x m i=1 n j=1 f 2 ij (19)
The pre-defined real variable values x r are compared with the identified ones x * , to evaluate the identification accuracy:

δx k = x * k -x r,k , k = 1, 2, ..., 3m + n (20)
where δx k is the difference between the real and identified variable values.

4 Simulation of 3-Dof, 2-cable planar CDPR

Results of the simulated identification problem

Simulations with 6 to 50 measurement poses are performed. For each pose, the simulation is repeated 500 times with different sensor measurement errors. Based on the accuracy test results given in Sec.2, the applied accuracy and repeatability are generated as independent normally distributed random values, while assuming the sensor accuracy and repeatability ranges to be three times the standard deviation, and their mean to be zero. The dispersion or standard deviations of the obtained results δx k is defined as: Figure 14 shows the standard deviation of the x-coordinate of the 2nd cable exit point, the initial cable length of the 1st cable and the x-coordinate of the 4th MP pose. It is apparent that the larger the number of poses, the lower the standard deviation σ, the better the calibration accuracy. It should be noted that σ does not decrease after 50 measurement poses.

Simulation of sensor effects on identification errors

The two sensors provide different types of quantities (lengths and angles), which makes the effect of each sensor on the overall calibration quality difficult to discern. To resolve this problem, the simulations are carried out where the sensor measurement errors are eliminated one at a time. If one sensor is eliminated, its measurement error will be set to zero, and the simulation will be processed by considering the errors in the other sensor only. An example of this set of results for the 2nd exit point x-coordinate a 2x is shown in Fig. 19. When the SICK sensor measurement errors are eliminated, the identification errors are much reduced, the dispersion of the 500 results are significantly reduced, as shown in Fig. 19(c). The identification errors are much more concentrated around 0 mm, and within ± 6 mm. Fig. 20 shows that the errors for the initial cable lengths are even more affected, with identification errors inside ± 1 mm range, and the σ values decreased from 7.5 mm to 0.2 mm. It should be noted that the errors in the x-coordinates of the MP positions are less affected, reduced from ±10 mm to ±5 mm, as seen in Fig. 21.

On the other hand, when the WitMotion inclinometer is not considered, the identification errors do not change significantly, as shown in Fig. 19(b). It means that the identification quality of the geometric parameters of the planar CDPR, the initial cable lengths and the moving-platform poses are not very sensitive to measurement errors in WitMotion inclinometer. We conclude that one could need a more accurate distance sensor to improve the calibration accuracy. Conversely, one could probably afford to use a less accurate tilt sensor without affecting the calibration.

Overall simulation

The previous sections focused on the effects of the number of poses and the sensors over the overall identification accuracy. In this section, we replicate in simulation the experiment that was preformed in the laboratory, for the purpose of validating our calibration model. The sensor measurement error models are the same as those used in previous simulations and defined in sections 2.3.1 and 2.3.2. One single calibration using all the 50 measurement poses is performed. The identification errors δx are calculated, being defined as the differences between the identified and the true system variables. These system variables are listed in Eq. ( 16). The MP movement is also simulated, the control scheme code that will be able to move the real robot is used in the simulation. The overall simulation results showing the probability density of each individual system variable identification error are summarized in Fig. 22. All identification errors fall within the range of ±9 mm, which is consistent with the previous detailed simulation. The results give a standard deviation σ of 3.9 mm, and the 3σ interval value being ±11.6 mm. Among the identification errors, 89.8% are less than ±5 mm.

CDPR experiments and methodology

The calibration results may be affected by plenty of uncertainties. Among them, the sensor measurement errors and the floor altitude are taken into account and examined. Another uncertainty is the cable elongation during the experiment, which can be hard to predict [START_REF] Choi | Integrated and nonlinear dynamic model of a polymer cable for low-speed cable-driven parallel robots[END_REF]. Because of the structure of CDPRs, the changes in cable lengths directly affect the MP poses [START_REF] Piao | Open-loop position control of a polymer cabledriven parallel robot via a viscoelastic cable model for high payload workspaces[END_REF]. In order to eliminate the effect of cable elongation these a Leica AT-901 laser tracker system (LTS, seen in Fig. 23(a)) used to measure the actual MP poses accurately.

The laser tracker can work with either Spherically Mounted RetroReflector (SMRs, shown in Fig. 23(b)) or T-probe to measure 3D Cartesian coordinates. All the measurements have the accuracy of at least ±0.1 mm. Apart from point measurements, the Metrolog software allows line, plane and cylindrical surface measurements, all of which consist of sets of points. The software is also able to construct a Cartesian coordinate system with the former and express the subsequent measurements directly in it. Figure 24 represents the general experimental setup. The MP workspace plane is noted as π 2 . The laser tracker is placed on a line passing through the center of the rectangular workspace and perpendicular to π 2 . This minimizes the incidence angle of the laser with the workspace The position of the laser tracker next to the robot needs to be calibrated every time before an experiment session, to ensure the measurement accuracy. The base frame F 1 is defined by the laser tracker, with respect to datum marks arbitrarily fixed on the CRAFT prototype frame. Six SMR brackets are fixed on the CRAFT frame, their coordinates p r in F 1 are measured. The position of the laser tracker are known at the same time. Every time the laser tracker is moved again, the vectors p r are measured again. By comparing the old and new coordinates p r , the actual laser tracker position is known. Therefore the laser tracker is able to give measurements in F 1 accurately. After calibrating the laser tracker, the rest of the measurement procedures can be carried out. Section 5.1 introduces the construction of the MP frame F p by measuring necessary MP features. Afterwards, section 5.2 details the method to obtain MP pose measurements through the transformation between F p and F 1 .

Construction of the MP frame F p

The MP pose measurements are based on the three SMRs attached on the platform, whose accurate Cartesian coordinates in F p are needed. The methodology to construct the MP frame F p in Metrolog software is introduced in this section. The SMR coordinates can then be measured directly. Moreover, the Cartesian coordinates of the anchor points B 1 and B 2 are also measured.

In order to carry out the measurements, the platform is steadily placed on the stand, and the measurements are done using the laser tracker T-probe, as shown in Fig. 25(a). First, 10 points on the revolute joint shaft cylinder surface π s are measured, to obtain the parameters of axis y p . Then, 10 more points on the face of the MP strut that is perpendicular to y p . This plane is named π y . π y is obtained by offsetting π y by 20 mm along the -y p direction, which represents the vertical symmetry plane of the MP. The plane π y is perpendicular to y p and its intersection with y p is P , the origin of F p . 10 more points are measured on the MP top plane π z , and the projection of P on π z is noted as O pz . The vector from P to O pz is the z p axis of F p . Then the x p axis can be constructed as the line through P and perpendicular to y p and z p . This fully defines frame F p . The coordinates of the three SMRs on the MP, S 1 , S 2 and S 3 , can be directly measured.

Figure 26 shows the datums defined to measure the MP anchor points. First, the probe is placed on B i to obtain the probe center coordinate P bi . This represents the center of the probe, not yet the exact anchor point. Then 10 points are recorded on the flat bracket surface to obtain plane π bi . Then, the projection of P bi on π bi is assumed to be the actual anchor point B i . With the help of the laser tracker, the anchor points coordinates are measured accurately, which provides an accurate geometric model of the robot.

Transformation between F p and F 1

During the measurements, the MP is successively moved to the 50 different poses. For each pose, all the three SMRs are measured and their coordinates are expressed in F 1 . The Kabsch algorithm takes the three SMR coordinates in both F p and F 1 frame to obtain the homogeneous transformation matrix between the two frames, T

1 is also known. Therefore, the point coordinates in F p can be transferred into F 2 by:

p 1 h = T 1 p p p h (22)
p 2 h = T 2 1 p 1 h (23)
where p p h , p 1 h and p 2 h are the homogeneous coordinates of P in different frames. The MP orientations θ in F 2 are calculated from the two end points of the straight MP strut. These points are noted as C 1 and C 2 with coordinates [-b, 0, 0] T and [b, 0, 0] T in F p , respectively. Based on the previous transformation, their coordinates are also expressed in F 2 . Let the vector pointing from C 1 to C 2 be v 2 , then the MP rotational angle θ can be simply calculated as:

θ = atan2( v 2 z2 , v 2 x2 ) (24)
where v 6 Experimental results and analysis

Comparison between sensor and laser tracker measurements

Once the laser tracker measurements are obtained, they are compared with those from the inclinometer and laser displacement sensor. Figures 27 and28 show the differences on MP height δp jz and on MP inclination δθ j between the measured and desired values, which are calculated as follows:

δp jz = p zm -p zt , (25)
δθ j = θ zm -θ zt , (26)
where p zm is the sensor measured MP heights, from both the laser tracker and DT50 laser displacement sensor, and p zt is the desired values generated with the method described in Sec. 2.2. Similarly, θ zm and θ zt are the measured and desired MP inclination, respectively. The results reported in Fig. 27 show that both the laser tracker and the DT50 displacement sensor measure higher values of p z than the nominal ones, with error ranges from 0-35mm. On the other hand, the measured MP inclination shown in Fig. 28 are quite close to the nominal values, and are roughly evenly distributed on the positive and negative sides of the origin. The DT50 and inclinometer measurements are then compared in details with those from the laser tracker, as it gives very accurate measurements. the differences δp jz and δθ j are calculated as:

δp jz = p DT 50 -p LT S ( 27
)
δθ j = θ Inc -θ LT S (28)
where p DT 50 and p LT S are the MP height measurements from the DT50 and the laser tracker, respectively. The results are summarized in Figs. 29 and30. For almost all the poses, except for poses 24 and 31, the DT50 gives errors lower than 6 mm, which is the magnitude of the Although the DT50 measurements are within the expected measurement error range, the measured MP heights show relatively larger differences to the theoretical values. The possible reason for that is the cable elongations during experiments. As the cable winds on the winch, the winch exit point moves along the rotational axis direction of the winch, the real cable length therefore changes. And as for the exceptions of poses 24 and 31 on δp jz values, the positions of these two poses are high in the workspace (2 m and 1.6 m), because of the suspended configuration of the CDPR, the cable lengths are short and the cable tensions are relatively high at these two poses. In this situation, the displacement of the winch exit points, as well as the cable elongations under high cable tension may combine to affect the real MP poses.

Real cable lengths

The real cable lengths are calculated through the inverse kinematics from the MP pose measurements of the laser tracker, and compared with the theoretical cable lengths. The results are shown in Fig. 31. The errors on real cable lengths δl t do vary among a relatively large range, from -38 mm to 12 mm. The cable with a shorter real length tends to have a negative elongation. On the contrary, a longer cable tends to have a positive elongation.

Calibration results and prediction of calibration accuracy

A calibration process based on NLLS method, using experimental measurements is carried out. Among all the inputs of the problem, p jz and θ j take the measurements from the laser displacement sensors and the inclinometer, respectively. Because of the uncertainties on the actual cable lengths, the measurements from the laser tracker are used to estimate these quantities, instead of those from the cable encoders. The system variables are then identified in three separate , where the characterized measurement error models of both the laser displacement sensor and the inclinometer are applied. Compared with the simulation results, the errors on all the system variables fall within the estimated ranges. The experimental results validate the current model of CDPR calibration therefore validate the method of predicting the calibration accuracy. From the figure, among all the system variables, the simulation yields reasonable ranges for the exit point coordinates and for the pose x-coordinates. But for the initial cable lengths, the identification errors are much less than estimated, the simulation results seem overly pessimistic.

Conclusions and future work

In summary, this paper presents the simulation and experimental validation of the calibration method of a 3-DoF, 2-cable, planar CDPR. This method of simulation and its underlying model can be reused by others to predict the calibration accuracy of different CDPRs and measurement techniques. It is hoped that it can help refine the sensor selection process for prescribed accuracies of the robot geometry. The calibration method proposed here relies on the combination of a laser displacement sensor and an inclinometer embedded on the movingplatform, as well as the motor encoders. Detailed pulley geometric modelling is considered. The actual accuracies of the sensors are examined beforehand in order to predict the calibration Figure 32: Calibrated system variable identification errors results. The simulations show that with more measurement poses used, the identification errors of the exit point, pose Cartesian coordinates and the initial cable lengths are reduced, and are less dispersed. It turns out that the laser displacement sensor has a much larger influence than the inclinometer on the identification errors. Furthermore, amongst all system variables, the errors in initial cable lengths are the most sensitive to the laser displacement sensor measurement errors whereas the errors in the MP position coordinates are the least sensitive ones. Based on the sensor considered in this work, the system variable errors are all below ±9 mm, and most are below ±5 mm.

Experiments are carried out following the calibration method proposed. Sensor measurements from laser displacement sensor and inclinometer are recorded for calibration. Because of the complex cable elongations during experiments, the accurate cable length variations are obtained with the measurements from a laser tracker system, instead of cable encoders. The differences between the real and theoretical cable lengths are also calculated. The laser tracker also provides the ground truth of the actual MP poses to validate the calibrated system variable values. As a result, all of the calibrated robot system variable errors fall within the simulated ranges.

The cable elongations have significant effects on the MP poses, and consequently on the CDPR calibration accuracy as well [START_REF] Nanthacoumarane | Polymer Cable Characterization in Cable-Driven Parallel Robots. 25me Congrs Franais De Mcanique Nantes[END_REF][START_REF] Baklouti | Dynamic and oscillatory motions of cabledriven parallel robots based on a nonlinear cable tension model[END_REF]. In future work, the cable elongations will be compensated with elastic modelling, creep modelling, and the modelling of the winch cable exit point movement. In addition to elongations, considering cable mass will all improve the calibration quality. Future work will also deal with the determination of optimal measurement poses for CDPR calibration and other calibration methods with different sensor combinations.

  (a) MP in F 1 and F 2 (b) Detail of the real MP

Figure 1 :Figure 2 :

 12 Figure 1: CDPR under study equipped with a laser displacement sensor and an inclinometer

Figure 3 :Figure 4 :

 34 Figure 3: The ith loop of the CDPR under study

Figure 5 :

 5 Figure 5: Forces acting on the MP

Figure 6 :

 6 Figure 6: Randomly generated calibration poses

  (a) Inclinometer test setup (b) Laser displacement sensor test setup

Figure 7 :

 7 Figure 7: WitMotion BWT61CL inclinometer and SICK DT50-2 laser displacement sensor accuracy test setup

Figure 8 :

 8 Figure 8: Inclinometer working principle

  (a) Inclinometer accuracy (b) Laser displacement sensor accuracy

Figure 9 :

 9 Figure 9: WitMotion BWT61CL inclinometer and SICK DT50-2 laser displacement sensor accuracy curve fitting

Figure 10 :

 10 Figure 10: Laser displacement sensor accuracy test schematic

Figure 11 :

 11 Figure 11: Floor elevation measurement with laser tracker T-probe

Figure 12 :

 12 Figure 12: Floor elevation along the x 2 direction measured with the laser tracker

Figure 14 :

 14 Figure 14: The standard deviation results of 3 examples of the system variables

Figures 15 to 17

 17 Figure14shows the standard deviation of the x-coordinate of the 2nd cable exit point, the initial cable length of the 1st cable and the x-coordinate of the 4th MP pose. It is apparent that the larger the number of poses, the lower the standard deviation σ, the better the calibration accuracy. It should be noted that σ does not decrease after 50 measurement poses. Figures 15 to 17 summarize three examples of the identification errors of the system variables with different numbers of measurement poses used, and standard deviation σ. The results are

Figure 15 :Figure 16 :

 1516 Figure 15: Identification error results in the x-coordinates of the 2nd cable exit point, with different numbers of measurement poses and σ values

Figure 17 :

 17 Figure 17: Identification error results of the x-coordinates of the 4th MP position vector, with different numbers of measurement poses and σ values

Figure 18 :

 18 Figure 18: Simulated identification errors for some parameters

  (a) Both sensors, σ=4.1mm (b) SICK sensor, σ=4.1mm (c) Inclinometer, σ=1.6mm

Figure 19 :

 19 Figure 19: Identification error results in the x-coordinates of the 2nd exit point when 50 poses are used, comparison of sensor effects

  (a) Both sensors, σ=7.7mm (b) SICK sensor, σ=7.5mm (c) Inclinometer, σ=0.2mm

Figure 20 :

 20 Figure 20: Identification error results in the initial cable length of the 1st cable when 50 poses are used, comparison of sensor effects

Figure 21 :

 21 Figure 21: Identification error results in the x-coordinates of the 4th MP position vector when 50 poses are used, comparison of sensor effects

Figure 22 :

 22 Figure 22: Identification errors in all system variables

  (a) MP and laser tracker (b) MP equipped with SMRs

Figure 23 :

 23 Figure 23: MP pose measurement with laser tracker system setup

Figure 24 :

 24 Figure 24: Schematic of a MP pose measurement with laser tracker system

1 p

 1 . As the frames F 1 and F 2 are both known,

  (a) MP on stand (b) CAD model of the measured entities

Figure 25 :Figure 26 :

 2526 Figure 25: MP geometry entity measurements

2 z2 and v 2 x2 are the projection of v 2 on axes z 2

 22 and x 2 , respectively. The MP poses are thus obtained in the form of the MP center P and tilt angle θ.

Figure 27 :

 27 Figure 27: Differences between measured and desired MP height

Figure 28 :

 28 Figure 28: Differences between measured and desired MP inclination

Figure 29 :

 29 Figure 29: DT50 MP height measurements compared with laser tracker

Figure 31 :

 31 Figure 31: Errors on real cable lengths

Acknowledgements This work was supported by the ANR CRAFT project, grant ANR-18-CE10-0004, https://anr.fr/Project-ANR-18-CE10-0004. The first author of the paper is grateful for the support of China Scholarship Council (CSC Grant No.202008070051).

Funding: This study was funded by the ANR CRAFT project, grant ANR-18-CE10-0004, https://anr.fr/Project-ANR-18-CE10-0004 and by the support of China Scholarship Council (CSC Grant No.202008070051).

Conflict of Interest Statement

Compliance with Ethical Standards: Conflict of Interest: The authors declare that they have no conflict of interest.