S Briot

P Martinet

F Chaumette
email: francois.chaumette@inria.fr

Determining the Singularities for the Observation of Three Image Lines

Sébastien Briot1 , Philippe Martinet 2 and Franc ¸ois Chaumette 3 Abstract-The determination of the singularity cases in the observation of image features is a complicated problem which is still open, apart from image points. For the first time, we provide the singularity cases in the observation of three image lines. We show that a concept named the "hidden robot", which was formerly used for understanding the singularities of a vision-based controller dedicated to parallel robots and, more recently which proved to be efficient for finding the singularity cases in the observation of image points, can be used for interpreting the singularities in the observation of three image lines. The hidden robot concept considerably simplifies the analysis by using geometric interpretations of the mapping degeneracy and tools provided by the mechanical engineering community. We prove that in the most complicated case where three general lines in space are observed, singularities appear when the origin of the observed object frame is either on a quadric or a cubic surface whose location in space depends on the configuration of the observed lines. In simpler cases where at least two lines belong to the same plane, these two surfaces can degenerate into simpler geometrical loci (e.g. planes, cylinders, lines).

I. INTRODUCTION

Finding the singularity cases in the observation of image features is crucial in visual servoing and pose estimation [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]. However, this problem is still unsolved for most of the cases due to its considerable complexity. The singularities were found in the case of three image points after rather complicated mathematical computations [START_REF] Michel | Singularities in the determination of the situation of a robot effector from the perspective view of 3 points[END_REF], but they are still unknown for other image features.

In order to avoid the controllability issues near singularities, authors usually propose to observe additional features. However this leads to a non-minimal representation of the interaction between the system and the camera, resulting in the potential appearance of local minima [START_REF] Chaumette | Visual Servoing and Visual Tracking[END_REF] which are also difficult to determine. Additionally, even adding visual features may not ensure the absence of singularity cases [START_REF] Rosenzveig | A method for simplifying the analysis of leg-based visual servoing of parallel robots[END_REF]. Therefore, finding the singularities is crucial but it is usually prevented by the complexity of the mathematical problem to solve.

Recently, a concept named the "hidden robot" was introduced in [START_REF] Briot | Minimal representation for the control of Gough-Stewart platforms via leg observation considering a hidden robot model[END_REF], [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 concept was first used to determine the singularity cases of a vision-based controller dedicated to parallel robots [START_REF] Andreff | Vision-based control of a Gough-Stewart parallel mechanism using legs observation[END_REF] but it was recently proven that it can be extended to more general cases, such as the observation of n image points [START_REF] Briot | Revisiting the determination of the singularity cases in the visual servoing of image points through the concept of "hidden robot[END_REF].

The basic idea 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], [START_REF] Briot | Revisiting the determination of the singularity cases in the visual servoing of image points through the concept of "hidden robot[END_REF] is that the singularity cases of the mappings involved in the visual observation of geometric primitives for estimating the pose of an object or for visual servoing are equivalent to mappings representing the geometric and kinematic properties of given hidden parallel robots. By geometric property, we mean that the solutions of the Forward Geometric Model (FGM) of the parallel robots under consideration are also the solutions of the 3-D localization problem associated with the considered observations. By kinematic property, we mean that the singularities of the inverse kinematic Jacobian matrix of the robot are the same as the singularities of the interaction matrix that defines the link between the time variation of the visual features and the camera velocity. The hidden robot is indeed a tangible visualization of the mapping between the observation space and the Cartesian space. A methodology is proposed 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] in order to define the hidden robot models associated with visual servoing dedicated to parallel robots while [START_REF] Briot | Revisiting the determination of the singularity cases in the visual servoing of image points through the concept of "hidden robot[END_REF] presents an approach for finding the virtual robot models associated with more general cases based on the observation of geometric visual primitives.

By finding this correlation, it is then possible to study the singularities of the interaction matrix, by using advanced tools coming from the mechanical engineering community (e.g. the Grassmann-Cayley algebra [START_REF] Ben-Horin | Singularity analysis of a class of parallel robots based on Grassmann-Cayley algebra[END_REF] and/or the Grassmann geometry [START_REF] Merlet | Parallel Robots[END_REF]). The interest in using these tools is that they are (most of the time) able to provide simple geometric interpretations of the singularity cases. However, it should be noted that these tools still require an experienced user.

In the present paper, we show that a virtual parallel robot is hidden within the mapping used in the observation of image lines [START_REF] Andreff | Visual servoing from lines[END_REF] and, thanks to this correlation, we analyze the singularity cases when three lines are observed. To the best of our knowledge, this is the first time that this problem is solved, except for obvious cases such as when the three 3-D lines are parallel (leading to an unconstrained translation along the direction of the lines).

As a result, the paper is organized as follows. Next section shows that the geometric / kinematic mapping involved in the observation of image lines is the same as the mapping required to control a particular 3-UPRC parallel robot 1 . Then, in Section III, the singularities of the mapping are analyzed and validated through simulations. Finally, conclusions are drawn in Section IV.

II. THE ROBOT MODEL HIDDEN IN THE OBSERVATION OF LINES

Before presenting the architecture of robot hidden in the observation of three lines, we make some brief recalls on the computation of the related interaction matrix.

A. Interaction matrix of three image lines

In the following section of the paper, we use the standard pin hole model with focal length equal to 1 for the representation of the camera model. However, any other model based on projective geometry could be used.

A 3-D line L i of Plücker coordinates [U T i L T i ] T in the camera frame (U i being a unit vector characterizing its direction while L i = U i ×X, X being the coordinates of any point M i belonging to L i ) is projected in the image plane on a 2-D line i of Plücker coordinates [u T i l T i ] T (u i being a unit vector characterizing its direction while l i = u i × x, x being the coordinates of any point m i belonging to i ) (Fig. 1). From [START_REF] Chaumette | La relation vision-commande : théorie et application à des tâches robotiques[END_REF], we know that

u i =   u xi u yi u zi   =   L yi /∆ -L xi /∆ 0   , l i =   l xi l yi l zi   =   L xi /∆ L yi /∆ L zi /∆   (1)
where

L i = [L xi L yi L zi ] T and ∆ = L 2 xi + L 2 yi .
By differentiating these equations, the classical equations linking the velocities ˙li (and as a result the velocities ui ) to the twist τ T c = [υ T c ω T c ] of the camera in its relative motion with respect to the observed object frame (υ c being the translational velocity and ω c the rotational velocity expressed in the camera frame) are

˙li = M i τ c (2)
where

M i = [M T i1 M T i2 M T i3 ] T with M i1 =         l xi l yi U zi /∆ l 2 yi U zi /∆ -l yi (l xi U xi + l yi U yi )/∆ l xi l yi l zi l 2 yi l zi -l yi         T (3)
M i2 =         -l 2 xi U zi /∆ -l xi l yi U zi /∆ l xi (l xi U xi + l yi U yi )/∆ -l 2 xi l zi -l xi l yi l zi l xi         T (4)
M i3 =         (U yi + l yi l zi U zi )/∆ -(U xi + l xi l zi U zi )/∆ l zi (l xi U yi -l yi U xi )/∆ l yi (l 2 zi + 1) -l xi (l 2 zi + 1) 0         T (5)
In order to fully servo the relative motion between an object and a camera, at least three lines fixed on the object must be observed [START_REF] Andreff | Visual servoing from lines[END_REF] (Fig. 2). Thus, considering the observation of three lines L 1 , L 2 and L 3 , the interaction matrix linking the velocities ˙li of the lines i (i = 1, 2, 3) grouped in the vector ṡ = [ ˙lT

1 ˙lT 2 ˙lT 3 ] T to the camera twist τ c by the relation ṡ = Mτ c (6)
is thus given by

M = M T 1 M T 2 M T 3 T (7)
Singularities appear when the matrix M is rank deficient. Looking at the analytical form of M, it is clear that determining the singularities of M through the determinant of M T M seems to be out of reach. Note that we have used above the Plücker representation of lines, but following results are not dependent of the choice of the representation (Cartesian, cylindrical) [START_REF] Chaumette | La relation vision-commande : théorie et application à des tâches robotiques[END_REF].

B. Hidden robot model

In this section, we show that a virtual parallel robot is hidden within the mapping used in the observation of image lines by considering what follows. First, from the single measure of the location of a line i in the image plane, it is impossible to know the location of its corresponding 3-D line L i . The only information that we can extract from

P i L i α i β i x i z i ξ i1 ξ i2 a i , U i y i y i , L i x i i γ Fig. 3. A U P RC kinematic chain
this measure is that 3-D line L i lies on the interpretation plane P i passing through i and the optical center O. As a result, the measure of i gives the location of the plane P i . Furthermore, from the same single measure, the orientation of the observed body B cannot be deduced (we need the three measures of the positions of lines 1 , 2 and 3 , see Fig. 2).

