Damien Six
email: damien.six@ls2n.fr

Sébastien Briot
email: sebastien.briot@ls2n.fr

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

Philippe Martinet
email: phlippe.martinet@ls2n.fr

The kinematics, dynamics and control of a flying parallel robot with three quadrotors

This paper deals with a novel flying mechanism. Inspired from parallel manipulators, this flying robot is composed of three quadrotors linked by a rigid articulated architecture composed of three legs and a platform. Associating quadrotor co-manipulation and rigid bodies, this new design offers novel possibilities for aerial robotics and manipulation. Previous work leads to the design of a flying parallel robot with two quadrotors. However, this robot did not allow the full control of the six degrees of freedom of its end-effector. With an additional quadrotor, this study seeks to obtain a full control of the platform position and orientation. To prove this property, the kinematic constraints are verified through screw theory. Then, the dynamic model is established and a decoupling property leads to the design of a specific controller for the platform and legs configurations. ADAMS/SIMULINK co-simulations validate the theoretical developments.

I. INTRODUCTION

The last decade has seen the development of a new generation of aerial robots combining a flying UAV (Unnamed Autonomous Vehicle) with a robotic manipulator. Several options have been explored, from one DoF (Degree of Freedom) manipulator [START_REF] Lindsey | Construction of Cubic Structures with Quadrotor Teams[END_REF] toward more complex serial [START_REF] Kim | Aerial manipulation using a quadrotor with a two DOF robotic arm[END_REF] [3] or parallel [START_REF] Danko | A parallel manipulator for mobile manipulating UAVs[END_REF] [START_REF] Kamel | Design and modeling of dexterous aerial manipulator[END_REF] robotic devices. The drawbacks of such solutions are a limited payload and autonomy due to the additional embedded motors and also a manipulability impacted by the quadrotor underactuation. To overcome the limits in payload, some studies investigated co-manipulation between several quadrotors in cable-load transportation applications [START_REF] Sreenath | Dynamics, Control and Planning for Cooperative Manipulation of Payloads Suspended by Cables from Multiple Quadrotor Robots[END_REF]. However, cables cannot apply any pushing forces to the environment, limiting the manipulation area under the drones. Such flying device can also be seen as a cable-driven parallel robot where the actuators have been replaced by quadrotors. Inspired by this analogy, a new type of flying robot can be developed from rigid parallel robots. This new type of flying robots, composed of quadrotors associated to a rigid articulated passive architecture, presents the following advantages [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF] • An enhanced payload, by spreading the efforts over several drones; • No additional embedded motors, the robot is controlled only by the quadrotors; • Reconfiguration of the passive architecture potentially allows to perform tasks under and over the drone; Fig. 1: A flying parallel robot with three quadrotors

• The effector can be placed away from the drone, limiting the perturbations due to rotor wash; • A large choice of leg topology, issued from the study of parallel robots, can lead to several properties of interest. A first design with two quadrotors, based on a planar parallel mechanism is described in [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF]. The end-effector of this robot was not fully controllable in the six DoFs of space: its evolution was limited to a vertical plane, thus restraining its applicability. This study proposes a new solution based on a spatial architecture with three drones, allowing control of all possible DoFs of the end-effector. The kinematic study in section II shows that this architecture offers full control of the platform position and orientation and it describes its singularity locus. Its dynamic model is developed in section III. Decoupling properties in the dynamic model are exploited to design a cascaded controller in section IV. The simulation results, illustrating the theoretical developments are given in section V.

II. KINEMATICS

Figure 1 illustrates the flying parallel robot under study. It consists of a moving platform connected to three quadrotors by means of three identical legs. F p (O p , x p , y p , z p ) is a frame attached to the moving platform. Its origin O p is the centroid of the platform. z p axis is normal to the platform top plane. Let AB be the distance between any set of points A and B. Three points R i (i = 1, 2, 3) are defined by their coordinates ---→ O p R i = [ OpRi cos γ i , OpRi sin γ i , 0] T in frame F p with γ i = {0, 2π/3, -2π/3} rad. A leg i is linked to the platform at point R i by a revolute joint of direction r i = [-sin γ i , cos γ i , 0] T expressed in frame F p . Each leg i is attached to a drone i by a spherical joint situated at O i (i = 1, 2, 3) (see Fig. 2). O i is also the center of mass of Fig. 2: Side view of the flying parallel robot drone i. The length of leg i is defined by RiOi . The actuation of the robot is provided by the three quadrotors.

Attaching the leg at the drone center of mass can be technically challenging, but this configuration is necessary to obtain the full decoupling property on the robot dynamics described in section III. On an experimental device, the centers of the spherical joint may not be exactly situated at the drones center of mass. Then, the centrifugal and Coriolis forces caused by the drones rotation motions will affect the dynamics of the passive architecture and generate additional coupling. With a large distance between the center of mass and the center of rotation, the perturbations induced by this coupling may be higher than the controller stability margin. In [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF], the expression of this coupling term is given for a flying robot with two quadrotors and an example illustrates its effect on the tracking precision. Thus, the robustness of the control law against the perturbations induced by this coupling requires to be studied. This point is not addressed in this paper.

A. Parameterization

Let F i (O i , x i , y i , z i ) be the frame attached to drone i (i = 1, 2, 3). Its origin O i is at the drone center of mass which is also the center of leg i spherical joint. x i axis bisects the angle defined by two arms of quadrotor i. z i axis is aligned with the quadrotor i propellers axes (see Fig. 1).

q is a vector of the robot generalized coordinates composed of (in the following order)

• x p , y p , z p , φ p , θ p , ψ p the position and orientation coordinates of the moving platform frame F p with respect to the world frame; • q 1 , q 2 , q 3 the relative angle between the platform plane and leg 1, 2 and 3, respectively (see Fig. 2);

• φ 1 , θ 1 , ψ 1 , φ 2 , θ 2 , ψ 2 , φ 3 , θ 3 , ψ 3 the orientation coordi-
nates of the drone frames F 1 , F 2 and F 3 with respect to the world frame. Note that in this paper, the rotation sequence of the Bryant angles [START_REF] Diebel | Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors[END_REF] is chosen such that

v = R x (φ i )R y (θ i )R z (ψ i )v i = 1, 2, 3 or p (1)
where v is a vector expressed in world frame, v is the same vector expressed in local frame F i and R u (α) is the elemental rotation about axis u though the angle α.

B. Mobility and singularity analyses

Joining quadrotors with a parallel kinematic architecture combines several challenges. One of them is to be ensured that the platform remains fully controllable. Another challenge is the drones' underactuation: they can only provide an upward thrust force and three independent torques for six DoFs. Classically, for the control of quadrotors, the underactuation is handled through an appropriate controller [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. In section IV, a cascaded control law is designed. Constraining each drone orientation coordinates by a desired thrust force, this control law allows each drone to apply a force f i at the associated leg tip (point O i ) in any desired direction. The torques provided by the quadrotors are not transmitted to the parallel architecture through the passive spherical joints. Assuming that such a control law is designed, the kinematic study of the flying robot can be reduced to the study of a passive kinematic architecture, composed of the end-effector and three legs attached to it by the three revolute joints. This architecture is actuated by a force f i applied at each leg tip in any desired direction.

The kinematic study will first focus on the end-effector mobility and its singularity analysis. Then, the results obtained are extended to the whole passive architecture, including the mobility of each leg.

1) End-effector mobility: Screw theory [10] [11] is suitable for the study of the instantaneous motion of parallel manipulators [START_REF] Kumar | Instantaneous kinematics of parallel-chain robotic mechanisms[END_REF]. A twist and a wrench are screws that represent the instantaneous motion of a rigid body and a system of forces and moments applied on a rigid body, respectively. Screw theory allows to compute the set of constraints, named actuation wrenches, provided by each leg to the platform when the actuators are locked. The rank of the overall wrench system, spanned by the actuation wrenches of all the legs, defines the mobility of the platform. If this system is of dimension six, then the platform if fully controlled. Otherwise, at least one degree of freedom of the platform is not controlled.

For each leg i (i = 1, 2, 3), four unit twists ζ ki define the instantaneous motions of the leg, one for the revolute joint motion (k = 1) and three for the spherical joint motions (k = 2, 3, 4). Any actuation wrench ξ of the leg is reciprocal to all the unit twists denoted ζ ki characterizing the displacements of the passive joints [START_REF] Kumar | Instantaneous kinematics of parallel-chain robotic mechanisms[END_REF], i.e. ξ T ζ ki = 0 for any k. This means that the virtual power developed by the wrench ξ along the direction of motion ζ ki is null. Let us consider one robot leg i. To simplify the expressions of the unit twists for this leg, let us define a new frame F li (O i , S1 i , S2 i , z li ) attached to this leg with: the origin O i is at the center of the spherical joint; axis S1 i pointing toward R i (along the leg); axis S2 i in the direction of the revolute joint axis (see Figs. 1 and2). The unit twist defining the motion of the passive revolute joint in frame F li at point O i is expressed as [START_REF] Hunt | Kinematic Geometry of Mechanisms[END_REF]

ζ 1i = 0 0 RiOi 0 1 0 T (2)
while the three unit twists defining the motion of the passive spherical joint in frame F li at point O i are given by

ζ 2i = 0 0 0 1 0 0 T (3) ζ 3i = 0 0 0 0 1 0 T (4) ζ 4i = 0 0 0 0 0 1 T (5)
In these twists, the first three components represent the direction of the translation velocity while the three last components represent the direction of the rotational velocity. For each leg i, only two unit actuation wrenches are reciprocal to all twists ζ ki . Their expressions in the frame F li at point O i are

ξ 1i = 1 0 0 0 0 0 T (6) ξ 2i = 0 1 0 0 0 0 T (7)
in which the first three components represent the direction of the force exerted on the platform and the three last components the direction of the moment. For each leg i (i = 1, 2, 3), ξ 1i and ξ 2i are two forces, exerted on the platform by the leg, applied at O i and directed along S1 i and S2 i , respectively. ξ 1i and ξ 2i can be seen as the two components of the force f i , applied by the drone i on the passive architecture, that are transmitted to the platform. The third component of f i along z li is not transmitted to the endeffector and do not cause any motion of the end-effector. This last component available for each leg i will be used to control the relative orientation between the leg and the platform (see section II-B.3). Now considering the three robot legs, the overall actuation wrench is given by Ξ = [ξ 11 , ξ 21 , ξ 12 , ξ 22 , ξ 13 , ξ 23 ], with all wrenches expressed in the same frame (e.g. the world frame). This overall wrench system is of dimension 6, then the endeffector is fully controllable. The next section introduces the singularity analysis, i.e. it studies the singular configurations for which the actuation wrench system locally loses its rank.

2) Singularity analysis of the end effector: There are some tools that define the conditions of degeneracy of a wrench system among which the Grassmann geometry [START_REF] Merlet | Parallel robots[END_REF] and the Grassmann-Cayley algebra [START_REF] Ben-Horin | Application of GrassmannCayley Algebra to Geometrical Interpretation of Parallel Robot Singularities[END_REF] [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF]. Regarding our particular case, the Grassmann-Cayley algebra was already used to prove the condition of singularities if the system of wrenches is composed of a triplet of two (independent) forces ξ 1i and ξ 2i (i = 1, 2, 3), ξ 1i and ξ 2i being applied at the same point O i [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF]. Let us define four planes

• P i (i = 1, 2, 3) is the plane passing through point O i
and containing the two forces being the resultants of the wrenches ξ 1i and ξ 2i (so it contains the vectors S1 i and S 2i );

• P 4 is the plane containing the three points O 1 , O 2 and O 3 .
Conditions of singularities appear if and only if all the four planes P i (i = 1, 2, 3, 4) intersect in (at least) a point, that can be at infinity [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF] (see Fig. 3). Let us illustrate this singularity condition for a case study ( OpRi = 0.1 m and RiOi = 0.8 m). Each plane P i is Fig. 3: Example of a singular configuration for the flying parallel robot. P 1 , P 2 , P 3 and P 4 are coplanar. defined by its equation

a i x + b i y + c i z + d i = 0 (8)
Let A be the matrix composed of the four plane coefficients

A =     a 1 b 1 c 1 d 1 a 2 b 2 c 2 d 2 a 3 b 3 c 3 d 3 a 4 b 4 c 4 d 4     (9)
The geometric condition of singularity, i.e. all planes intersect in at least a point, is equivalent to det A = 0.

The relative positions of the planes P 1 , P 2 , P 3 and P 4 depends only on the mobile robot internal DoFs. Then, the robot singular configurations are not affected by translation or rotation of the whole flying robot in space (i.e. x p , y p , z p , φ p , θ p and ψ p coordinates). The internal degrees of freedom, i.e. the legs orientation q 1 , q 2 and q 3 , are sufficient to study the singularity locus. To avoid auto-collision, each coordinates q i is restricted to range [-π/2; π/2]. The 2dimensional representations as function of q 2 and q 3 of the singularity locus for q 1 = 0, -π/6, -π/3 and -π/2 are shown in Fig. 4. The singularity locus is obtained through the numerical computation of | det A| for which we retain the values "close to zero". For this numerical analysis, a threshold is set at 5e-4. The results for positive values of q 1 can be obtained by a central symmetry in the origin of the figures. From those figures, we can see that a workspace with each q i restricted to range [-π/2; 0] is singularity free. However, a change of configuration with a platform under the quadrotors (all q i > 0) will require the crossing of a singularity. Still, the effective workspace of the robot could be extended using singularity crossing techniques already applied on standard parallel robots [START_REF] Six | A Controller Avoiding Dynamic Model Degeneracy of Parallel Robots During Singularity Crossing[END_REF].

3) Singularity analysis of the passive architecture (endeffector and legs): Only two components of each force provided by the drones are transmitted to the platform (see section II-B.1). Then, one component by leg remains available to control an additional degree of freedom: one pure force along z li . This remaining inputs on leg i will allow the control of the leg configuration q i (see Fig. 2). In section II-B.2, the singularity conditions have been studied for the end-effector alone. In this section, by expressing the Jacobian matrix for the whole passive structure (end-effector and legs), it will be proven that including those additional coordinates do not change the conditions of singularity. Let q p = [x p , y p , z p , φ p , θ p , ψ p , q 1 , q 2 , q 3 ] T be the generalized coordinates of the passive architecture only. Let v i (i = 1, 2, 3) be the velocity of the point O i (leg tip) expressed in world frame and J be the Jacobian matrix such that

a) q 1 = 0 rad (b) q 1 = -π/6 rad (c) q 1 = -π/3 rad (d) q 1 = -π/2 rad
v = J qp (10)
with v = v T 1 , v T 2 , v T 3
T . A singularity of this (9×9) matrix corresponds to a robot configuration where an infinitesimal motion of the passive architecture coordinates q p is possible even if the leg tips are fixed.

To study the singularities of this matrix, its is convenient to change the reference frame. For each i = (1, 2, 3), let v 1i (respectively v 2i and v 3i ) be the coordinate of v i along S1 i (respectively S2 i and z li ). The vector v l is defined by

v l = v 11 v 12 v 12 v 21 v 22 v 23 v 31 v 32 v 33 T (11)
Let J l be a Jacobian matrix such that

v l = J l qp (12)
J and J l are related through a full rank matrix defined by the rotations from the world frame to the local frames F li and the coordinate ordering defined in [START_REF] Hunt | Kinematic Geometry of Mechanisms[END_REF]. Skipping the mathematical calculations, J l takes the form

J l =     J r 0 6×3 R1O1 0 0 R2O2 0 0 0 R3O3     (13)
with J r linking the velocity of the leg tips along S1 i and S2 i only (i = 1, 2, 3) to the derivative of the end-effector coordinates, i.e.

v 11 v 12 v 13 v 21 v 22 v 23 T = J r qe (14)
q e = [x p , y p , z p , φ p , θ p , ψ p ] T is the vector of the endeffector coordinates. J r is thus the inverse Jacobian matrix of the passive architecture in the case where the drones are able to exert motions (but also forces) along S1 i and S2 i (i = 1, 2, 3) only. It is known from [START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF] that the rank of this Jacobian matrix is equal to the rank of the overall wrench system. This overall wrench, composed of units forces applied on the end-effector by the legs along S1 i and S2 i , is the one studied in section II-B.2. If the endeffector is not in a singular configuration, shown in section II-B.2, then J r is of full rank. By [START_REF] Merlet | Parallel robots[END_REF], the rank of the matrix J l is equal to the rank of J r plus 3, both matrices have the same conditions of rank deficiency. Then, the singularity locus of the whole passive architecture is equivalent to the singularity locus of the endeffector given in section II-B.2.

III. DYNAMICS