From a mechanical engineer point of view, these geometric properties can be obtained by the kinematic architecture depicted in Fig. 3. This architecture is made of an actuated cardan (or universal (U )) joint rotating around point O and fixed on the camera frame at that point. The U joint is followed by a passive prismatic (P ) joint whose direction z i is reciprocal to the axes x i and y i of the cardan joint. Then, the passive P joint is attached at its other extremity to a passive revolute (R) joint whose axis is directed along y i . Finally, the last link (moving platform, equivalent to the observed body B) is connected to the R joint by a passive cylindrical (C) joint whose axis a i is orthogonal to the axis of the R joint. Thus, we have a leg with a U P RC architecture linking the camera frame to the observed object frame. To this leg, we associate a vector q i = [α i β i ] T , where α i and β i represent the rotation angles of each revolute joints composing the U joint (Fig. 3). With such a leg, the observed body B is constrained to translate along L i and to rotate around it. As a result, the vectors U i and L i associated with the Plücker coordinates of the line L i are collinear with the axes a i and y i , respectively (Fig. 3).

It should be mentioned that the rotational component of the C joint in this leg is necessary because we know that line L i is attached to a rigid element (the observed body B) whose orientation cannot be defined by considering a single measure i . Without taking into account this information, the displacement of the line L i can be obtained by a leg with a U P R architecture, i.e. a kinematic chain with two passive degrees of freedom (instead of four for the U P RC leg).

Considering now that:

• the three observed 3-D lines L 1 , L 2 and L 3 have mo- tions whose geometric properties can be parameterized, for each of them, by a U P RC mechanical architecture depicted above, • L 1 , L 2 and L 3 are fixed on the same body B,

• there is a global diffeomorphism between the configuration of the line i and the U joint rotations q i ; as a result, l i is a singularity-free observation of the U joint motions q i , then the relative motion between the body B and the camera frame has the same geometric and kinematic properties as the motion of a 3-U P RC parallel robot (Fig. 4), i.e.

• the solutions of the FGM of the 3-U P RC parallel robot are also the solutions of the 3-D localization problem when three 3-D lines are observed by a perspective camera.

• the singularities of the inverse kinematic Jacobian matrix of the robot are the same as the singularities of the interaction matrix (7). To prove the last item, let us consider what follows. The inverse kinematic Jacobian matrix J inv of the robot links the actuators velocities q = [ qT

1 qT 2 qT
3 ] T to the platform twist τ p through the relation

q = J inv τ p ( 8
)
J inv is of dimension (6×6). Moreover, as mentioned above, there is a global diffeomorphism between the measure of the image line configuration i and the U joint rotations q i . As a result, q = Gṡ (9

)
where G is a (6 × 9) block-diagonal matrix which cannot be rank-deficient. Finally, by combining ( 6) and ( 9) with the screw transformation relation linking the camera twist τ c

(expressed with respect to the object frame) to the observed object twist τ p (expressed with respect to the camera frame) by [START_REF] Chaumette | Visual Servoing and Visual Tracking[END_REF]

τ c = H τ p , where H = -I 3 [ c t p ] × 0 3×3 -I 3 (10)
(in which I 3 is the (3 × 3) identity matrix, 0 3×3 is a (3 × 3) zero matrix, and [ c t p ] × is the cross-product matrix associated to the 3-D vector c t p characterizing the translation between the camera and the observed object frames), we found:

J inv = GMH (11)
Then, as matrices G and H are never rank-deficient, we obligatorily have a loss of rank of M if and only if J inv is rank-deficient. Therefore, the conditions of singularity of the hidden robot inverse kinematic Jacobian matrix J inv are the same as the singularity conditions of the interaction matrix M.

III. SINGULARITY ANALYSIS A. General conditions of singularity

Singularities on the inverse Jacobian matrix J inv of a parallel robot (also called Type 2 or parallel singularities [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF]) appear when at least two solutions of the FGM are identical [START_REF] Merlet | Parallel Robots[END_REF]. As mentioned in Section II, these singularities are analogous to the singularities of the interaction matrix.

In Type 2 singularities, parallel robots gain one (or more) uncontrollable motion. Kinematically speaking, there exists a non-null vector t s defined such that J inv t s = 0 while q = 0, i.e. the actuators are fixed (which means that t s is in the null space of J inv ). As known in mechanics, if a rigid body got an uncontrollable motion, this means that it is not fully constrained by the system of wrenches applied on it, i.e. the static equilibrium is not ensured. As this uncontrollable motion appears only in a singularity, this means that locally the system of actuation wrenches, i.e. wrenches transmitted from the actuators to the platform by the legs, is degenerated [START_REF] Merlet | Parallel Robots[END_REF].

For a given leg i, any actuation wrench denoted by ξ ij is reciprocal to the unit twists denoted ζ ik characterizing the displacements of the passive joints [START_REF] Bonev | Geometric analysis of parallel mechanisms[END_REF], i.e. ξ T ij ζ ik = 0 for any j and k. This means that the virtual power developed by the wrench ξ ij along the direction of motion ζ ik is null; in other words, the actuator j of the leg i cannot transmit a wrench ξ ij to the platform along the direction ζ ik .

Let us consider a U P RC leg belonging to our 3-U P RC hidden robot. In the frame F i : (M i , x i , y i , z i ) (M i being the point of intersection between the axis of the prismatic and revolute joints) attached to the leg, the unit twist defining the motion of the passive P joint is expressed as (see Fig. 3):

ζ i1 = 0 0 1 0 0 0 T (12) while ζ i2 = 0 0 0 0 1 0 T ( 13
)
is the unit twist defining the motion of the passive R joint and

ζ i3 = cos γ i 0 -sin γ i 0 0 0 T ( 14
)
ζ i4 = 0 0 0 cos γ i 0 -sin γ i T (15)
are the unit twists defining the motions of the passive C joint, γ i being the angle between x i and a i (and thus U i ), both axes being contained in P i (Fig. 3). In these twists, the three first components represent the direction of the translation velocity while the three last components represent the direction of the rotational velocity.

As a result, the unit actuation wrenches

ξ ij = [f T ij m T ij ] T expressed in the frame F i are ξ i1 : f i1 = [0 1 0] T , m i1 = 0 3×1 (16)
ξ i2 : f i2 = 0 3×1 , m i2 = [sin γ i 0 cos γ i ] T (17)
in which f ij represents the direction of the force exerted on the platform and m ij the direction of the moment. As a result, ξ i1 is a pure force directed along y i and ξ i2 is a pure moment reciprocal to ξ i1 and U i which is included in the plane P i (Fig. 3). Thus, for the three robot legs, the system of actuation wrenches is given by ξ = [ξ 11 ξ 21 ξ 31 ξ 12 ξ 22 ξ 32 ]. There exists some tools able to define the conditions of degeneracy of a wrench system among which are the Grassmann geometry [START_REF] Merlet | Parallel Robots[END_REF] and the Grassmann-Cayley algebra [START_REF] Ben-Horin | Singularity analysis of a class of parallel robots based on Grassmann-Cayley algebra[END_REF], [START_REF] Ben-Horin | Application of Grassmann-Cayley algebra to geometrical interpretation of parallel robot singularities[END_REF]- [START_REF] Kanaan | Singularity analysis of lower mobility parallel manipulators using Grassmann-Cayley algebra[END_REF]. In what follows, we use the Grassmann-Cayley algebra in order to find the singularity conditions related to our problem. Indeed, as any wrench of the wrench system ξ can be seen as the Plücker representation of a line2 [START_REF] Merlet | Parallel Robots[END_REF], the Grassmann-Cayley algebra makes it possible to compute the determinant of the wrench system ξ as an expression involving twelve points selected on the six lines (corresponding to the six wrenches). This expression is a linear combination of 24 monomials [START_REF] Ben-Horin | Singularity analysis of a class of parallel robots based on Grassmann-Cayley algebra[END_REF], each monomial representing the volume of a tetrahedron whose extremities correspond to four of the considered twelve points. By a smart selection of the twelve points, due to the experience of the user, it is possible to vanish a large number of the monomials and thus to simplify the expression of the determinant of the wrench system.

Regarding our particular case, the Grassmann-Cayley algebra was used in [START_REF] Kanaan | Singularity analysis of lower mobility parallel manipulators using Grassmann-Cayley algebra[END_REF] to prove that, if the system of wrenches is composed of a pair of three forces ξ 11 , ξ 21 , ξ 31 and of three moments ξ 12 , ξ 22 , ξ 32 , conditions of singularities appear if and only if :

f 1 = f T 11 (f 21 ×f 31 ) = 0 or f 2 = m T 12
(m 22 ×m 32 ) = 0 (18) which means that the vectors f 11 , f 21 and f 31 (m 11 , m 21 and m 31 , resp.) lie in the same plane (or are collinear to the same vector c). As a result, the mechanism is not able to resist to forces f (moments m, resp.) orthogonal to this plane (or in any directions orthogonal to c) and thus translations (rotations, resp.) are gained along the direction of f (around the axis of m, resp.).

From these two conditions, it is possible to find the configurations of the end-effector (the observed body) in the camera frame leading to singularities. Examples of singular configurations, depending on the arrangement of the lines on the observed body, are detailed below.

B. Singularity loci

Thanks to the particular geometric properties of the system to analyze (invariance of robot leg configurations for any rotation around O), it is possible to simplify the singularity analysis by fixing the platform orientation. Thus, all expressions below will be given for the "zero" platform orientation. For another given platform orientation (parameterized by the rotation matrix R), singularity loci will be found by parameterizing the variables X, Y and Z characterizing the position of the origin of the object frame F b in the camera frame (see Figs. 5 and6) when considering the "zero" platform orientation thanks to new variables X , Y and Z representing the position of the origin of the object frame for the considered "non-zero" platform orientation such that

X Y Z T = R X Y Z T (19)
1) Three coplanar lines with no common intersection point: In such a case, the three lines intersect in points M 1 , M 2 and M 3 (Fig. 5(a)) (we disregard the case where at least two lines are parallel, which is considered in Section III-B.5). We define the frame F b : (Q, x b , y b , z b ) attached to the observed body B such that Q is the center of the circumcircle of the triangle ∆M 1 M 2 M 3 and x b is collinear to ---→ QM 1 . Moreover, for the analysis, we fix the orientation of the body so that F b can be obtained from the camera frame by a translation of vector --→ OQ. As a result,

--→ OQ = [X Y Z] T , ---→ QM 1 = [ρ 0 0] T , ---→ QM 2 = ρ[cos φ sin φ 0] T , ---→ QM 3 = ρ[cos ψ sin ψ 0] T (20)
where φ, ψ are two angles defined in Fig. 5(a) and ρ is the radius of the circle. Then, from section III-A, we know that

f i1 ∝ U i × --→ OM i , m i2 ∝ U i × -→ f i1 (21)
leading to, from (18), and after some tedious developments,

f 1 = 0⇔Z = 0 f 2 = 0⇔Z(X 2 + Y 2 -ρ 2 ) = 0 (22)
Condition Z = 0 means that the plane containing the lines L 1 , L 2 and L 3 also contains the camera center O, which is not critical in practice. In such a case, all three lines project in a single line in the image plane and the pose of the object cannot be determined. Condition X 2 + Y 2 = ρ 2 means that the camera center O lies on the cylinder whose axis is perpendicular to the plane containing all three points M 1 , M 2 and M 3 and which includes the three points. These conditions are the same as the singularity cases where three image points are observed [START_REF] Michel | Singularities in the determination of the situation of a robot effector from the perspective view of 3 points[END_REF], which is indeed not a surprise due to the equivalence of both problems.

1 M 2 M 3 M b x b y 1 U 2 U 3 U 1 L 3 L 2 L ) X Y Z ( Q φ ψ ρ (a) with no common intersection point b x b y ) X Y Z ( Q 1 U 2 U 1 L 3 L 2 L φ ψ
(b) with a common intersection point Fig. 5. Observation of three coplanar lines 2) Three coplanar lines with a common intersection point: We define the frame F b : (Q, x b , y b , z b ) attached to the observed body B such that Q is the intersection point of all three lines and x b is collinear to U 1 (Fig. 5(b)). Again, we fix the orientation of the body so that F b can be obtained from the camera frame by a translation of vector --→ OQ. Thus,

--→ OQ = [X Y Z] T , U 1 = [1 0 0] T , U 2 = [cos φ sin φ 0] T , U 3 = [cos ψ sin ψ 0] T (23)
where φ, ψ are two angles defined in Fig. 5(b).

Then, from section III-A,

f i1 ∝ U i × --→ OQ, m i2 ∝ U i × -→ f i1 (24)
leading to, from (18),

f 1 = 0 for any object configuration f 2 = 0⇔Z(X 2 + Y 2 ) = 0 (25)
Thus, if three coplanar lines intersect in a single point, the rank of the interaction matrix is never full, i.e. at least one degree of freedom (a translation towards the intersection point) of the object cannot be controlled, as already known [START_REF] Andreff | Visual servoing from lines[END_REF]. Moreover, additional degrees of freedom become uncontrollable if

• Z = 0, which means as previously that the plane containing the lines L 1 , L 2 and L 3 also contains the camera center O. • X 2 + Y 2 = 0, which means that the camera center O lies on the line which passes through Q and which is perpendicular to all vectors U i , which is the same condition as when f 1 = 0. 3) Three lines in space with a common intersection point: Here, the three lines are spatial, not orthogonal and intersect in a common point. We define the frame F b : (Q, x b , y b , z b ) attached to the observed body B such that x b is collinear to U 1 , y b is in the plane containing U 1 and U 2 , while Q is the intersection point of the lines. As a result,

--→ OQ = [X Y Z] T , U 1 = [1 0 0] T , U 2 = [a b 0] T , U 3 = [c d e] T (26)
where a, b, c, d and e are variables parameterizing the direction of the lines.

) X Y Z ( Q 1 U 2 U 3 U 1 L 3 L 2 L (a) no common intersection point ) X Y Z ( Q b x b y 1 U 2 U 3 U b z a b 1 L 3 L 2 L
(b) two lines are parallel Fig. 6. Observation of lines in space f i1 and m i2 can be found thanks to expression (24), and from (18) we obtain, f 1 = 0 for any object configuration

f 2 = 0⇔b(adeY 3 + ((-ad 2 + bcd + ae 2 )Z + (ac -bd)eX)Y 2 -e(bcX 2 + (ad -bc)Z 2 + 2beXZ)Y + ((-ad 2 + bcd -ae 2 )X 2 Z + (bd + ac)eXZ 2 )) = 0 (27)
Thus, if three spatial lines intersect in a single point, the rank of the interaction matrix is never full, i.e. at least one degree of freedom (a translation) of the object cannot be controlled, as already known [START_REF] Andreff | Visual servoing from lines[END_REF]. Moreover, additional degrees of freedom become uncontrollable if the origin of the body frame belongs to a cubic surface parameterized by f 2 = 0.

Note also that in the case where the three lines are orthogonal, the equation f 2 = 0 degenerates into XY Z = 0.

4) Three orthogonal lines in space: Here, the three lines are orthogonal and do not intersect. We define the frame F b : (Q, x b , y b , z b ) attached to the observed body B such that x b is collinear to U 1 , y b is collinear to U 2 and z b is collinear to U 3 . Moreover, the distance between L 1 and L 2 is equal to 2a, between L 2 and L 3 is equal to 2b, and between L 1 and L 3 is equal to 2c, while Q is the barycenter of the lines (Fig. 6(a)). As a result,

--→ OP 1 = [X (Y -c) (Z -a)] T , U 1 = [1 0 0] T --→ OP 2 = [(X -b) Y (Z + a)] T , U 2 = [0 1 0] T --→ OP 3 = [(X + b) (Y + c) Z] T , U 3 = [0 0 1] T (28)
where P i is a point belonging to L i . Then, we have

f i1 ∝ U i × --→ OP i , m i2 ∝ U i × -→ f i1 (29)
which leads to, from (18),

f 1 = 0⇔aXY + bY Z -cXZ -abc = 0 f 2 = 0⇔acX -abY + bcZ -XY Z = 0 (30)
Expression f 1 represents a quadric surface while expression f 2 is a cubic surface (Fig. 7).

Z [m]    - - X [m]    - - Y [m]    - - = 0 1 f = 0 1 f = 0
--→ OP 1 = [X (Y -a) Z] T , U 1 = [1 0 0] T --→ OP 2 = [X (Y + b) Z] T , U 2 = [1 0 0] T --→ OP 3 = [X Y Z] T , U 3 = [c d e] T (31)
where c, d, e are variables parameterizing the direction of the line L 3 .

Then, f i1 and m i2 being found from (29), by using the singularity conditions (18), we obtain

f 1 = 0⇔Z(dZ -eY ) = 0 f 2 = 0⇔Z(X(d 2 + e 2 ) -cY d -cZe) = 0 (32)
Here, singularities are equations of planes:

• Z = 0, which occur when the plane P containing L 1 and L 2 also contains the optical center, • eY -dZ = 0 is the plane containing U 1 , U 3 and the optical center, • X(d 2 + e 2 ) -cdY -ceZ = 0 is the plane containing (U 1 × U 3 ), U 3 and the optical center.

The two first cases are degenerated cases where two 3-D lines are contained in the same interpretation plane (i.e. the plane passing through the camera center and the image line), leading to the same image line for the two 3-D lines. 6) Three general lines in space: This is the most general case. We do not provide here the expressions which are more complex than in the previous cases but they are still exploitable and can be found in a technical report [18].

As in the case of the three perpendicular lines, condition f 1 = 0 provides the expression of a quadric surface while f 2 = 0 leads to a cubic surface. Inverse of the condition number of the interaction matrix M obtained when observing three lines in the general case.

C. Illustrative examples

1) Singularities when observing three general lines in space: In order to show the exactness of our results, we first perform simulations in the case where the three lines have a general configuration (Section III-B.6). The three general lines are parameterized as follows. We define the frame F b : (Q, x b , y b , z b ) attached to the observed body B such that x b is collinear to U 1 , y b is lying in the plane P containing L 1 and L 2 . Then, the line L i is parameterized by its direction U i and a point P i lying on it by taking the following expressions:

--→ OP 1 = [(X -b) (Y -c) (Z -a)] T , U 1 = [1 0 0] T --→ OP 2 = [(X -b) (Y -c) (Z + a)] T , U 2 = [d e 0] T --→ OP 3 = [(X + b) (Y + c) Z] T , U 3 = [f g h] T (33)
where d, e, f , g and h are variables parameterizing the direction of the lines L 2 and L 3 , and X, Y , Z are the coordinates of Q in the camera frame. For simulation purpose, we set a = b = c = 1 m and d = e = 1, f = -1, g = h = 1. Then, a relative motion between the object frame and camera frame origins is imposed and is parameterized by the following functions characterizing the location of Q in the camera 8) which is a point lying on the quadric surface whose expression can be found in the technical report [18]. As a result, the rank-deficiency of the matrix M appears for s = 0, which is shown on Fig. 9 in which the inverse of the condition number is null only for s = 0.