The major difficulty in controlling quadrotor is underactuation. A quadrotor can only produce three independent torques and an upward thrust force to control its six DoF in space [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. In a single quadrotor dynamic model, the attitude dynamics is independent from its position dynamics and is controlled only through torques, while the translational dynamics only depends on the thrust force and the drone orientation. This property is exploited to design cascaded controllers for quadrotors. Because the spherical joint centers are situated at the drones center of mass, the attitude dynamics of the quadrotors are not linked to the passive architecture dynamics, allowing to find the same kind of decoupling. In this section, the dynamic equations of the passive architecture and the quadrotor attitude are given. The following assumptions are made. On each quadrotor, four control inputs are available, one thrust input u 1i and three torques

τ i = [u 2i , u 3i , u 4i ]
T . The rotor dynamics is considered fast before the body dynamics and is neglected in the scope of this study. Disturbances from aerodynamic effects are also neglected in the model. In this section, the dynamic equations of the passive architecture alone will first be established. Then, the Newton-Euler equations, applied at the drones center of gravity, will give two set of dynamic equations: one linking the drones attitude dynamics to the input torques and one linking the passive architecture dynamics to the drone thrust forces.

A. Passive architecture dynamics

From the quadrotors, only one force f i per quadrotor i (i = 1, 2, 3) is transmitted to the passive architecture through the passive spherical joints at each leg tip. Let L(q p , qp ) the Lagrangian of the passive architecture. The Lagrangian equation of the architecture with three forces (f 1 , f 2 , f 3 ) applied on the leg tips are

d dt ∂L ∂ qp T - ∂L ∂q p T = J T   f 1 f 2 f 3   (15)
with J a Jacobian matrix defined in (10) that remains full rank out of singularities as described in section II-B.2.

The robot is now considered to remain away from singularities. Equation ( 15) can be written as [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]

M p qp + c p = J T   f 1 f 2 f 3   ( 16
)
where M p is the (9 × 9) generalized inertia matrix and c p the 9-dimensional vector of Coriolis, centrifugal and gravitational effects.

B. Quadrotor dynamics

Only reaction forces (-f 1 , -f 2 and -f 3 ) are transmitted to the drones from the passive architecture through the spherical joints. The Newton-Euler laws applied at the drone i center of mass give

m vi = mg -f i + t i (17) τ i = I i ωi + ω i × I i ω i (18)
where

• m the mass of one drone1 ;

• t i the thrust force of drone i;

• g the gravity vector

• τ i is the torque vector of quarotor i in local frame;

• ω i is the angular velocity of quarotor i in local frame ;

• I i is the inertia matrix of quadrotor i. The angular velocity of a quadrotor i is related to the derivative of its coordinates through the relation

ω i = T i ηi (19)
with

η i = [φ i , θ i , ψ i ] T and
T i =   cos θ i cos ψ i sin ψ i 0 -cos θ i sin ψ i cos ψ i 0 sin θ i 0 1   (20)
This matrix depends on the choice in the orientation representation which is defined in (1) for this study. Introducing [START_REF] Voos | Nonlinear control of a quadrotor micro-uav using feedbacklinearization[END_REF] in [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] gives

τ i = I i (T i ηi + Ṫi ηi ) + T i ηi × I i T i ηi (21)
Equation ( 21) is the attitude coordinate dynamic model for each drone as function of its input torques. This equation is independent from the thrust force and the passive architecture dynamics.

C. Whole robot dynamics

Equations ( 17) and ( 16) link the translational dynamics of the drones and the passive architecture dynamics through the reaction forces. Introducing [START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF] and the time derivative of [START_REF] Ball | A treatise on the theory of screws[END_REF] in [START_REF] Six | A Controller Avoiding Dynamic Model Degeneracy of Parallel Robots During Singularity Crossing[END_REF] gives, skipping the mathematical calculations,

Mq p + c = J T t (22) with • t = t T 1 , t T 2 , t T 3 T • M = M p + mJ T J • c = c p + mJ T J qp -mJ T g T , g
T , g T T Another way to compute directly this expression is to consider an equivalent passive architecture that has an additional mass corresponding to one quadrotor mass at each leg tip with only the thrust force of each drone applied to the leg tip. Equation (22) gives the dynamic equations of the passive structure coordinates q p as function of the drones thrust forces only. The thrust force of a drone i in world frame is related to its thrust input u 1i by

t i = R i 0 0 u 1i T (23) with R i = R x (φ i )R y (θ i )R z (ψ i ) the rotation matrix of the local frame F i .
Equations ( 21) and ( 22) are the two set of equations expressing the dynamics of the whole robot. They show a decoupling between the passive structure dynamics and the drones attitude dynamics. However, the decoupling is not complete as the thrust forces t in the passive structure dynamic equations (22) depends on the thrust inputs u 1i but also on the orientation coordinates of the drones (23). Section IV-A shows how this coupling can be handled by binding the desired drone coordinates to the desired thrust forces.

IV. CONTROL

The main challenge to control the full architecture is linked to the quadrotors underactuation. However, imposing a condition on the quadrotors orientation coordinates allows the definition of a virtual input that decouples the dynamics of the passive architecture from the drone attitude coordinates. Based on this property, a two steps controller is established. The global scheme of this controller is given in Fig. 5.

A. Decoupling for control of the passive architecture

The design of the passive architecture control law is based on feedback linearisation of the robot dynamics. Such controller is well known for the control of robot manipulators [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] but has also been used for the control quadrotors [START_REF] Voos | Nonlinear control of a quadrotor micro-uav using feedbacklinearization[END_REF].

From the robot dynamic model ( 22), let us define a virtual control input vector ν 1

ν 1 = M -1 (J T t -c) (24)
If such control vector exists, then introducing (24) in ( 22) gives the double integrator

qp = ν 1 (25)
Fig. 5: General controller scheme From (24), the desired thrust forces

t d = t d 1 T , t d 2 T , t d 3 T T
in world frame to obtain the virtual control input ν 1 are given by

t d = J -T (Mν 1 + c) (26)
From ( 23), the thrust force provided by each drone depends on its thrust input its orientation coordinates

t i =   sin θ i -sin φ i cos θ i cos φ i cos θ i   u 1i (27)
The desired thrust t d in (26) force may only be obtained for a given orientation of each quadrotor. From (27), this desired orientation for a quadrotor i is given by the following angles

φ d i = -atan2(t d yi , t d zi ) (28)
θ d i = asin t d xi ||t d i || 2 (29)
with

t d i = t d xi , t d yi , t d zi T . t zi > 0 is assumed in all flight
configurations (which is the case as long as the dynamic forces and the passive architecture reaction forces along z axis remain lower than the drone weight). The thrust input u 1i is defined by the projection of the desired thrust t d i in world frame onto the local z i axis of the drone

u 1i = sin θ i t d xi -cos θ i sin φ i t d yi + cos φ i cos θ i t d zi (30)
Suppose θ i = θ d i and φ i = φ d i , then introducing (28) and (29) in (30) gives, skipping the mathematical calculations

u 1i = ||t i d || 2 (31)
And introducing (28), ( 29) and (31) in (23)

t i = R i 0 0 u 1i T = t i d i = 1, 2, 3 (32)
Equation (32) proves that the desired thrust force t d as well as the virtual control input ν 1 (24) can be obtained as long as the desired attitudes angles defined by ( 28) and ( 29) are respected with thrust inputs given by (30) on each quadrotor.

B. Passive architecture control law

Suppose a trajectory q d p defined for the coordinates of the passive architecture. Let us design a control law for the virtual input ν 1 . As ν 1 is the control input of a double integrator (25), any control law that stabilizes a double integrator can be designed. A PD (Proportional-Derivative) controller is then designed

ν 1 = qd p -K p (q p -q d p ) -K d ( qp -qd p ) (33)
with K p , K d strict positive gains to ensure the closed-loop stability.

C. Quadrotors attitude control law

The control of the attitude coordinates of the quadrotors is based on Sliding Mode Control (SMC). SMC have already been applied for the control of a single quadrotor attitude in [START_REF] Xu | Sliding Mode Control of a Quadrotor Helicopter[END_REF]. The reference trajectories of the drone i (i = 1, 2, 3) attitude coordinates φ d i and θ d i are defined by ( 28) and (29). The reference trajectories for the last coordinate on each drone ψ d i remains free to be defined by the user. It can be exploited to avoid the collision between the drone arms and the legs. Those three reference trajectories for each drone forms the desired attitude vector

η d i = φ d i , θ d i , ψ d i T .
From (21), three auxiliary inputs (one for each drone) are defined

ν 2i = T -1 i I -1 i (τ i -T i ηi × I i T i ηi ) -T -1 i Ṫi ηi (34)
Introducing (34) in (21) gives

ν 2i = ηi (35)
Each auxiliary input ν 2i acts as a double integrator on the drone attitude coordinates η i . A fast convergence is required on the attitude closed-loop system as the desired angles are necessary to obtain the convergence on the passive architecture coordinates (see section IV-A). To ensure also robustness against model errors, a SMC is defined. Let σ i be the sliding variable for each quadrotor

σ i = ( ηi -ηd i ) + λ(η i -η d i ) (36)
with λ > 0. The sliding mode control law is defined by

ν 2i = ηd i -λ( ηi -ηd i ) -K s sign(σ i ) (37)
where K s > 0 and the sign(.) function of a n-dimensional vector σ is the n-dimensional vector s defined by

s j =      1 if σ j > 0 0 if σ j = 0 -1 if σ j < 0 (38)
where s j is the i th element of a vector s. Consider now the Lyapunov function candidate

V i = σ T i σ i (39)
Taking the first time derivative of (39) by considering (35), ( 36) and (37), we have

Vi = -2K s σ T i sign(σ i ) (40)
Vi is negative definite for K s > 0 therefore the stability and convergence of σ i toward zero is ensured. The stability is also ensured with any bounded disturbances by increasing the gain K s . The sliding mode control law ensures that the manifold σ i = 0 is reached in finite time. Once on the manifold, the convergence of the attitude coordinates η i towards the desired one η d i is guaranteed by (36) with λ > 0. One drawback of the sliding mode control is that it generates chattering on the control inputs. To reduce this chattering effect, the sign(.) function is replaced by the sat(.)

s j =      1 if σ j > σ j / if |σ j | ≤ -1 if σ j < - (41)
defines a thin boundary layer neighboring the sliding surface. The stability beyond this boundary layer is guaranteed, and the system remains in this boundary once entering it.

D. Overall control system

The overall control system consists in the inner attitude and the outer passive architecture control loops, forming the structure as shown in Fig. 5. The two control systems determine the twelve input variables u 1i , u 2i , u 3i and u 4i for i = (1, 2, 3). The feedback linearization for the passive architecture control law takes into account the dynamics nonlinear behavior of the robot. The sliding mode control for the attitude control law ensures a fast convergence of the attitude coordinates and robustness against disturbances.

V. SIMULATION

A simulator has been developed using co-simulation between the softwares ADAMS and MATLAB/SIMULINK. A discrete sample time of communication is set at 1kHz between the two softwares. The geometrical parameters for the simulation are OpRi = 0.1 m and RiOi = 0.8 m. The dynamics parameters are given in Table I. To test the robustness of the controller against model errors, an aerodynamic perturbation is implemented. On each drone a

a i = -Av i ( 42
)
where A is a diagonal friction matrix. The friction coefficients are the same for each direction (1 kg.s -1 ). Also, a white noise of amplitude 10 -2 (m or rad) is added to the robot state measurements. The coefficients of the control laws defined in section IV are K p = 4, K d = 4, λ = 15, K s = 10, = 1. Successive trajectories, defined by nine-degrees polynomials have been designed to test the several degrees-of-freedom of the robot (translation, rotation of the platform and legs orientation). Those trajectories and the tracking obtained are show in Fig. 6. Fig. 7 shows the drone 1 desired attitude coordinates and the tracking of those coordinates for the whole simulation. The controller robustly completes the mission under disturbances and sensor noise. Table II gives the mean of absolute value of the filtered error for each Figure 8 shows the thrust force and torque inputs on drone 2 during the whole simulation. The chattering effect is induced on torques by the high level of noise on the sensors. The parameter acts as a trade-off between stability and level of chattering effect. A video of the simulation is available in the attached media content. In this video content, an example of a task that can be performed by such device is also illustrated. For this task, the robot has to insert a tool into a hole and provide a rotation motion about the vertical axis while maintaining the effector in the hole. This task illustrates the ability of the system to perform a rotation of the platform along a direction without any coupling along the other coordinates.

VI. CONCLUSION

In this paper, a new type of flying robot was proposed. The robot is composed of three quadrotors linked by a rigid articulated passive architecture. This new flying robot offers the ability to control a platform position and orientation in space, with a high potential to perform aerial operations from this platform. The study of the dynamic model showed decoupling properties. Those properties have been exploited to design a cascaded controller adapted for this new flying robot. Simulations showed controller performance and robustness against noisy pose estimation and disturbances. With a white noise of amplitude 10 -2 (m or rad) on the robot state measurements, the controller showed a mean error lower than 7.5e -3 (m or rad) on the passive architecture coordinates and lower than 1.5e -3 rad on the orientation coordinates.

Future work includes exploring the application of reconfiguration techniques known for parallel robots on the flying robot to increase its workspace. Thus, more possibilities will be offered for the relative position of the drones and the platform. Physical interactions with the environment can also be explored, leading to applications in aerial manipulation. Finally, the extension of the theoretical aspects to a general architecture associated to any number of quadrotors will be studied. More than three quadrotors will offer redundancy to control the platform, enhancing the robustness of the whole system against quadrotor failures and facilitating reconfiguration tasks.

(

Fig. 4 :

 4 Fig. 4: Singularity locus (| det A| < 5e-4) for the passive architecture with OpRi = 0.1 m and RiOi = 0.8 m. The red square areas represent a singularity free workspace.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Desired trajectory and filtered tracked trajectory in simulation for successive motions of the passive architecture

along z 2 Fig. 8 :

 28 Fig. 8: Torque and forces input for drone 2 in simulation

TABLE II :

 II Mean of the absolute value of the filtered error on each coordinate obtained in simulation function. The sat(.) function of a n-dimensional vector σ is the n-dimensional vector s defined by

Masses of the three drones are considered equal to simplify notations, but the result can easily be extended to drones with various masses.


Damien Six
email: damien.six@ls2n.fr

Sébastien Briot
email: sebastien.briot@ls2n.fr

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

Philippe Martinet
email: phlippe.martinet@ls2n.fr

The kinematics, dynamics and control of a flying parallel robot with three quadrotors

This paper deals with a novel flying mechanism. Inspired from parallel manipulators, this flying robot is composed of three quadrotors linked by a rigid articulated architecture composed of three legs and a platform. Associating quadrotor co-manipulation and rigid bodies, this new design offers novel possibilities for aerial robotics and manipulation. Previous work leads to the design of a flying parallel robot with two quadrotors. However, this robot did not allow the full control of the six degrees of freedom of its end-effector. With an additional quadrotor, this study seeks to obtain a full control of the platform position and orientation. To prove this property, the kinematic constraints are verified through screw theory. Then, the dynamic model is established and a decoupling property leads to the design of a specific controller for the platform and legs configurations. ADAMS/SIMULINK co-simulations validate the theoretical developments.

I. INTRODUCTION

The last decade has seen the development of a new generation of aerial robots combining a flying UAV (Unnamed Autonomous Vehicle) with a robotic manipulator. Several options have been explored, from one DoF (Degree of Freedom) manipulator [START_REF] Lindsey | Construction of Cubic Structures with Quadrotor Teams[END_REF] toward more complex serial [START_REF] Kim | Aerial manipulation using a quadrotor with a two DOF robotic arm[END_REF] [3] or parallel [START_REF] Danko | A parallel manipulator for mobile manipulating UAVs[END_REF] [START_REF] Kamel | Design and modeling of dexterous aerial manipulator[END_REF] robotic devices. The drawbacks of such solutions are a limited payload and autonomy due to the additional embedded motors and also a manipulability impacted by the quadrotor underactuation. To overcome the limits in payload, some studies investigated co-manipulation between several quadrotors in cable-load transportation applications [START_REF] Sreenath | Dynamics, Control and Planning for Cooperative Manipulation of Payloads Suspended by Cables from Multiple Quadrotor Robots[END_REF]. However, cables cannot apply any pushing forces to the environment, limiting the manipulation area under the drones. Such flying device can also be seen as a cable-driven parallel robot where the actuators have been replaced by quadrotors. Inspired by this analogy, a new type of flying robot can be developed from rigid parallel robots. This new type of flying robots, composed of quadrotors associated to a rigid articulated passive architecture, presents the following advantages [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF] • An enhanced payload, by spreading the efforts over several drones; • No additional embedded motors, the robot is controlled only by the quadrotors; • Reconfiguration of the passive architecture potentially allows to perform tasks under and over the drone; Fig. 1: A flying parallel robot with three quadrotors

• The effector can be placed away from the drone, limiting the perturbations due to rotor wash; • A large choice of leg topology, issued from the study of parallel robots, can lead to several properties of interest. A first design with two quadrotors, based on a planar parallel mechanism is described in [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF]. The end-effector of this robot was not fully controllable in the six DoFs of space: its evolution was limited to a vertical plane, thus restraining its applicability. This study proposes a new solution based on a spatial architecture with three drones, allowing control of all possible DoFs of the end-effector. The kinematic study in section II shows that this architecture offers full control of the platform position and orientation and it describes its singularity locus. Its dynamic model is developed in section III. Decoupling properties in the dynamic model are exploited to design a cascaded controller in section IV. The simulation results, illustrating the theoretical developments are given in section V.

II. KINEMATICS

Figure 1 illustrates the flying parallel robot under study. It consists of a moving platform connected to three quadrotors by means of three identical legs. F p (O p , x p , y p , z p ) is a frame attached to the moving platform. Its origin O p is the centroid of the platform. z p axis is normal to the platform top plane. Let AB be the distance between any set of points A and B. Three points R i (i = 1, 2, 3) are defined by their coordinates ---→ O p R i = [ OpRi cos γ i , OpRi sin γ i , 0] T in frame F p with γ i = {0, 2π/3, -2π/3} rad. A leg i is linked to the platform at point R i by a revolute joint of direction r i = [-sin γ i , cos γ i , 0] T expressed in frame F p . Each leg i is attached to a drone i by a spherical joint situated at O i (i = 1, 2, 3) (see Fig. 2). O i is also the center of mass of Fig. 2: Side view of the flying parallel robot drone i. The length of leg i is defined by RiOi . The actuation of the robot is provided by the three quadrotors.

Attaching the leg at the drone center of mass can be technically challenging, but this configuration is necessary to obtain the full decoupling property on the robot dynamics described in section III. On an experimental device, the centers of the spherical joint may not be exactly situated at the drones center of mass. Then, the centrifugal and Coriolis forces caused by the drones rotation motions will affect the dynamics of the passive architecture and generate additional coupling. With a large distance between the center of mass and the center of rotation, the perturbations induced by this coupling may be higher than the controller stability margin. In [START_REF] Six | Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot[END_REF], the expression of this coupling term is given for a flying robot with two quadrotors and an example illustrates its effect on the tracking precision. Thus, the robustness of the control law against the perturbations induced by this coupling requires to be studied. This point is not addressed in this paper.

A. Parameterization

Let F i (O i , x i , y i , z i ) be the frame attached to drone i (i = 1, 2, 3). Its origin O i is at the drone center of mass which is also the center of leg i spherical joint. x i axis bisects the angle defined by two arms of quadrotor i. z i axis is aligned with the quadrotor i propellers axes (see Fig. 1).

q is a vector of the robot generalized coordinates composed of (in the following order)

• x p , y p , z p , φ p , θ p , ψ p the position and orientation coordinates of the moving platform frame F p with respect to the world frame; • q 1 , q 2 , q 3 the relative angle between the platform plane and leg 1, 2 and 3, respectively (see Fig. 2);

• φ 1 , θ 1 , ψ 1 , φ 2 , θ 2 , ψ 2 , φ 3 , θ 3 , ψ 3 the orientation coordi-
nates of the drone frames F 1 , F 2 and F 3 with respect to the world frame. Note that in this paper, the rotation sequence of the Bryant angles [START_REF] Diebel | Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors[END_REF] is chosen such that

v = R x (φ i )R y (θ i )R z (ψ i )v i = 1, 2, 3 or p (1)
where v is a vector expressed in world frame, v is the same vector expressed in local frame F i and R u (α) is the elemental rotation about axis u though the angle α.

B. Mobility and singularity analyses

Joining quadrotors with a parallel kinematic architecture combines several challenges. One of them is to be ensured that the platform remains fully controllable. Another challenge is the drones' underactuation: they can only provide an upward thrust force and three independent torques for six DoFs. Classically, for the control of quadrotors, the underactuation is handled through an appropriate controller [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. In section IV, a cascaded control law is designed. Constraining each drone orientation coordinates by a desired thrust force, this control law allows each drone to apply a force f i at the associated leg tip (point O i ) in any desired direction. The torques provided by the quadrotors are not transmitted to the parallel architecture through the passive spherical joints. Assuming that such a control law is designed, the kinematic study of the flying robot can be reduced to the study of a passive kinematic architecture, composed of the end-effector and three legs attached to it by the three revolute joints. This architecture is actuated by a force f i applied at each leg tip in any desired direction.

The kinematic study will first focus on the end-effector mobility and its singularity analysis. Then, the results obtained are extended to the whole passive architecture, including the mobility of each leg.

1) End-effector mobility: Screw theory [10] [11] is suitable for the study of the instantaneous motion of parallel manipulators [START_REF] Kumar | Instantaneous kinematics of parallel-chain robotic mechanisms[END_REF]. A twist and a wrench are screws that represent the instantaneous motion of a rigid body and a system of forces and moments applied on a rigid body, respectively. Screw theory allows to compute the set of constraints, named actuation wrenches, provided by each leg to the platform when the actuators are locked. The rank of the overall wrench system, spanned by the actuation wrenches of all the legs, defines the mobility of the platform. If this system is of dimension six, then the platform if fully controlled. Otherwise, at least one degree of freedom of the platform is not controlled.

For each leg i (i = 1, 2, 3), four unit twists ζ ki define the instantaneous motions of the leg, one for the revolute joint motion (k = 1) and three for the spherical joint motions (k = 2, 3, 4). Any actuation wrench ξ of the leg is reciprocal to all the unit twists denoted ζ ki characterizing the displacements of the passive joints [START_REF] Kumar | Instantaneous kinematics of parallel-chain robotic mechanisms[END_REF], i.e. ξ T ζ ki = 0 for any k. This means that the virtual power developed by the wrench ξ along the direction of motion ζ ki is null. Let us consider one robot leg i. To simplify the expressions of the unit twists for this leg, let us define a new frame F li (O i , S1 i , S2 i , z li ) attached to this leg with: the origin O i is at the center of the spherical joint; axis S1 i pointing toward R i (along the leg); axis S2 i in the direction of the revolute joint axis (see Figs. 1 and2). The unit twist defining the motion of the passive revolute joint in frame F li at point O i is expressed as [START_REF] Hunt | Kinematic Geometry of Mechanisms[END_REF]

ζ 1i = 0 0 RiOi 0 1 0 T (2)
while the three unit twists defining the motion of the passive spherical joint in frame F li at point O i are given by

ζ 2i = 0 0 0 1 0 0 T (3) ζ 3i = 0 0 0 0 1 0 T (4) ζ 4i = 0 0 0 0 0 1 T (5)
In these twists, the first three components represent the direction of the translation velocity while the three last components represent the direction of the rotational velocity. For each leg i, only two unit actuation wrenches are reciprocal to all twists ζ ki . Their expressions in the frame F li at point O i are

ξ 1i = 1 0 0 0 0 0 T (6) ξ 2i = 0 1 0 0 0 0 T (7)
in which the first three components represent the direction of the force exerted on the platform and the three last components the direction of the moment. For each leg i (i = 1, 2, 3), ξ 1i and ξ 2i are two forces, exerted on the platform by the leg, applied at O i and directed along S1 i and S2 i , respectively. ξ 1i and ξ 2i can be seen as the two components of the force f i , applied by the drone i on the passive architecture, that are transmitted to the platform. The third component of f i along z li is not transmitted to the endeffector and do not cause any motion of the end-effector. This last component available for each leg i will be used to control the relative orientation between the leg and the platform (see section II-B.3). Now considering the three robot legs, the overall actuation wrench is given by Ξ = [ξ 11 , ξ 21 , ξ 12 , ξ 22 , ξ 13 , ξ 23 ], with all wrenches expressed in the same frame (e.g. the world frame). This overall wrench system is of dimension 6, then the endeffector is fully controllable. The next section introduces the singularity analysis, i.e. it studies the singular configurations for which the actuation wrench system locally loses its rank.

2) Singularity analysis of the end effector: There are some tools that define the conditions of degeneracy of a wrench system among which the Grassmann geometry [START_REF] Merlet | Parallel robots[END_REF] and the Grassmann-Cayley algebra [START_REF] Ben-Horin | Application of GrassmannCayley Algebra to Geometrical Interpretation of Parallel Robot Singularities[END_REF] [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF]. Regarding our particular case, the Grassmann-Cayley algebra was already used to prove the condition of singularities if the system of wrenches is composed of a triplet of two (independent) forces ξ 1i and ξ 2i (i = 1, 2, 3), ξ 1i and ξ 2i being applied at the same point O i [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF]. Let us define four planes