= -1 ⇒ 1 L 2 L 3 L X O [m] Y O [m] Z O [m]
2) Singularities when observing three orthogonal lines in space: A second simulation is performed in order to show now the singularities when observing three orthogonal lines in space (see Section III-B.4). For simulation purpose, we take the following values of the parameters a, b and c in (28): a = b = 0.5 m, c = 1.5 m. Then, a relative motion between the object frame and camera frame origins is imposed and is parameterized by the following functions characterizing the location of Q in the camera frame: X = 1-0.5s, Y = 1 -0.5s,

Z = 2 -s (35)
with s ∈ [-1 1] a linearly increasing function. The relative motion of the camera with respect to the observed lines is shown in Fig. 10. In such simulation, when s = 0, the point O reaches the configuration X O = 0.5 m, Y O = 0.5 m and Z O = 1 m in the object frame (see Fig. 10) which is a point lying on the cubic surface whose expression is given in (30). As a result, the rank-deficiency of the matrix M appears for s = 0, which is shown on Fig. 11 in which the inverse of the condition number is null only for s = 0.

IV. CONCLUSION In this paper, we determined the singularity cases for the observation of three image lines thanks to a tool named the "hidden robot concept". We showed that the hidden robot concept allows for considerable simplification of the analysis, leading to the computation of new singularity cases. Indeed, the hidden robot is a tangible visualization of the mapping between the observation space and the Cartesian space. As a result, the singular configurations of the hidden robot corresponded to the singularities of the interaction matrix.

Indeed, the concept of hidden robot allowed to change the way we defined the problem. By doing so, we were able to replace the degeneracy analysis of the velocity transmission between inputs (velocity of the observed features) and outputs (camera twist), by its dual but fully equivalent problem which was to analyze the degeneracy in the transmission of wrenches between the inputs of a virtual mechanical system (virtual actuators of the hidden robot whose displacement was linked to the motions of the observed features) and its outputs (wrenches exerted on the virtual platform, i.e. the observed object).

Then, by using geometric interpretations of the mapping degeneracy and tools provided by the mechanical engineering community such as the Grassmann-Cayley algebra, we were able to find rather simple geometric interpretation of the interaction matrix degeneracy. We proved that in the most complicated cases where three general lines in space are observed, singularities appear when the origin of the observed object frame is either on a quadric or a cubic surface. In simpler cases where at least two lines belong to the same plane, these two surfaces can degenerate into simpler geometrical loci (e.g. planes, cylinders, lines).

Future works concern the analysis of the singularity cases in the observation of n lines, and research about the singularities when observing other primitives (spheres, circles, and eventually of combination of different primitives), or when using several cameras.

Fig. 1 .

 1 Fig. 1. Observation of a line

Fig. 2 .

 2 Fig. 2. Observation of three lines lying on a body B

2 L 3 LFig. 4 .

 234 Fig.4. The hidden robot model corresponding to the observation of three lines lying on a body B: a 3-U P RC parallel robot with all active cardan joints merged at the point O (for reason of clarity of drawings, the axes of the cardan joints are not represented)

2 fFig. 7 .

 27 Fig. 7. Example of singularity loci when three perpendicular lines are observed (computed for a = 3.13, b = 3.13, c = -3.13)

Fig. 8 .

 8 Fig. 8. Relative motion of the camera with respect to the observed lines in a general configuration.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Relative motion of the camera with respect to the three observed orthogonal lines.

In the following of the paper, P, U, R, C will stand for passive prismatic, universal, revolute and cylindrical joints, respectively. If the letter is underlined, the joint is considered active.

A pure force will be represented by a line in the

3-D space while a pure moment will be a line in the projective space at infinity.

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


Damien Six
email: damien.six@irccyn.ec-nantes.fr

Abdelhamid Chriette
email: abdelhamid.chriette@irccyn.ec-nantes.fr

Sébastien Briot
email: sebastien.briot@irccyn.ec-nantes.fr

Philippe Martinet
email: philippe.martinet@irccyn.ec-nantes.fr

Dynamic Modeling and Trajectory Tracking Controller of a Novel Flying Parallel Robot

Keywords: Flying robots, Nonlinear Control, Feedback linearization, Dynamics, Parallel robots

This paper presents a new flying robot, composed of two quadrotors linked by a rigid articulated passive chain. The robot obtained is similar to a parallel robot where the classic actuators have been replaced by flying drones. With its rigid structure presenting an internal degree of freedom, the robot presented is a step forward in flying robotics and presents new challenges for the design of its feedback control. In this paper the dynamic modelling of the robot is analyzed. From this analysis, a decoupling property is extracted leading to a cascaded controller based on feedback linearization. An application to trajectory tracking is developed. The effectiveness and robustness against noise in pose estimation of the proposed controller is verified through a realistic simulation.

INTRODUCTION

The past decade has seen the development of quadrotors Unmanned Aerial Vehicle (UAV) with high performances in constrained environments. They have gained increasing interest in a large panel of applications as observation, surveillance and logistics. Several studies have also explored drone aerial manipulation, through the use of rigid tool attached to the drone (Nguyen et al., 2015a) or serial manipulator associated to one drone as in [START_REF] Yang | Dynamics and Control of Quadrotor with Robotic Manipulator[END_REF]; [START_REF] Korpela | Dynamic stability of a mobile manipulating unmanned aerial vehicle[END_REF]. Such robots have low payload and, when serial manipulators are used, the additional motors of the arm reduce the energetic autonomy. To overcome those limitations, some studies have explored the collaborative work of several drones associated to a rigid body [START_REF] Nguyen | Aerial tool operation system using quadrotors as Rotating Thrust Generators[END_REF][START_REF] Kim | Visionbased collaborative lifting using quadrotor UAVs[END_REF] or cable suspended payload transportation [START_REF] Manubens | Motion Planning for 6-D Manipulation with Aerial Towed-cable Systems[END_REF][START_REF] Sreenath | Dynamics, Control and Planning for Cooperative Manipulation of Payloads Suspended by Cables from Multiple Quadrotor Robots[END_REF]. This last configuration can be seen as a cable-driven parallel robot where classic actuators have been replaced by quadrotors. The cables are an appropriate solution to lift a load but, due to the intrinsic characteristics of cables, such structure cannot apply a pushing force on the environment. This paper explores a new flying structure inspired from rigid parallel robots. Several drones are linked with a rigid articulated passive chain. The flying robot obtained present several advantages over the other explored solutions:

• the ability to perform tasks under but also over the drones, depending of the passive chain configuration, • the efforts applied on the tool are spread over the drones, enhancing the total payload of the robot, • the absence of additional embedded motors to actuate the effector reduces the load of the system itself and maintain the energetic autonomy of the drones, • the large choice of leg topology [START_REF] Gogu | Structural synthesis of parallel robots[END_REF] leading to a variety of physical properties of potential interest.

Those properties offer interesting potential future applications in maintenance, manipulation and inspection in constrained environments. The paper is structured as follows.

A kinematic description of the robot is presented in section 2. In the design of the flying robot, the authors selected quadrotors to lift the structure and perform the actuation (see Fig. 1). This solution was taken on due to the increasing availability of low cost and reliable quadrotors.

To the knowledge of the authors, such flying structures have not been already explored and offer new challenges in the design of a feedback control law. In section 3, an analysis of the dynamic model is then performed in order to highlight the possible decouplings to control the underactuated structure. From this analysis, a controller based on cascaded feedback linearization is developed in section 4. The stability of the designed controller is discussed. Simulations are then performed to check the controller performances in section 5.

PARAMETRIZATION AND KINEMATIC MODEL

The flying structure is composed of quadrotors, linked by a passive kinematic chain (Fig. 1). The underactuation To reduce the complexity of the study, the robot is designed to perform motions restricted to a vertical plane. Precautions taken in the control law to maintain the flying structure in this plane are discussed in section 4. The flying robot is then derived from a planar parallel mechanism, the biglide mechanism. This flying robot obtained is composed of two drones (named drone 1 and drone 2) and a two links passive kinematic chain (B 1 and B 2 ). Revolute joints are selected to constrain the motion of the passive kinematic chain to a planar motion.

The passive joints connecting the passive chain to the drone 1 and 2 are respectively named joint x and y. The passive joint corresponding to an end-effector in the middle of the passive kinematic chain is named joint z. Figures 1 and 2 give a general view of the designed structure and its kinematic parameters. The frame F with axis x, y and z is defined as world frame. To describe the position and orientation of the structure in space, a local frame F is attached to the end-effector. The orientation of this frame is defined with axis x aligned with the passive joint z axis and axis y aligned with the passive chain link attached to the drone 1. The orientation of the the local frame F is parametrized by three successive rotations about z (Yaw), y (Pitch) and x (Roll) axis. The revolute joints limit the internal motions of the structure to a plane orthogonal to the joint axes, the plane spanned by the vectors x and y . The full coordinate vector q = [x, y, z, ψ, θ, φ, q 21 , ϕ 1 , ϕ 2 ] T of the flying parallel robot is then composed of • x, y, z the origin of the local frame F expressed in world frame F, • ψ, θ, φ the representation in the global frame F of the local frame F orientation (Yaw/Pitch/Roll angles), • q 21 the relative angle between B 1 and B 2 ,

• ϕ 1 and ϕ 2 the roll angle of respectively drone 1 and drone 2.

Two local frames D 1 and D 2 are respectively attached to drone 1 and 2 center. Their respective orientations in the global frame F are given by (yaw/pitch/roll): (ψ, θ, ϕ 1 ) for drone 1, (ψ, θ, ϕ 2 ) for drone 2. The drones twists (velocity screws) t 1 and t 2 expressed in local drone frames are related to the coordinate velocities q through the Jacobian matrices J 1 and

J 2 t 1 t 2 = J 1 J 2 q = J q (1)
J is a (12 × 9) matrix. The matrix J is full rank, unless the two passive links of the structure are aligned. In this case, the Jacobian matrix J becomes rank deficient. Those singularities are the analogue of Type-2 singularities (Gosselin and [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF]) that appear on the biglide mech-anism when the two distal links (equivalents to bodies B 1 and B 2 ) are aligned. An uncontrolled motion of the effector appears in such singularities. They are not considered in this paper.

As previously said, the mechanism was designed to carry out planar vertical motions. Then, for the study of its dynamics, the robot is considered maintained close to a vertical plane x = 0, ψ = 0, θ = 0. In this vertical plane, the configuration of the passive chain is parametrized by q p = [y, z, φ, q 21 ] T . The velocity of joints x and y are related to the coordinates velocities qp through the Jacobian matrices J p1 and J p2

ṙ1 ṙ2 = J p1 J p2 qp = J p qp (2)
with r 1 and r 2 the respective positions of joint x and y in the vertical plane. J p is a (4 × 4) matrix. The other coordinates of the flying robot are regrouped in the attitude coordinate vector q a = [ψ, θ, ϕ 1 , ϕ 2 ] T . The translation along the axis x is treated separately (see section 3.3).

DYNAMIC MODEL

The dynamic model of a quadrotor presents a natural decoupling between the attitude and the translation dynamics [START_REF] Castillo | Stabilization of a mini rotorcraft with four rotors: Experimental implementation of linear and nonlinear control laws[END_REF]. This decoupling is preserved when the joint linking a rigid body and the drone is situated at the drone center of mass (CoM) as in [START_REF] Nguyen | Aerial tool operation system using quadrotors as Rotating Thrust Generators[END_REF]. The decoupling is an asset making backstepping control popular for drones. Unfortunately, placing a joint at the drone CoM for the studied structure would restrict the range for joint motion to avoid collision between the passive chain and the drone. Therefore, the joints are not situated at the drone CoMs and the dynamic decoupling is not obtained directly. However, expressing the dynamic equations for the passive chain in the vertical plane separately from the dynamic associated to the other coordinates will allow finding decoupling properties.

Drone input wrenches

Each drone provides an actuation wrench composed of three moments and one thrust force whose direction is normal to the drone propellers plane. Those inputs are related to the eight propeller angular velocities squared through a bijection as in [START_REF] Castillo | Stabilization of a mini rotorcraft with four rotors: Experimental implementation of linear and nonlinear control laws[END_REF]. The response time of the rotor dynamics is considered negligible compared to the drone and structure dynamics. The input wrenches in local frames D 1 and D 2 are given by

w 1 = [0 0 f 1z τ 1x τ 1y τ 1z ] T (3)
w 2 = [0 0 f 2z τ 2x τ 2y τ 2z ] T

Passive chain dynamic model

The passive chain motion restricted to the vertical plane (x = 0) is considered. The dynamic model of the passive chain will be computed via Euler-Lagrange approach. The Lagrangian of the passive chain is defined by L p (q p , qp ) = E p ( qp , q p ) -U p (q p ) (4) where E p is the kinetic energy of the passive chain and U p its potential energy. In the vertical plane (x = 0) two Fig. 3. Forces acting on one drone 1 planar forces f p1 and f p2 are exerted by the drones on the passive chain at respectively joints x and y. Euler-Lagrange equations applied to the passive chain gives

d dt ∂L ∂ qp T - ∂L ∂q p T = J T p1 f p1 + J T p2 f p2 (5)
with J p1 and J p2 the Jacobian matrices from (2). Using (4), equation ( 5) is expressed under the form

M p qp + c p = J T p1 f p1 + J T p2 f p2 = J T p f p1 f p2 (6) with
• M p the (4 × 4) definite positive generalized inertia matrix of the planar passive chain depending on q p , • c p the 4-dimensional vector of gravitational, Coriolis and centrifugal effects depending on q p and qp .

Away from singularities (see section 2), the Jacobian matrix J p in ( 6) is invertible and the inverse dynamic model of the passive chain is then

f p1 f p2 = J -T p (M p qp + c p ) (7)
The forces exerted by a drone on the passive chain can be linked to the drone inputs by expressing the translational dynamics of the drone i (i = 1, 2) in local frame (Fig. 3)

0 f iz = R i f ti + p i (8)
with

f ti = f pi + m di (r i -g) (9
)

p i = -m di d i φ2 i m di d i φi (10)
and

• g the gravitational acceleration vector,

• d i the distance between the joint i and the drone i CoM, • m di the mass of the drone i,

• R i = cos ϕ i sin ϕ i -sin ϕ i cos ϕ i the rotation matrix about -ϕ i
in the vertical plane.

Introducing ( 7) and ( 8) in (9) gives

R inv (Ψf -p) = M t qp + c t (11)
with

R inv = R -1 1 0 0 R -1 2 Ψ = 0 0 0 1 0 1 0 0 T p = p 1 p 2 c t = J -T p c p +     -m d1 l φ2 c 1 m d1 (-l φ2 s 1 + g) -m d2 lc 2 ( φ + q12 ) 2 m d2 (-ls 2 ( φ + q12 ) 2 + g)     M t = J -T p M p +    m d1 0 -m d1 ls 1 0 0 m d1 -m d1 lc 1 0 m d2 0 -m d2 ls 2 -m d2 ls 2 0 m d2 -m d2 lc 2 -m d2 lc 2   
and c 1 = cos φ, s 1 = sin φ, c 2 = cos(φ + q 21 ), s 2 = sin(φ + q 21 ). l is the passive link length (the same for the two arms). M t an invertible (4 × 4) matrix depending on q p .

The term p i in (8) generates a coupling between the passive chain dynamics and the drone dynamics. In the objective of designing a cascaded controller, this term will be considered as a perturbation of the passive chain dynamics. Then, the dynamic equations of the passive kinematic chain depends only on the input thrust forces (f 1z and f 2z ) and the drone roll coordinates (ϕ 1 and ϕ 2 ). This decoupling allows to consider the attitude coordinates as auxiliary inputs for those dynamics and lead to the control law established in section 4.

Attitude dynamic model

The dynamics of the attitude coordinates q a = [ψ, θ, ϕ 1 , ϕ 2 ] T is now studied. Let us now consider the Lagrangian of the complete robot L(q, q) = E( q, q) -U (q) (12) where E is the kinetic energy of the robot and U its potential energy. Euler-Lagrange equations applied to the robot gives

d dt ∂L ∂ q T - ∂L ∂q T = J T w 1 w 2 (13)
with J the Jacobian matrice from (1). Using ( 12), ( 13) is expressed under the form

Mq + c = J T w 1 w 2 (14) with
• M the (9 × 9) definite positive generalized inertia matrix of the structure, • c the 9-dimensional vector of gravitational, Coriolis and centrifugal effects.

The direct dynamic model of the complete structure is then obtained by inverting M in ( 14)

q = M -1 J T w 1 w 2 -c (15)
The acceleration of the attitude coordinates is selected in equation ( 15)

qa = M inv,a J T w 1 w 2 -c (16)
where M inv,a corresponds to the rows of the inverse of the matrix M restricted to the attitude coordinates. Introducing (3) in ( 16), the wrenches are expressed as function of drones input torques τ = [τ 1x , τ 1y , τ 1z , τ 2x , τ 2y , τ 2z ] T and drone input forces

f = [f 1z , f 2z ] T giving qa = M inv,a J T τ τ + J T f f -c = A τ τ + M inv,a J T f f -c (17)
Fig. 4. General scheme of the controller with J T τ and J T f the columns of the Jacobian J T restricted to the input torques τ and forces f respectively. The matrix A τ is a (4 × 6) matrix linking the attitude coordinates to the input torques. This matrix has the property to remain full rank in all configurations, allowing the feedback linearization of the attitude loop realized in section 4.3. Six torque inputs are available to generate four attitude angles. This is due to the fact that over-constraints are generated by the input moments of both drones about y and z axis. Those moments are acting along the same axis and are completely transmitted through the revolute joints of the passive structure. Then, several combinations of those input torques can generate the same acceleration on the attitude variables.

Planar motion and normal-to-plane dynamics

To maintain the vertical position of the structure, the dynamics of the system in the direction x is expressed. Considering ψ = 0 and |θ| << 1 the dynamic translation of the CoM of the flying structure is given by m t ẍG = sin(θ)(cos ϕ 1 f 1z + cos ϕ 2 f 2z ) θ(cos ϕ 1 f 1z + cos ϕ 2 f 2z ) (18) where m t is the total mass of the robot and x G is its CoM position along axis x.