• P i (i = 1, 2, 3) is the plane passing through point O i
and containing the two forces being the resultants of the wrenches ξ 1i and ξ 2i (so it contains the vectors S1 i and S 2i );

• P 4 is the plane containing the three points O 1 , O 2 and O 3 .
Conditions of singularities appear if and only if all the four planes P i (i = 1, 2, 3, 4) intersect in (at least) a point, that can be at infinity [START_REF]Singularity condition of six-degree-of-freedom three-legged parallel robots based on Grassmann -Cayley Algebra[END_REF] (see Fig. 3). Let us illustrate this singularity condition for a case study ( OpRi = 0.1 m and RiOi = 0.8 m). Each plane P i is Fig. 3: Example of a singular configuration for the flying parallel robot. P 1 , P 2 , P 3 and P 4 are coplanar. defined by its equation

a i x + b i y + c i z + d i = 0 (8)
Let A be the matrix composed of the four plane coefficients

A =     a 1 b 1 c 1 d 1 a 2 b 2 c 2 d 2 a 3 b 3 c 3 d 3 a 4 b 4 c 4 d 4     (9)
The geometric condition of singularity, i.e. all planes intersect in at least a point, is equivalent to det A = 0.

The relative positions of the planes P 1 , P 2 , P 3 and P 4 depends only on the mobile robot internal DoFs. Then, the robot singular configurations are not affected by translation or rotation of the whole flying robot in space (i.e. x p , y p , z p , φ p , θ p and ψ p coordinates). The internal degrees of freedom, i.e. the legs orientation q 1 , q 2 and q 3 , are sufficient to study the singularity locus. To avoid auto-collision, each coordinates q i is restricted to range [-π/2; π/2]. The 2dimensional representations as function of q 2 and q 3 of the singularity locus for q 1 = 0, -π/6, -π/3 and -π/2 are shown in Fig. 4. The singularity locus is obtained through the numerical computation of | det A| for which we retain the values "close to zero". For this numerical analysis, a threshold is set at 5e-4. The results for positive values of q 1 can be obtained by a central symmetry in the origin of the figures. From those figures, we can see that a workspace with each q i restricted to range [-π/2; 0] is singularity free. However, a change of configuration with a platform under the quadrotors (all q i > 0) will require the crossing of a singularity. Still, the effective workspace of the robot could be extended using singularity crossing techniques already applied on standard parallel robots [START_REF] Six | A Controller Avoiding Dynamic Model Degeneracy of Parallel Robots During Singularity Crossing[END_REF].