CONTROL LAW

The control law designed for the flying robot is based on nonlinear feedback linearization already applied on classic quadrotors [START_REF] Voos | Nonlinear control of a quadrotor microuav using feedback-linearization[END_REF]. The control strategy will be the combination of two cascaded loops. A slow loop is dedicated to the trajectory tracking of the passive chain coordinates. The control inputs for this loop are the drone thrust forces and the desired attitude coordinates. Then a fast control loop ensures the convergence of the robot attitude coordinates toward the desired ones. A feedback linearization is applied in each loop and time scale separation of the two loops is ensured. An additional corrective term is also implemented to maintain the structure in the vertical plane x = 0. Controller scheme is shown in Fig. 4 4.1 Lateral control law From (18), θ can be used as an auxiliary input to maintain the CoM of the flying robot on the vertical plane x G = 0. The desired angle is chosen in order to obtain the closedloop equation ẍG

+ k px x G + k dx ẋG = 0 (19)
with gains k px > 0 and k dx > 0. Identifying (18) with the desired system closed-loop equation ( 19) gives the desired auxiliary input for θ

u θ = -k px x G -k dx ẋG m t (cos ϕ 1 f 1z + cos ϕ 2 f 2z ) (20)
The input thrust in ( 20) are those obtained in the computation of the passive chain control law ( 22). In order that the hypothesis of the flying robot maintained on a vertical plane remains valid, the gains k px and k dx are chosen sufficiently low. Then, the variations of the angular acceleration θ will not affect the dynamic model established for the passive chain. This control law is not established for tracking purpose but just to counter small perturbations along the direction x.

Passive chain control law

Consider a desired trajectory for the passive chain coordinates q d p . Consider an auxiliary control u p corresponding to those coordinate accelerations. A PD control law is defined for this auxiliary input u p = qd p -K pp e p -K dp ėp ( 21) with

• e p = q p -q d p the position tracking error, • K pp and K pd two definite positive matrices.

Ignoring the perturbation term p, from (11) the drone input thrust forces f = [f 1z , f 2z ] T required to produce the desired auxiliary inputs u p are computed by

f = -sin ϕ 1 cos ϕ 1 0 0 0 0 -sin ϕ 2 cos ϕ 2 f u1 f u2 (22) f u1 f u2 = M t u p + c t
The closed-loop equation of the passive chain is obtained by applying the computed input forces ( 22) to its dynamic model ( 11).

(M t u p + c t ) -δ -R inv p = M t qp + c t (23) where δ = R inv    f u1y cos ϕ 1 + f u1z sin ϕ 1 0 f u2y cos ϕ 2 + f u2z sin ϕ 2 0    (24)
is a perturbation term due to the drones underactuation.

The term δ can be canceled by choosing appropriate drone roll angles. Two auxiliary inputs u ϕi corresponding the drone roll angles are defined u ϕi = atan2(f uiy , -f uiz ) (25) From ( 24) and ( 25), we have δ → 0 when ϕ i → u ϕi for i = 1, 2. M t being an invertible matrix, the closed loopequation obtained by introducing ( 21) in ( 23) is ëp + K pp e p + K dp ėp + M -1 t (δ + δ p ) = 0 (26) The perturbation δ p = R inv p associated to the drones angular velocities and accelerations will be treated in section 4.4. Considering this perturbation negligible, the attitude loop will be controlled to make each drone roll angle ϕ i converge toward the desired input value u ϕi . Then, the closed-loop of the passive chain behaves as a secondary ordinary differentiation equation ensuring the convergence of e p and ėp towards zero.

Attitude control law

Similarly to the passive chain control loop, the attitude control loop uses dynamic inversion to converge toward the desired configuration. The desired configuration vector q d a is contains the terms • ψ d = 0 to maintain the robot on its vertical plane, • u θ from (20) to generate the correction along the direction x, • u ϕ1 and u ϕ2 from ( 25) to obtain the convergence on the passive chain loop.

Through dynamic inversion, a PD control law is designed on the auxiliary input qa u a = qd a -K pa e a -K da ėa ( 27) with e a = q a -q d a the position tracking error and K pa and K da two definite positive matrices. Taking into account the over-actuation discussed in section 3.3, the input torques τ are then computed with ( 17)

τ = A + τ (u a + M inv,a c -M inv,a J T f f ) (28
) with f the input thrust vector issued from ( 22). Using the pseudo-inverse A + τ , the solution minimizing the input torques norm is chosen among the possible solutions. As the inputs, f , τ are applied to the robot, the attitude closed-loop converges as a second order ode.

Perturbation handling

A perturbation term δ p has been identified in the closedloop equation of the passive structure control (26). From equation ( 8), this perturbation can be partly rejected by replacing f iz the thrust inputs computed in equation ( 22)

f = -sin ϕ 1 cos ϕ 1 0 0 0 0 -sin ϕ 2 cos ϕ 2 f u1 f u2 + m d1 d 1 φ1 m d2 d 2 φ2 (29)
In this case, the perturbation term becomes

δ p = R inv -m di d 1 φ2 1 0 -m di d 2 φ2 2 0 T (30)
The perturbation associated to the drone acceleration have been rejected. This modification of the thrust input forces do not impact the attitude loop. The remaining perturbation δ p that could not be rejected converges toward zero as the drone angular velocity converges toward zero. In static equilibrium, simulations showed that this term does not impact the stability. However, it affects the tracking performance for fast motions.

SIMULATIONS

Simulation design

A simulator was designed using MATLAB-SIMULINK and ADAMS co-simulation. The flying structure dynamic simulation was ensured by ADAMS while the control loop was implemented on SIMULINK. A discrete communication time (1 kHz) is set between the two software to represent the measurement sampling. The simulator considers that the coordinate vector q is fully measured. Such measures could be obtained on a real prototype using a vision-based pose estimation as in [START_REF] Achtelik | Visual tracking and control of a quadcopter using a stereo camera system and inertial sensors[END_REF] and encoders on passive joints.

To test the robustness of the controller against noisy pose estimation, a uniformly distributed noise was added to Fig. 6. Stabilization of the passive kinematic chain from an unstable initial position each variables of the pose estimation and encoder resolution was also taken into account. The interval of the uniform distribution in the simulation results presented is -1e -3 ; 1e -3 rad for the orientation coordinate estimations and -1e -3 ; 1e -3 m for the position estimations. The robustness of the controller was observed with noise range up to -1e -2 ; 1e -2 . As the feedback linearizations lead to equivalent second order ordinary differential equations, the PD gains were set to ensure a unitary relative dampling. The tuning parameter is then the natural pulsation ω n leading to proportional gain ω 2 n and derivative gain 2ω n . The natural pulsation selected for each controlled coordinate are for drone roll: 20 rad/s, for attitude (ψ, θ): 7 rad/s, for the passive chain coordinates: 3 rad/s, for lateral control: 2 rad/s. The natural pulsation associated the variables used to maintain the flying structure in a vertical plane are voluntary chosen low to limit perturbations over the dynamic model of the passive chain. The gains are also selected to ensure a time scale separation between the loops. The time derivatives of the desired configuration q d a for the attitude loop are numerically computed and filtered with lowpass Butterworth filters.

Stabilization

A first simulation shows the stabilization of the flying structure from an initial position. The structure is placed in an unstable initial position as in Fig. 1. The two drones are in a horizontal position. As the drones cannot compensate the lateral forces generated by the passive chain, the position is unstable. With a constant input on the passive chain desired position, the controller will drive the drones angles to a constant position assuring the stabilization of the structure. Fig. 5 shows the desired drone roll angles and the convergence toward those desired angles. Fig. 6 shows the error obtained on the passive kinematic chain parameters q p . As the initial position is not stable, those parameters diverge from their initial position but are driven back by the controller when the convergence is obtained on the attitude loop.

Trajectory tracking

From a stabilized position, a square trajectory was designed in the plane (x = 0). At each corner of the square, different angles φ and q 21 are desired on the passive chain.

The trajectory and its first to forth time derivatives are designed to be continuous on all controlled variables. The desired duration for the robot to travel each edge of the square is defined by the time t e . Fig. 7 shows the desired trajectory and obtained trajectory with t e = 3 s. In simulation, the controller shows good performances in maintaining the passive chain on the desired trajectory and remains robust against noisy measures. Figure 8 compares the evolution of the position tracking error and drone angular velocity along the trajectory. The perturbation term associated to the drone angular velocity clearly impact the trajectory tracking precision, with phase shift as it acts on the passive structure acceleration. Nevertheless, it remains sufficiently low to not affect the controller stability. The inputs wrenches (3) are linearly related to the force provided by each rotor. In this simulation, the force provided by each propeller is in range [4.8,8.5] N.

CONCLUSION AND FUTURE WORK

In this paper, a new type of flying parallel robot was studied. The robot is composed of two quadrotors linked by a rigid articulated passive chain, that can be assimilated to a flying parallel robot. The study of the dynamic model showed decoupling properties. Those properties have been exploited to design a controller dedicated to this new flying robot. Simulations showed the performance and the robustness against noisy pose estimation of the controller designed. The promising results obtained provide a good base for the future development of flying parallel robots. Future work on the architecture presented in this paper includes the design of optimal trajectories with respect to the drones payload; the consideration of the impact of the drone centrifugal force to avoid unstability and to increase accuracy; application of reconfiguration techniques known for parallel robots on the flying robot; design of a position/force controller to interact with the environment. Such technique will allow positioning the tool tip over or under the drones, in function of the task objective. A prototype of the robot is also under construction to test the controller developed in a real environment.

Fig. 1 .

 1 Fig. 1. Snapshot of the flying parallel robot simulator realized with the software ADAMS

Fig. 2 .

 2 Fig.2. Parametrisation of the flying parallel robot of quadrotors combined with the coupling induced by the rigid chain is a new challenge for the modelling and control of the flying robot. To reduce the complexity of the study, the robot is designed to perform motions restricted to a vertical plane. Precautions taken in the control law to maintain the flying structure in this plane are discussed in section 4. The flying robot is then derived from a planar parallel mechanism, the biglide mechanism. This flying robot obtained is composed of two drones (named drone 1 and drone 2) and a two links passive kinematic chain (B 1 and B 2 ). Revolute joints are selected to constrain the motion of the passive kinematic chain to a planar motion. The passive joints connecting the passive chain to the drone 1 and 2 are respectively named joint x and y. The passive joint corresponding to an end-effector in the middle of the passive kinematic chain is named joint z. Figures1 and 2give a general view of the designed structure and its kinematic parameters. The frame F with axis x, y and z is defined as world frame. To describe the position and orientation of the structure in space, a local frame F is attached to the end-effector. The orientation of this frame is defined with axis x aligned with the passive joint z axis and axis y aligned with the passive chain link attached to the drone 1. The orientation of the the local frame F is parametrized by three successive rotations about z (Yaw), y (Pitch) and x (Roll) axis. The revolute joints limit the internal motions of the structure to a plane orthogonal to the joint axes, the plane spanned by the vectors x and y . The full coordinate vector q = [x, y, z, ψ, θ, φ, q 21 , ϕ 1 , ϕ 2 ] T of the flying parallel robot is then composed of

  Fig. 5. Convergence on drone roll angles

Fig. 8 .

 8 Fig. 7. Desired and obtained motion for t e = 3 s


S Briot

A Goldsztejn
email: alexandre.goldsztejn@ls2n.fr

Topology Optimization of a Reactionless Four-bar Linkage

Keywords: Four-bar linkage, shaking force and shaking moment balancing, optimal design, topology optimization

Most of existing works on the optimal design of balanced four-bar linkages deal essentially with the minimization of their inertia or input torques under balancing constraints. These approaches are not suitable to include constraints based on the elastic behavior of the mechanism. In order to solve this issue, we propose in this paper to perform the topology optimization of a reactionless four-bar linkage. Conditions for balancing the mechanism are first recalled and a topology optimization algorithm is run in order to maximize the first natural frequency while ensuring the balancing and constraining the mechanism compliance. We show that in order to obtain an interesting design solution, it is necessary to modify the balancing constraints in order to penalize them. Interesting design solutions are obtained in a rather short computational time.

Introduction

Transmitting no reactions to the ground is very appealing in many applications such that space robotics or high-speed manipulation [START_REF] Jiang | Dynamic optimization of reactionless four-bar linkages[END_REF]. However, complete shaking force and shaking moment balancing is usually obtained by using both counterweights and counter-rotations, thus leading to an increase in the design complexity and to noise and backlash issues due to the use of gears [START_REF] Jiang | Dynamic optimization of reactionless four-bar linkages[END_REF].

In order to avoid the drawbacks in using counter-rotations, a reactionless four-bar linkage without them was proposed in [START_REF] Ricard | On the development of reactionless parallel manipulators[END_REF]. For obtaining this property, conditions on both the geometry and the mass distribution of the linkage must be respected. It was later shown in [START_REF] Gosselin | Synthesis and design of reactionless three-degreeof-freedom parallel mechanisms[END_REF] that this reactionless linkage can be used as an elementary block in order to design reactionless robots In [START_REF] Jiang | Dynamic optimization of reactionless four-bar linkages[END_REF], the optimal design of the reactionless four-bar linkage was carried out for minimizing its input torques under balancing conditions constraints only. In this work, the shape of the links and counterweights is fixed and the authors focus on the optimal positioning of the counterweights. Other works deal with the optimal design of balanced four-bar linkages for allowing the full [START_REF] Farmani | Multiobjective optimization for force and moment balance of a four-bar linkage using evolutionary algorithms[END_REF] or partial [START_REF] Chaudhary | Shape optimization of dynamically balanced planar four-bar mechanism[END_REF][START_REF] Erkaya | Investigation of balancing problem for a planar mechanism using genetic algorithm[END_REF] dynamic balancing while optimizing energy consumption or input torque amplitude. More complete list of reference can be found in [START_REF] Arakelian | Balancing of Linkages and Robot Manipulators -Advanced Methods with Illustrative Examples[END_REF].

The main issue with the aforementioned methods comes from the fact that the shape of the links is already fixed (the design variables are their dimension) and deformation or vibration constraints may lead to an unfeasible design in practice (bulky mechanism to resist to the external efforts while ensuring the balancing conditions). For avoiding this issue, it is necessary to optimize the shape of the links. This was done in [START_REF] Chaudhary | Shape optimization of dynamically balanced planar four-bar mechanism[END_REF] for the partial balancing of the four-bar linkage, however the approach proposed allows for finding only the external shape of the links. No voids can be included, which is not optimal w.r.t. the minimization of the link deformation while ensuring the decrease of the link mass.

ψ 1 ψ 2 x 0 ψ 3 y 0 S 1 S 3 S 2 x 1 y 1 x 2 y 2 x 3 y 3 O 1 O 0 O 2 O 3 O' 2 θ 1 (a) General scheme x 0 y 0 O 0 x i y i O i M ij m ij Body B i y ij x ij (b) Location of the finite element i j in the body B i
Performing a Topology Optimization (TO) of the linkage can solve this issue [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]. This technique aims at optimizing the material distribution in a link in order to satisfy performance criteria including deformations or vibration constraints.

Our contribution in the present paper is to perform the TO of the reactionless four-bar linkage for ensuring its full balancing conditions while constraining elastic performance criteria. To the best of our knowledge, this is the first time that the TO of a fully force and moment balanced mechanism is performed. Furthermore, we show that balancing constraints must be partially penalized in order to obtain meaningful designs.

Problem formulation

The general scheme of the four-bar linkage is given in Fig. 1 • r i is the distance between O i and S i (i = 1, 2, 3) • for body i, S i is the center of mass, m i is its mass, zz i is the moment of inertia expressed at O i , mx i (my i , resp.) is the static moment around x i (y i , resp.), i.e

mx i my i = B i ---→ O i M i j dm = m i r i cos ψ i sin ψ i . (1)

Shaking force and shaking moment balancing conditions

The conditions given in [START_REF] Ricard | On the development of reactionless parallel manipulators[END_REF] for achieving the full shaking force and shaking moment balancing of the four-bar linkage without counter-rotations involve both contraints on the mechanism geometry and mass distribution. Three different set of links lengths are possible: Set 1: 1 = 3 and 2 = 4 ; Set 2: 1 = 2 and 3 = 4 ; Set 3: 1 = 4 and 2 = 3 .

In what follows, we focus only on mechanisms designed with the first set of geometric constraints, which is the set most often studied in the papers (see for instance [START_REF] Gosselin | Synthesis and design of reactionless three-degreeof-freedom parallel mechanisms[END_REF][START_REF] Jiang | Dynamic optimization of reactionless four-bar linkages[END_REF]) and corresponds to the anti-parallelogram linkage. In this case, the conditions on the mass distribution given in [START_REF] Ricard | On the development of reactionless parallel manipulators[END_REF] for the full balancing can be rewritten as: my 1 = 0 and my 2 = 0 and my 3 = 0 (2)

mx 1 / 1 + m 2 -mx 2 / 2 = 0 and mx 3 / 3 + mx 2 / 2 = 0 (3) zz 1 -mx 1 1 + zz 2 -mx 2 2 = 0 and zz 3 -mx 3 3 + zz 2 -mx 2 2 = 0 (4)

Modeling of the linkage elastic behavior

Topology optimization uses the same physical model as in the finite element method for modeling bodies and linkages: each body is meshed using finite elements. The presence or absence of an element i j (i.e. the element j of the body B i ) is parameterized with a material density variable ρ i j which can take values between 0 (which represents the absence of material) and 1 (which represents the presence of material).

Based on these density variables, we use an interpolation scheme 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 or ρ i j = 0 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 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] (5)
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].

The computation of the body and linkage mass matrices (necessary for the computation of the elastodynamic performance) is also based on the use of elementary mass matrices M i j equal to

M i j = ρ i j M 0 i j (6)
where M 0 i j is the mass matrix of the element i j computed for a density ρ i j = 1. Once the linkage stiffness and mass matrices are obtained, the elastic performance of the mechanism can be defined, such that the deformations or natural frequencies [START_REF] Shabana | Dynamics of Multibody Systems[END_REF].

It should be mentioned that, in order to decrease the time for computing the elastic performance of the linkage, a Craig-Bampton model reduction technique [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] is applied on each body independently, as done in [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF].

Optimization problem

The general formulation of a mono-objective TO problem is given by:

min ρ f (ρ) subject to ρ ∈ [0, 1] n , g(ρ) ≤ 0, h(ρ) = 0, (7)
where

• ρ is the decision variable vector containing all element densities ρ i j ,

• f , g and h are functions of ρ characterizing performance indices or structure constraints.

Several optimization algorithms can be used, among which we can cite the Optimality Criteria method [START_REF] Yin | Optimality criteria method for topology optimization under multiple constraints[END_REF], the Method of Moving Asymptotes [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF] or the Linearization Method [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF] (LM) only recently used in the context of TO in [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF]. We used the latter in this work. All these techniques have the same specificity: they require the symbolic computation of the functions and their gradients.

Here, we decided to maximize the first natural frequency for the linkage under the following constraints: • the balancing equalities (2) to (4) must be respected. As shown in [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF], the inertial parameters used in these equalities are all linear with respect to the decision variables in ρ. • the compliance (i.e. twice the potential elastic energy or also the dot product of the nodal wrenches by the nodal displacements) must be lower than a given threshold under a given loading (as often done in TO, see for instance [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]).

For computing the compliance and natural frequency, we consider that the body B 1 is actuated in O 1 and that the computation of these performances is made for θ 1 = π/2 (Fig. 1(a)).

3 Topology optimization

Initial domain

The initial design domain for the proximal and distal links is represented in Figs. Fig. 2 Initial design domain and first temptative optimal design for all bodies (black elements correspond to ρ i j = 1, white elements to ρ i j = 0, and gray elements to 0 < ρ i j < 1)

Links are considered to be made of steel with Young's modulus E 0 = 210 GPa, Poisson's ratio ν = 0.

Optimization results

We run the TO algorithm with a fixed threshold for the compliance value of 1.6 • 10 -3 Nm for the assembled linkage under the following loading:

• at O 2 : a force of 15 N along x 0 , a force of 15 N along y 0 , a moment of 2 Nm around z 0 on body B 1 , a moment of 1 Nm around z 0 on body B 2 , • at O 2 : a force of 15 N along x 0 , a force of 15 N along y 0 , a moment of 2 Nm around z 0 on body B 3 , a moment of 1 Nm around z 0 on body B 2 , • at O 3 : a moment of 2 Nm around z 0 on body B 3 . As usually done in TO, in order to obtain a smoother layout without checkerboards problem, we modified the density variables assigned to the elements with the information of its neighborhoods as was proposed in [START_REF] Bourdin | Filters in topology optimization[END_REF].

All functions have been encoded with Matlab in the Windows 7 environment. The LM optimization algorithm is available under request to the second author.

The optimization process is run and we considered that the algorithm converged when the maximal change between two sequential iterations for any component of the density vector ρ is lower than 0.01. First results of optimization are shown in Fig. 3 Optimal design of bodies B 1 , B 2 and B 3 in two cases; case 1: the initial design domain of body B 2 is the one depicted in Fig. 2; case 2: the initial design domain of body B 2 is is changed: the body's height is now of 30 mm instead of 20 mm as in case 1 (black elements correspond to ρ i j = 1, white elements to ρ i j = 0, and gray elements to 0 < ρ i j < 1)

Figs. 2(b), 2(d) and 2(f). The algorithm stopped after 524 iterations, with a maximal constraint violation of 5 • 10 -3 %1 . In totality, the optimization procedure took 28 minutes with a Pentium 4 2.70 GHz, 16 GB of RAM. However, the algorithm had difficulty to converge (large oscillations in the value of the objective function, not displayed for reasons of page limitations) and finally attained a first natural frequency of 67 Hz (which is quite low, as will be shown below).

Obtained results showed that on the left-hand area of the points O 1 and O 3 , the material density for bodies B 1 and B 3 is between 0 and 1 (gray elements), thus leading to bodies which are not easy to design by engineers [START_REF] Sigmund | A 99 line topology optimization code written in MATLAB[END_REF]. We increased the number of allowed optimization iterations and obtained no improvement: These portions of materials are mostly here to fulfill the balancing constraints and have less impact on the compliance constraint or the frequency of the full linkage.

We propose here a way to improve this solution which is based on a partial penalization of the balancing constraint equations. As said in Sect. 2.3, the equality constraints (2), ( 3) and ( 4) are linear, i.e they can be written under the form:

Aρ = A l ρ l + A r ρ r = 0 (8)
in which ρ l contains the decision variables associated with the elements on the lefthand side of the gray lines in Figs. 2(a) and 2(e) for bodies B 1 and B 3 while ρ r contains all other variables, including those of the body B 2 . Thus, the vector ρ l contains the variables associated with the portions of materials which are almost here to fulfill the balancing constraints, which have few impacts on the linkage elastic performance, and which takes intermediary values for density.

In order to force the values of the variables in ρ l to be only 0 or 1, we modify the balancing equality (8) by raising the variables ρ l at the power of q as follows:

A l ρ q l + A r ρ r = 0 (9)
In our code, we put q = 2 or 3. To the best of our knowledge, this is the first time that penalization method is applied in order to achieve balancing constraints. Based on this new formulation, condition (9) does well represent the balancing equality [START_REF] Farmani | Multiobjective optimization for force and moment balance of a four-bar linkage using evolutionary algorithms[END_REF] if all elements in ρ l are equal to 0 or 1. Thanks to this penalization of the variables ρ l , a small removal of material has a considerable impact on the balancing of the system, thus forcing the algorithm to impose 0 or 1 values to the variables ρ l in order to counterbalance the effect of the variables ρ r . The optimization results by taking into account the constraints (9) instead of ( 8

Discussion

This work was made in order to show that TO can be used in order to obtain solutions in a rather short time of a complicated design problem which was to balance a fourbar linkage while ensuring deformation, compliance or frequency constraints.

However, in this paper, some issues have not been solved, which should be investigated later. First, the optimization is performed for compliance and frequency computed at θ 1 = π/2. Thus, our optimization will be local by default. A more global optimization procedure ensuring the performance for any linkage configuration could be used (see [START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF] for instance) and some other optimization problems could be envisaged: objective and constraints could be changed.

It should be finally mentioned that our simulations have shown that the convergence of the algorithm is considerably sensitive by the threshold on the inequality constraints. In case all constraints are not achievable (i.e. there is no solution to the problem), the algorithm may become unstable.

Conclusion

Most of existing works on the optimal design of balanced four-bar linkages deal essentially with the minimization of their inertia or input torques under balancing constraints. These approaches are not suitable to include constraints based on the elastic behavior of the mechanism. In order to solve this issue, we performed in this paper the topology optimization of a reactionless four-bar linkage.

In our paper, a topology optimization algorithm was run in order to maximize the first natural frequency while ensuring the balancing and constraining the mechanism compliance. We showed that in order to obtain an interesting design solution, it was necessary to modify the balancing constraints in order to penalize them. Interesting design solutions were obtained in a rather short computational time.

Future works include solving the problem in 3-D and also carrying out multimaterial topology optimization in order to ensure the balancing conditions by partially using materials with higher density leading thus to smaller mechanism footprint. The design of a prototype is also envisaged.

Fig. 1 A

 1 Fig.1A general four-bar linkage.

2 3 y 3 O 3 O' 2 (

 23332 2(a), 2(c) and 2(e). Four-bar geometric parameters are taken at 1 = 3 = 60 mm and 2 = 4 = 200 mm. Each link has got two holes (voids) in which the joints will be inserted. Bodies B 1 and B 3 external shapes are rectangles of dimension (150 × 40) mm and thickness of 10 mm, while body B 2 external shape is a rectangle of size (220 × 20) mm and thickness of 1 mm. For the meshing of the links, QUA4 finite elements (i.e. four-nodes rectangular planar elements) of size 1×1 mm. Elements associated with ρ l Elements associated with ρ r (a) Initial design domain for body B 1 First version of optimal design for body B Elements associated with ρ l Elements associated with ρ r (e) Initial design domain for body B 3 x f) First version of optimal design for body B 3

  3 and density of 7800 kg/m 3 . As a result, 6000 elements are used for meshing the bodies B 1 and B 3 while the body B 2 is made of 4400 elements (Figs. 2(a), 2(c) and 2(e)).

2 (

 2 f) Optimal design of body B 3 : case 2

  ) are shown in Figs. 3(a), 3(c) and 3(e) (results obtained in 20 minutes, objective: first natural frequency of 646 Hz, constraint violation of 1.4 • 10 -4 %). Results were obtained without any instability of the optimization algorithm and it can be observed that gray elements have been removed from the design solution. Figures. 3(b), 3(d) and 3(f) show the same optimization problem but the difference comes from a change in the initial domain for body B 2 which was increased (body's height is now of 30 mm instead of 20 mm). Final objective was of 763 Hz and was attained in around 10 hours. It can be observed that a slight change in the design domain may lead to a totally different design solution.

1 is the distance between O 1 and O 2 , 2 is the distance between O 2 and O 2 and 3 is the distance between O 3 and O 2 , 4 is the distance between O 1 and O 3 ,

  (a). In what follows, • body B 1 is the body between joints at O 1 and O 2 , body B 2 is the body between joints at O 2 and O 2 and body B 3 is the body between joints at O 3 and O 2 , •

Constraints are normalized using their values computed for the initial design, except for Eq. (2) whose initial values are null.

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