3) Singularity analysis of the passive architecture (endeffector and legs): Only two components of each force provided by the drones are transmitted to the platform (see section II-B.1). Then, one component by leg remains available to control an additional degree of freedom: one pure force along z li . This remaining inputs on leg i will allow the control of the leg configuration q i (see Fig. 2). In section II-B.2, the singularity conditions have been studied for the end-effector alone. In this section, by expressing the Jacobian matrix for the whole passive structure (end-effector and legs), it will be proven that including those additional coordinates do not change the conditions of singularity. Let q p = [x p , y p , z p , φ p , θ p , ψ p , q 1 , q 2 , q 3 ] T be the generalized coordinates of the passive architecture only. Let v i (i = 1, 2, 3) be the velocity of the point O i (leg tip) expressed in world frame and J be the Jacobian matrix such that

a) q 1 = 0 rad (b) q 1 = -π/6 rad (c) q 1 = -π/3 rad (d) q 1 = -π/2 rad
v = J qp (10)
with v = v T 1 , v T 2 , v T 3
T . A singularity of this (9×9) matrix corresponds to a robot configuration where an infinitesimal motion of the passive architecture coordinates q p is possible even if the leg tips are fixed.

To study the singularities of this matrix, its is convenient to change the reference frame. For each i = (1, 2, 3), let v 1i (respectively v 2i and v 3i ) be the coordinate of v i along S1 i (respectively S2 i and z li ). The vector v l is defined by

v l = v 11 v 12 v 12 v 21 v 22 v 23 v 31 v 32 v 33 T (11)
Let J l be a Jacobian matrix such that

v l = J l qp (12)
J and J l are related through a full rank matrix defined by the rotations from the world frame to the local frames F li and the coordinate ordering defined in [START_REF] Hunt | Kinematic Geometry of Mechanisms[END_REF]. Skipping the mathematical calculations, J l takes the form

J l =     J r 0 6×3 R1O1 0 0 R2O2 0 0 0 R3O3     (13)
with J r linking the velocity of the leg tips along S1 i and S2 i only (i = 1, 2, 3) to the derivative of the end-effector coordinates, i.e.

v 11 v 12 v 13 v 21 v 22 v 23 T = J r qe (14)
q e = [x p , y p , z p , φ p , θ p , ψ p ] T is the vector of the endeffector coordinates. J r is thus the inverse Jacobian matrix of the passive architecture in the case where the drones are able to exert motions (but also forces) along S1 i and S2 i (i = 1, 2, 3) only. It is known from [START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF] that the rank of this Jacobian matrix is equal to the rank of the overall wrench system. This overall wrench, composed of units forces applied on the end-effector by the legs along S1 i and S2 i , is the one studied in section II-B.2. If the endeffector is not in a singular configuration, shown in section II-B.2, then J r is of full rank. By [START_REF] Merlet | Parallel robots[END_REF], the rank of the matrix J l is equal to the rank of J r plus 3, both matrices have the same conditions of rank deficiency. Then, the singularity locus of the whole passive architecture is equivalent to the singularity locus of the endeffector given in section II-B.2.

III. DYNAMICS

The major difficulty in controlling quadrotor is underactuation. A quadrotor can only produce three independent torques and an upward thrust force to control its six DoF in space [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. In a single quadrotor dynamic model, the attitude dynamics is independent from its position dynamics and is controlled only through torques, while the translational dynamics only depends on the thrust force and the drone orientation. This property is exploited to design cascaded controllers for quadrotors. Because the spherical joint centers are situated at the drones center of mass, the attitude dynamics of the quadrotors are not linked to the passive architecture dynamics, allowing to find the same kind of decoupling. In this section, the dynamic equations of the passive architecture and the quadrotor attitude are given. The following assumptions are made. On each quadrotor, four control inputs are available, one thrust input u 1i and three torques

τ i = [u 2i , u 3i , u 4i ]
T . The rotor dynamics is considered fast before the body dynamics and is neglected in the scope of this study. Disturbances from aerodynamic effects are also neglected in the model. In this section, the dynamic equations of the passive architecture alone will first be established. Then, the Newton-Euler equations, applied at the drones center of gravity, will give two set of dynamic equations: one linking the drones attitude dynamics to the input torques and one linking the passive architecture dynamics to the drone thrust forces.

A. Passive architecture dynamics

From the quadrotors, only one force f i per quadrotor i (i = 1, 2, 3) is transmitted to the passive architecture through the passive spherical joints at each leg tip. Let L(q p , qp ) the Lagrangian of the passive architecture. The Lagrangian equation of the architecture with three forces (f 1 , f 2 , f 3 ) applied on the leg tips are

d dt ∂L ∂ qp T - ∂L ∂q p T = J T   f 1 f 2 f 3   (15)
with J a Jacobian matrix defined in (10) that remains full rank out of singularities as described in section II-B.2.

The robot is now considered to remain away from singularities. Equation ( 15) can be written as [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]

M p qp + c p = J T   f 1 f 2 f 3   ( 16
)
where M p is the (9 × 9) generalized inertia matrix and c p the 9-dimensional vector of Coriolis, centrifugal and gravitational effects.

B. Quadrotor dynamics

Only reaction forces (-f 1 , -f 2 and -f 3 ) are transmitted to the drones from the passive architecture through the spherical joints. The Newton-Euler laws applied at the drone i center of mass give

m vi = mg -f i + t i (17) τ i = I i ωi + ω i × I i ω i (18)
where

• m the mass of one drone1 ;

• t i the thrust force of drone i;

• g the gravity vector

• τ i is the torque vector of quarotor i in local frame;

• ω i is the angular velocity of quarotor i in local frame ;

• I i is the inertia matrix of quadrotor i. The angular velocity of a quadrotor i is related to the derivative of its coordinates through the relation

ω i = T i ηi (19)
with

η i = [φ i , θ i , ψ i ] T and
T i =   cos θ i cos ψ i sin ψ i 0 -cos θ i sin ψ i cos ψ i 0 sin θ i 0 1   (20)
This matrix depends on the choice in the orientation representation which is defined in (1) for this study. Introducing [START_REF] Voos | Nonlinear control of a quadrotor micro-uav using feedbacklinearization[END_REF] in [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] gives

τ i = I i (T i ηi + Ṫi ηi ) + T i ηi × I i T i ηi (21)
Equation ( 21) is the attitude coordinate dynamic model for each drone as function of its input torques. This equation is independent from the thrust force and the passive architecture dynamics.

C. Whole robot dynamics

Equations ( 17) and ( 16) link the translational dynamics of the drones and the passive architecture dynamics through the reaction forces. Introducing [START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF] and the time derivative of [START_REF] Ball | A treatise on the theory of screws[END_REF] in [START_REF] Six | A Controller Avoiding Dynamic Model Degeneracy of Parallel Robots During Singularity Crossing[END_REF] gives, skipping the mathematical calculations,

Mq p + c = J T t (22) with • t = t T 1 , t T 2 , t T 3 T • M = M p + mJ T J • c = c p + mJ T J qp -mJ T g T , g
T , g T T Another way to compute directly this expression is to consider an equivalent passive architecture that has an additional mass corresponding to one quadrotor mass at each leg tip with only the thrust force of each drone applied to the leg tip. Equation (22) gives the dynamic equations of the passive structure coordinates q p as function of the drones thrust forces only. The thrust force of a drone i in world frame is related to its thrust input u 1i by

t i = R i 0 0 u 1i T (23) with R i = R x (φ i )R y (θ i )R z (ψ i ) the rotation matrix of the local frame F i .
Equations ( 21) and ( 22) are the two set of equations expressing the dynamics of the whole robot. They show a decoupling between the passive structure dynamics and the drones attitude dynamics. However, the decoupling is not complete as the thrust forces t in the passive structure dynamic equations (22) depends on the thrust inputs u 1i but also on the orientation coordinates of the drones (23). Section IV-A shows how this coupling can be handled by binding the desired drone coordinates to the desired thrust forces.

IV. CONTROL

The main challenge to control the full architecture is linked to the quadrotors underactuation. However, imposing a condition on the quadrotors orientation coordinates allows the definition of a virtual input that decouples the dynamics of the passive architecture from the drone attitude coordinates. Based on this property, a two steps controller is established. The global scheme of this controller is given in Fig. 5.

A. Decoupling for control of the passive architecture

The design of the passive architecture control law is based on feedback linearisation of the robot dynamics. Such controller is well known for the control of robot manipulators [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] but has also been used for the control quadrotors [START_REF] Voos | Nonlinear control of a quadrotor micro-uav using feedbacklinearization[END_REF].

From the robot dynamic model ( 22), let us define a virtual control input vector ν 1

ν 1 = M -1 (J T t -c) (24)
If such control vector exists, then introducing (24) in ( 22) gives the double integrator

qp = ν 1 (25)
Fig. 5: General controller scheme From (24), the desired thrust forces

t d = t d 1 T , t d 2 T , t d 3 T T
in world frame to obtain the virtual control input ν 1 are given by

t d = J -T (Mν 1 + c) (26)
From ( 23), the thrust force provided by each drone depends on its thrust input and orientation coordinates

t i =   sin θ i -sin φ i cos θ i cos φ i cos θ i   u 1i (27)
The desired thrust t d in (26) force may only be obtained for a given orientation of each quadrotor. From (27), this desired orientation for a quadrotor i is given by the following angles

φ d i = -atan2(t d yi , t d zi ) (28)
θ d i = asin t d xi ||t d i || 2 (29)
with

t d i = t d xi , t d yi , t d zi T . t zi > 0 is assumed in all flight
configurations (which is the case as long as the dynamic forces and the passive architecture reaction forces along z axis remain lower than the drone weight). The thrust input u 1i is defined by the projection of the desired thrust t d i in world frame onto the local z i axis of the drone

u 1i = sin θ i t d xi -cos θ i sin φ i t d yi + cos φ i cos θ i t d zi (30)
Suppose θ i = θ d i and φ i = φ d i , then introducing (28) and (29) in (30) gives, skipping the mathematical calculations

u 1i = ||t i d || 2 (31)
And introducing (28), ( 29) and (31) in (23)

t i = R i 0 0 u 1i T = t i d i = 1, 2, 3 (32)
Equation (32) proves that the desired thrust force t d as well as the virtual control input ν 1 (24) can be obtained as long as the desired attitudes angles defined by ( 28) and ( 29) are respected with thrust inputs given by (30) on each quadrotor.

B. Passive architecture control law

Suppose a trajectory q d p defined for the coordinates of the passive architecture. Let us design a control law for the virtual input ν 1 . As ν 1 is the control input of a double integrator (25), any control law that stabilizes a double integrator can be designed. A PD (Proportional-Derivative) controller is then designed

ν 1 = qd p -K p (q p -q d p ) -K d ( qp -qd p ) (33)
with K p , K d strict positive gains to ensure the closed-loop stability.

C. Quadrotors attitude control law

The control of the attitude coordinates of the quadrotors is based on Sliding Mode Control (SMC). SMC have already been applied for the control of a single quadrotor attitude in [START_REF] Xu | Sliding Mode Control of a Quadrotor Helicopter[END_REF]. The reference trajectories of the drone i (i = 1, 2, 3) attitude coordinates φ d i and θ d i are defined by ( 28) and (29). The reference trajectories for the last coordinate on each drone ψ d i remains free to be defined by the user. It can be exploited to avoid the collision between the drone arms and the legs. Those three reference trajectories for each drone forms the desired attitude vector

η d i = φ d i , θ d i , ψ d i T .
From (21), three auxiliary inputs (one for each drone) are defined

ν 2i = T -1 i I -1 i (τ i -T i ηi × I i T i ηi ) -T -1 i Ṫi ηi (34)
Introducing (34) in (21) gives

ν 2i = ηi (35)
Each auxiliary input ν 2i acts as a double integrator on the drone attitude coordinates η i . A fast convergence is required on the attitude closed-loop system as the desired angles are necessary to obtain the convergence on the passive architecture coordinates (see section IV-A). To ensure also robustness against model errors, a SMC is defined. Let σ i be the sliding variable for each quadrotor

σ i = ( ηi -ηd i ) + λ(η i -η d i ) (36)
with λ > 0. The sliding mode control law is defined by

ν 2i = ηd i -λ( ηi -ηd i ) -K s sign(σ i ) (37)
where K s > 0 and the sign(.) function of a n-dimensional vector σ is the n-dimensional vector s defined by

s j =      1 if σ j > 0 0 if σ j = 0 -1 if σ j < 0 (38)
where s j is the i th element of a vector s. Consider now the Lyapunov function candidate

V i = σ T i σ i (39)
Taking the first time derivative of (39) by considering (35), ( 36) and (37), we have

Vi = -2K s σ T i sign(σ i ) (40)
Vi is negative definite for K s > 0 therefore the stability and convergence of σ i toward zero is ensured. The stability is also ensured with any bounded disturbances by increasing the gain K s . The sliding mode control law ensures that the manifold σ i = 0 is reached in finite time. Once on the manifold, the convergence of the attitude coordinates η i towards the desired one η d i is guaranteed by (36) with λ > 0. One drawback of the sliding mode control is that it generates chattering on the control inputs. To reduce this chattering effect, the sign(.) function is replaced by the sat(.)

s j =      1 if σ j > σ j / if |σ j | ≤ -1 if σ j < - (41)
defines a thin boundary layer neighboring the sliding surface. The stability beyond this boundary layer is guaranteed, and the system remains in this boundary once entering it.

D. Overall control system

The overall control system consists in the inner attitude and the outer passive architecture control loops, forming the structure as shown in Fig. 5. The two control systems determine the twelve input variables u 1i , u 2i , u 3i and u 4i for i = (1, 2, 3). The feedback linearization for the passive architecture control law takes into account the dynamics nonlinear behavior of the robot. The sliding mode control for the attitude control law ensures a fast convergence of the attitude coordinates and robustness against disturbances.

V. SIMULATION

A simulator has been developed using co-simulation between the softwares ADAMS and MATLAB/SIMULINK. A discrete sample time of communication is set at 1kHz between the two softwares. The geometrical parameters for the simulation are OpRi = 0.1 m and RiOi = 0.8 m. The dynamics parameters are given in Table I. To test the robustness of the controller against model errors, an aerodynamic perturbation is implemented. On each drone a

a i = -Av i ( 42
)
where A is a diagonal friction matrix. The friction coefficients are the same for each direction (1 kg.s -1 ). Also, a white noise of amplitude 10 -2 (m or rad) is added to the robot state measurements. The coefficients of the control laws defined in section IV are K p = 4, K d = 4, λ = 15, K s = 10, = 1. Successive trajectories, defined by nine-degrees polynomials have been designed to test the several degrees-of-freedom of the robot (translation, rotation of the platform and legs orientation). Those trajectories and the tracking obtained are show in Fig. 6. Fig. 7 shows the drone 1 desired attitude coordinates and the tracking of those coordinates for the whole simulation. The controller robustly completes the mission under disturbances and sensor noise. Table II gives the mean of absolute value of the filtered error for each Figure 8 shows the thrust force and torque inputs on drone 2 during the whole simulation. The chattering effect is induced on torques by the high level of noise on the sensors. The parameter acts as a trade-off between stability and level of chattering effect. A video of the simulation is available in the attached media content. In this video content, an example of a task that can be performed by such device is also illustrated. For this task, the robot has to insert a tool into a hole and provide a rotation motion about the vertical axis while maintaining the effector in the hole. This task illustrates the ability of the system to perform a rotation of the platform along a direction without any coupling along the other coordinates.

VI. CONCLUSION

In this paper, a new type of flying robot was proposed. The robot is composed of three quadrotors linked by a rigid articulated passive architecture. This new flying robot offers the ability to control a platform position and orientation in space, with a high potential to perform aerial operations from this platform. The study of the dynamic model showed decoupling properties. Those properties have been exploited to design a cascaded controller adapted for this new flying robot. Simulations showed controller performance and robustness against noisy pose estimation and disturbances. With a white noise of amplitude 10 -2 (m or rad) on the robot state measurements, the controller showed a mean error lower than 7.5e -3 (m or rad) on the passive architecture coordinates and lower than 1.5e -3 rad on the orientation coordinates.

Future work includes exploring the application of reconfiguration techniques known for parallel robots on the flying robot to increase its workspace. Thus, more possibilities will be offered for the relative position of the drones and the platform. Physical interactions with the environment can also be explored, leading to applications in aerial manipulation. Finally, the extension of the theoretical aspects to a general architecture associated to any number of quadrotors will be studied. More than three quadrotors will offer redundancy to control the platform, enhancing the robustness of the whole system against quadrotor failures and facilitating reconfiguration tasks.

(

Fig. 4 :

 4 Fig. 4: Singularity locus (| det A| < 5e-4) for the passive architecture with OpRi = 0.1 m and RiOi = 0.8 m. The red square areas represent a singularity free workspace.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Desired trajectory and filtered tracked trajectory in simulation for successive motions of the passive architecture

along z 2 Fig. 8 :

 28 Fig. 8: Torque and forces input for drone 2 in simulation

TABLE II :

 II Mean of the absolute value of the filtered error on each coordinate obtained in simulation function. The sat(.) function of a n-dimensional vector σ is the n-dimensional vector s defined by

Masses of the three drones are considered equal to simplify notations, but the result can easily be extended to drones with various masses.


Lila Kaci

Sébastien Briot

Clément Boudaud

Philippe Martinet

Design of a Wooden Five-bar Mechanism

Eco-design of robots is a field of research which has been rarely explored in the past. In order to considerably decrease the environmental impact of robot during the design phase, metal or carbon composite parts can be replaced by bio-sourced materials, such as wood. Indeed, wood has interesting mechanical properties, but its performance / dimensions will vary with the atmospheric conditions / external solicitations and with the conditions in which trees have grown. This paper deals with the design of a stiff and accurate wooden five-bar mechanism. First, a control-based design problem is formulated. This problem aims at finding the optimal parameters of the robot, taking into account the nature of the desired control (sensor-based control). Then, reliable topology optimization methodology is proposed, taking into account the variability of the wood performance and that will allow the definition of robot architectures (shape of the links) for which the impact of this variability in terms of deformation is minimal. Finally, the optimal design variables are given and are used for the realisation of industrial prototype of a wooden five-bar mechanism.

Introduction

The Climate Change Mitigation (CCM) has become a priority in the world, as shown during the COP21. It becomes urgent to decrease the Environmental Impact (EI) of the human activities. Robotics, which takes an increasing place in our everyday life, must also be part of this effort. A recent study (Fiz) shows that, contrary to what could be imagined, a large percentage (around 50 %) of the EI during the entire robot life is due to the use of metallic materials (steel or aluminum, even if they are recycled) for the design of their links. The percentage is even worth if carbon composites are used. The rest of the EI is globally due to the robot energy consumption during use. A method to reduce EI of robot design, is to replace the polluting metals and carbon composites with bio-sourced materials (BSM) that have a little (even no) EI [START_REF] Thakur | Green Composites from Natural Resources[END_REF]). Moreover, BSM have good stiffness performance, and their quantity is only limited by their growing time [START_REF] Kretschmann | Ch. 5: Mechanical properties of Wood[END_REF]). In this work, developped in the scope of the French project RobEcolo (Rob), we are interested in using wood instead of metals to design the robot links.

Using wood in machines is not a new idea. For instance, it is used in chassis of cars (Mor). It found also applications in the design of buildings due to its stiffness and low cost [START_REF] Kretschmann | Ch. 5: Mechanical properties of Wood[END_REF]). It is also used in Robotics for the design of mock-ups and proof-of-concept prototypes [START_REF] Quigley | A low-cost compliant 7-DOF robotic manipulator[END_REF]). However, a detailed study of these works shows that wood is never used in critical parts ensuring accuracy. Indeed, the wood performance / dimensions will vary with the atmospheric conditions / external solicitations and with the conditions they have grown [START_REF] Kretschmann | Ch. 5: Mechanical properties of Wood[END_REF]); thus, new robot design issues appear: How to be sure that an industrial robot made with wood can be accurate and stiff even if wood properties vary? A first attempt to introduce wood in industrial robot design was presented in [START_REF] Laurent | Vers un robot en bois: Première partie[END_REF]. The results showed that the approach was valid enough to compete with usual materials. However, this study did not deal with the aforementioned issues.

In this paper we introduce a design methodology leading to the manufacture of a stiff and accurate wooden five-bar mechanism. This methodology is based on two optimization problems that can be solved in cascade. First, control-based design approach is used. This approach allows to take into account the specificities of the controller used to servo the robot motions, and aims at finding the length of the links of the five-bar mechanism. Then the shape of the robot links using the reliable topology optimisation is developed. This problem is subject to a set of constraints related to the elastostatic and dynamic performance of the robot. Finally, the optimal length and shape of the wooden five-bar mechanism links are given and are used for the development of an prototype.

Design Process

The planar five-bar mechanism is shown in Fig. 1a. It is a two degrees-of-freedom (dof ) parallel robot able to achieve two translations in the plane (x 0 Oy 0 ). The endeffector is located at point C. The mechanism, is composed of five revolute joints, two of them located at point A 1 and A 2 are active while the others are passive. The links at points A i and B i (i = 1,2) are called the proximal links and between points B i and C the distal links. In this work, our objective is to design a wooden five-bar mechanism, the specifications that the robot should satisfy are summed up in Tab. 1. The robot should be as compact as possible due to some industrial constraints. Moreover, in order to ensure the accuracy of the robot even of the dimensions of the wooden links may vary with the atmospheric conditions, sensorbased controllers will be used. The most common approach consist of the direct In some cases however it may be difficult or unwise to observe the end-effector of the robot, e.g. in the case of a five-bar mechanism the end-effector is a point. For this reason we decided to define a controller based on the observation of the distal links, as in [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF]. However, as shown in [START_REF] Briot | The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques[END_REF] this controller is not free of singularities, that effect its accuracy performance which that not only depend on the camera localisation and specifications but also on the robot dimensions as shown in [START_REF] Kaci | Control-based design of a five-bar mechanism[END_REF]. Thus the optimization process must take all this parameters into account.

Link Length Optimization

The first design optimization problem aims at finding the optimal lengths of the robot links that minimize the footprint size of the five-bar mechanism in the plane of motion for a prescribed regular dexterous workspace, taking into account the controller accuracy performance. Four cameras with 1920×1200 pixels resolution are supposed to observe the two links attached to end-effector to control the robot end-effector pose. We introduced in Kaci et al. ( 2016) a procedure to define the error model between a single camera and the end-effector and these to define the robot optimal design variables. This method is adapted to the case of four cameras to optimally design a five-bar mechanism. Based on [START_REF] Kaci | Control-based design of a five-bar mechanism[END_REF], we solved (1)

Where:

x = [ 0 1 2 x c1 x c2 x c3 x c4 ] with 0 = A 1 A 2 , 1 = A 1 B 1 = A 2 B 2 and 2 = B 1 C = B 2 C
, and x ci (i = 1 to 4 ) is a vector of dimension 6 defining the position and orientation of the cameras. A = LH is the robot footprint (Figure 1a). W x and W y are the dimensions along x 0 and y 0 of the rectangular dexterous workspace [START_REF] Merlet | Parallel Robots[END_REF]. The robot must have a rectangular dexterous workspace of rectangular shape with dimensions along x 0 bigger than W x0 = 800 mm and along y 0 of W y0 = 200 mm in which the following properties are certified:

• no Type 2 singularities of the robot,

• the end-effector is within the image frame and thus the legs can be observed,

• knowing the resolution of the camera, the resolution of the end-effector position should be lower than 0.5 mm. • the static forces exerted into the passive joints are proportional to 1/ sin ξ , ξ being the angle between the distal links [START_REF] Briot | Investigation on the effort transmission in planar parallel manipulators[END_REF]). Consequently, it is decided that sin ξ should be higher than 0.1 to avoid excessive efforts in the joints. The optimal design variables of problem (1) are given in Tab. 2 and Tab. 3. The architecture of the optimized five-bar mechanism and the position of the dexterous regular workspace is illustrated in Fig. 1b. It should be mentioned that, in order to reduce the number of decision variables in the optimization problem, only the position/orientation of cameras 3 and 4 are optimized. Then, the position/orientation of the two other cameras is obtained by symmetry with respect to the y 0 axis.

Link Topology Optimization

As mentioned in Tab. 1, the deformation of the wooden five-bar mechanism should be lower then 0.5 mm under external forces. To reach the desired stiffness of the robot, topology optimization of the robot proximal links is used. Moreover, in order to decrease the dimensional variability of wood due to humidity, we decided to design our links with a novel type of wood named acetylated beech [START_REF] Hill | Wood Modification: Chemical, Thermal and Other Processes[END_REF]). This wood is known to have many advantages over raw wood species in terms of durability and dimensional stability. However, this acetylated beech, as any other type of wood, presents a normal distribution of its Young's modulus: here, the mean value is equal to 1277.2 MPa and the standard deviation is equal to 2043.5 MPa. To deal with this variation of the mechanical properties of acetylated beech, reliable topology optimization is developed. It should be mentioned that in what follows, only the shape of the robot proximal link are optimized. Indeed, the controller [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF]) impose that the distal links are cylindrical, for reason of facility of observation.

Modeling of the linkages elastic behavior Topology optimization uses the same physical model as in the FEM for modelling of a link, except that an interpolation scheme is used in order to define an artificial material. This method is called the Solid Isotropic Material with Penalization (SIMP, [START_REF] Bendsoe | Material interpolation schemes in topology optimization[END_REF]) and is known to be the most effective and the most widely used material interpolation scheme. This interpolation scheme is adopted in order to avoid having optimization results with too much intermediate material density, i.e. in order to have a final design solution with densities only equal to ρ i j = 1 (which represent the presence of materials)or ρ i j = 0 (which represent the absence of materials) without too many intermediate values (0 < ρ i j < 1) that are difficult to manage by the designer. The SIMP scheme is used to parameterize the Young's modulus associated with the stiffness matrix of the finite element i j and it is defined as follows:

E i j = E min + ρ p i j (E 0 -E min ), with ρ i j ∈ [0, 1] (2)
where E 0 is the Young's modulus of the material, E min is a very small stiffness value assigned to void regions in order to prevent the stiffness matrix from becoming singular, p (typically p = 3) is the penalization factor, and E i j is the Young's modulus of element j of the body i corresponding to the density variable ρ i j . Then, based on this definition of the Young's modulus for the element i j, we build its stiffness matrix. Finally, once all elementary matrices are defined, the computation of the body and linkage stiffness matrices is the same as in the traditional methodology [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]). Once the linkage stiffness is obtained, the elastic performance of the mechanism can be defined, in our case we will consider the deformation of the five-bar mechanism end-effector. For dealing with reliability of the robot performance due to the uncertainty of the wood mechanical properties, reliable topology optimisation approach must be used. We take advantages of the work [START_REF] Asadpoure | Robusttopology optimization of structures with uncertainties in stiffness application to truss structures[END_REF] to define the following constraint included in the design process:

E( u e ) + kσ ( u e ) ≤ u max (3)
Where: E(.) is the expectation operator, σ (.) is the standard deviation operator, k is a positive real and u e is the deformation vector at the node located at the center of the end-effector under a giving loading and u max is the maximum desired deformation. This constraint, is derivated from the well-known Bienaymé-Tchebichev theorem, which states that the probability, for a stochastic variable X, to have

| X -E(X) | ≥ k, is lower than σ (X) 2 k 2 .
Thus, by increasing the value of k in (3) we reduce the probability to have a robot whose deformation under the loading given in Tab. 1 will be bigger than 0.5 mm. For instance, with k = 3, the probability to have a robot whose deformations are acceptable is almost 90 %. Approximate formulas are provided in [START_REF] Asadpoure | Robusttopology optimization of structures with uncertainties in stiffness application to truss structures[END_REF] to compute E( u e ) and σ ( u e ) and their derivatives, which are necessary in topology optimization process.

Optimization Problem

The topology optimisation problem can be formalized as: min

ρ F = τ 2 (T * ) under g 1i = (δ 2 C (f i , q * ) -δ 2 max )/δ 2 max ≤ 0, i = 1, 2 g 2i = (θ 2 C (f i , q * ) -θ 2 max )/θ 2 max ≤ 0, i = 1, 2 (4)
Where q * are the selected exciting configurations for the computation of the translational/rotational deformation when the loading f i is applied (i=1,2). Configuration q * can be found using the methodology proposed in [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF]. The objective is to minimize the RMS of the actuators torques τ along optimized trajectory T * selected based on the approach shown in [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF]. Concerning the constraints, we must ensure that wherever in the dexterous workspace defined in section 2.1, the elastic translational displacement δ C of the end-effector under external forces f 1 = [50 N 0 N 1 Nm] T and f 2 = [0 N 50 N 1 Nm] T is lower than δ max = 0.5 mm and that the rotational displacement θ C of the distal links at the end-effector is lower than θ max = 0.25 mrad.

Here, δ C = E( u de ) + kσ ( u de ) and θ C = E( u re ) + kσ ( u re ), with u de is the transitional deformation of the node at the end-effector and u re is the rotational one. Indeed, we recall that in our case k = 3 in order to have a design reliability of 90 %.

Result of the topology optimization

The initial design domain for the proximal links is represented in Fig. 2a. Each link has two empty holes at the extremely for the joints insertion (motor shaft at points A i , and passive rotary joints for all others. For the meshing of the links, 24132 planar elements of size 1 × 1 mm are used and with 50 mm of thickness. The obtained shape of proximal links is shown Based on the previous optimisation results, a prototype of a wooden five-bar mechanism was designed (Fig. 3). Our future works will concern its visual servoing and the experimental validations of the obtained theoretical results.

Conclusion

Robot eco-design is a field of research in Robotics that should be deeper investigated: this could contribute to the Climate Change Mitigation. In order to considerably decrease the environmental impact of robot during the design phase, we proposed to replace metal or carbon composite parts by bio-sourced materials, such as wood. Wood has interesting mechanical properties but also drawbacks: its mechanical performance encounters a large variability. In order to be able to design a stiff industrial robot, robust design approaches must be used. In this paper, we developed a wooden parallel robot, two optimization problems have been formulated and solved in cascade in order to find the robot optimal design architecture for which the impact of the wood properties variability is minimal. Finally, the CAD design and prototyping of the wooden five-bar mechanism and its key technological solutions have been described.

Our future works, will concern the implementation of sensor-based controllers. Then the experimental validations of the theoretical results obtained in the framework of the RobEcolo project will be made soon to characterize the robot properties.

  Figure 1: Five-bar mechanism and the optimized one

Figure 2 :Figure 3 :

 23 Figure 2: Evolution of the design of the five-bar links: the links are shown in grayscale (black elements correspond to ρ i = 1, white elements to ρ i = 0, and gray elements to 0 < ρ i < 1)

Table 2 :

 2 Optimal Design Solution of the Wooden Five-bar Mechanism

	A [m 2 ]	l 0 [m]	l 1 [m]	l 2 [m]
	0.1372	0.125	0.280	0.400

Table 3 :

 3 Optimal Position and Orientation of Cameras W x0 and W y > W y0

	Pos/Ori	x c [m]	y c [m]	z c [m]	φ [rad]	θ [rad]	ψ [rad]
	Camera3	0.01	0.2291	0.75	0	π	0
	Camera4	0.02	0.100	0.100	0	π	0
	the following optimization problem:
		minimize	A = LH
			over	x
		subject to	W x >

Acknowledgement

This work was supported by the project RobEcolo funded by the French Région Pays de la Loire (Convention No. 2015-10773).