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

Vigen Arakelian

ON THE DYNAMIC PROPERTIES OF FLEXIBLE PARALLEL MANIPULATORS IN THE PRESENCE OF PAYLOAD AND TYPE 2 SINGULARITIES

, the condition of passing through a Type 2 singularity for parallel manipulators with flexible joints has been studied.

In the present paper, we expand information about the dynamic properties of parallel manipulators in the presence of Type 2 singularity by including in the studied problem the link flexibility and the payload. The suggested technique is illustrated by a 5R parallel manipulator with flexible elements (actuated joints and moving links) and a payload. The obtained results are validated by numerical simulations carried out using the software ADAMS.

1

INTRODUCTION

There are many papers dealing with the singularity analysis of parallel manipulators. We present some of them in the remaining.

The analysis of singular configurations has been first discussed from a kinematic point of view [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF]- [START_REF] Briot | Singularity Analysis of PAMINSA Manipulator[END_REF]. However, it is also known that, when parallel manipulators have Type 2 singularities [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF], they lose their stiffness and their quality of motion transmission, and as a result, their payload capability. Therefore, the singularity zones in the workspace of manipulators may be analyzed not only in terms of kinematic criterions, from the theoretically perfect model of manipulators, but also in terms of kinetostatic approaches [START_REF] Gosselin | The Optimum Design of Robotic Manipulators using Dexterity Indices[END_REF]- [START_REF] Hubert | Singularity Analysis through Static Analysis[END_REF].

The further study of singularity in parallel manipulators has revealed an interesting problem that concerns the path planning of parallel manipulators under the presence of singular positions, i.e. the motion feasibility in the neighbourhood of singularities. In this case the dynamic conditions can be considered in the path planning process. One of the most evident solutions for the stable motion generation in the neighbourhood of singularities is to use redundant sensors and actuators [START_REF] Alvan | On the Control of the Spatial Parallel Manipulators with Several Degrees of Freedom[END_REF]- [START_REF] Hesselbach | Aspects on the Design of High Precision Parallel Robots[END_REF]. However, it is an expensive solution to the problem because of the additional actuators and the complicated control of the manipulator caused by actuation redundancy. Another approach concerns with motion planning to pass through singularity [START_REF] Bhattacharya | Comparison of an Exact and Approximate Method of Singularity Avoidance in Platform Type Parallel Manipulators[END_REF]- [START_REF] Perng | Inverse Kinematic Solutions a Fully Parallel Robot with Singularity Robustness[END_REF],

i.e. a parallel manipulator may track a path through singular poses if its velocity and acceleration are properly constrained. This is a promising way for the solution of this problem. However, only a few research papers on this approach have addressed the path planning for obtaining a good tracking performance. But they have not adequately addressed the physical interpretation of dynamic aspects.

In [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF], optimal force generation in parallel manipulators for passing through the singular positions has been studied. It was shown that any parallel manipulator can pass through the singular positions without perturbation of motion if the wrench applied on the end-effector by the legs and external efforts of the manipulator are orthogonal to the twist along the direction of the uncontrollable motion. This approach has been generalised in the case of rigid-link flexible-joints parallel manipulators [START_REF] Briot | On the Dynamic Properties of Rigid-Link Flexible-Joint Parallel Manipulators in the Presence of Type 2 Singularities[END_REF].

This study is the continuation of works [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF], [START_REF] Briot | On the Dynamic Properties of Rigid-Link Flexible-Joint Parallel Manipulators in the Presence of Type 2 Singularities[END_REF]. The purpose of this paper is to study the dynamic properties of parallel manipulators with flexible links and joints in the presence of preponderant payload, which can be the main source of elastic deformations.

The paper is organized as follows. The next section presents theoretical aspects of the examined problem, using the Lagrangian formulation. The condition of force distribution is defined, that allows the passing of any parallel manipulator through the Type 2 singular positions. In section 3, the suggested solution is illustrated via 5R planar parallel manipulator. In section 4, the conclusions are given.

OPTIMAL DYNAMIC CONDITIONS FOR PASSING THROUGH TYPE SINGULARITY

Let us consider a parallel manipulator of m links, n degrees of freedom and driven by n actuators, for which the effects of the payload and of the geared motors inertia are preponderant from the others.

The general Lagrangian dynamic formulation for a nonrigid manipulator can be expressed as [START_REF] Bouzgarrou | New Approach for Dynamic Modelling of Flexible Manipulators[END_REF]:

a a L L q q τ                dt d , ( 1a
) e e L L q q 0                dt d , (1b)
where, -L is the Lagrangian of the manipulator; L = T-V, where T is the kinetic energy and V the potential energy due to gravitational forces, friction and elasticity;

- T a n a a a q q q ] ,..., , [ 2 1  q and T a n a a a q q q ] ,..., , [ 2 1      q
represent the vectors of position and velocity of the actuators, respectively; -T e n e e e q q q ] ,..., , [

 q and T e n e e e q q q ] ,..., , [

2 1      q
represent the vectors of position and velocity of the elastic coordinates (deformations of links and joints); -  is the vector of the actuators efforts.

In the case where the load on the end-effector and the inertia of actuators are preponderant, the expressions of the kinetic energy T and potential energy V are given by:

a a T a T T q I q v M v     2 , ( 2a
) e g V V V   , z g m V g  and ) ( e e f V q  (2b)
where V g corresponds to the potential energy due to gravity, and V e to the energy due to elastic deformations. M is the mass matrix of the payload (comprising both mass and rotational inertia terms), I a the diagonal matrix containing the axial moments of inertia of the geared motors, m the mass of the load, g the gravity, f(q e ) being a function depending of the deformations q e . Moreover,

T z y x ] , , , , , [     x and T z y x ] , , , , , [           v
are vectors containing the trajectory parameters and their derivatives, respectively; x, y, z represent the position of the controlled point in the global frame and  and  the rotation of the platform about three axes a   a  and a  . Vector x depends on both rigid coordinates q a and elastic coordinates q e . Analyzing expression (2), the potential and kinetic energies do not explicitly depend both of the actuated variables q a and elastic coordinates q e , but also from the positions x and velocities v of the payload. Therefore it is preferable to rewrite Eq. ( 1) using the Lagrange multipliers [START_REF] Bouzgarrou | New Approach for Dynamic Modelling of Flexible Manipulators[END_REF], as follows:

λ B W τ b T   , a a a a L L q I q q W b                   dt d (3a) λ C W 0 c T   , e e L L q q W c                dt d ( 3b
)
where  is the Lagrange multipliers vector, which is related to the wrench W p applied on the platform by:

p W A   T ,                     x v W p L L dt d (4)
and,  - A, B and C are three matrices relating the vectors v, e q  and a q  according to

e a q C q B Av    
. They can be found by differentiating the closure equations f i (x, q a , q e ) = 0 (taking into account the rigid as well as the elastic coordinates [START_REF] Bouzgarrou | New Approach for Dynamic Modelling of Flexible Manipulators[END_REF]) with respect to time. In the hypothesis of small elastic displacements ( 0 q  e ), matrices A and B may be found assuming that the robot is composed of rigid links only.

-W p is the wrench applied on the platform by the legs and external forces expressed along axes a   a  and a   [START_REF] Khalil | A Novel So-lution for the Dynamic Modeling of Gough-Stewart Manipulators[END_REF]. Expressing W p in the base frame, one can obtain:

p R q W J q I τ 0 T a a a     , ( 5a
) p R q c W J W 0 0 T e   , (5b)
where

  B A J R q 1 0   a
is the square Jacobian matrix between twist t of the platform (expressed in the base frame) and the vector a q  of actuators velocities,

  C A J R q 1 0   e
is the nonsquare Jacobian matrix between twist t of the platform (expressed in the base frame) and the vector e q  of deformations velocities, D A A R  0 is the expression of matrix A in the base frame, where D is a transformation matrix, of which expression is given in [START_REF] Merlet | Parallel Robots[END_REF].

For any prescribed trajectory x(t), the values of vector q a can be found using the inverse kinematics and dynamics. Thus, taking into account that the manipulator is not in a Type 1 singularity [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF], the terms W c and p R W 0 can be computed (using, or not, some recursive algorithm [START_REF] Boyer | Simulation of Flexible Manipulators Using Newton-Euler Inverse Dynamic Model[END_REF]). However, for a trajectory passing through a Type 2 singularity, the determinant of matrix A R 0 vanishes. Numerically, the values of the efforts applied by the actuators become infinite. In practice, the manipulator either is locked in such a position of the end-effector or it can not follow the prescribed trajectory.

As previously mentioned, in a Type 2 singularity, the determinant of matrix A R 0 vanishes. In other words, at least two of its columns are linearly dependant [START_REF] Merlet | Parallel Robots[END_REF]. So, one may obtain such a relationship:

   6 1 j j j 0 A  , ( 6
)
where A j represents the j-th column of matrix A R 0 and  j are coefficients, which in general can be functions of q a and q e . It should be noted that the vector t s = [   2 , …,  6 ] T represents the direction of the uncontrollable motion of the platform in a Type 2 singularity.

By substituting ( 6) into (4), we obtain

j T j W  λ A , j = 1, …, 6 (7)
where W j is the j-th row of vector p R W 0 . Then, from ( 6) and ( 7) the following conditions are derived:

          6 1 6 1 0 j j j j T j j W   λ A . ( 8
)
The right term of Eq. ( 8) corresponds to the scalar product of vectors t s and Otherwise, the dynamic model is not consistent. Obviously, in the presence of a Type 2 singularity, the displacement of the end-effector of the manipulator has to be planned to satisfy [START_REF] Zlatanov | Con-straint Singularities of Parallel Mechanisms[END_REF]. Therefore, our task will be to achieve a trajectory which will allow the manipulator passing trough the Type 2 singularities, i.e. which will allow the manipulator respecting condition [START_REF] Zlatanov | Con-straint Singularities of Parallel Mechanisms[END_REF].

In the next section, an example illustrates the obtained results discussed above. This example presents a planar 5R parallel manipulator.

ILLUSTRATIVE EXAMPLE

In the planar 5R parallel manipulator, as shown in Fig. 1, the output axis is connected to the base by two legs, each of which consists of three revolute joints and two links. In of the two legs, revolute joint connected to the base is actuated. Thus, such a manipulator is able to position its output axis in a plane.

As shown in Fig. 1, the input joints are denoted as A and E. The orientation of elements 1 and 4 are denoted e q 1 and e q 2 , respectively. The common joint of the two legs is denoted as C, which is also the output axis with controlled parameters x and y. A fixed global reference system xOy is located at the middle of segment AE with the y-axis normal to AE and the xaxis directed along AE. The lengths of the links AB, BC, CD, DE are respectively denoted as L 1 , L 2 , L 3 and L 4 . Actuators 1 and 2 are connected to links 1 and 4, respectively, via Harmonic Drive® systems which are represented by a model similar to that given in [START_REF] Spong | An Integral Manifold Approach to the Feedback Control of Flexible Joint Robots[END_REF]. The position of actuator i is denoted as a i q . It is assumed that the actuator i is capable to deliver a couple  i to the motor shaft, which is elastically coupled to the link j of the robot (i = 1 or 2, j = 1 or 4). The flexibility of the drive system is modeled by a torsion spring with stiffness k 1 . The gear ratio is denoted n. I a is the axial moment of inertia of the motor i plus the Harmonic drive system.

In order to obtain relatively simple symbolic model for demonstrating the expected results, we propose to modelize the deformations of the robot by adding virtual joints on the links. We assume here that the link deformations of elements 2 and 3 are preponderant. Therefore, we add on them linear springs of stiffness k 2 , that are directed along directions BC and CD, respectively. The displacements of the spring linked to element 2 (3, resp.) will be denoted as  1 ( 2 , resp.).

The singularity analysis of this manipulator [START_REF] Liu | Kinematics, Singularity and Workspace of Planar 5R Symmetrical Parallel Mechanism[END_REF] shows that the Type 2 singularities appear when links 2 and 3 are parallel (see also Fig. 2 in [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF]). In both cases, the gained degree of freedom is an infinitesimal translation perpendicular to the links 2 and 3. However, if L 2 = L 3 , the gained degree of freedom may become a finite rotary motion.

Taking into account that the gravity is directed along z axis (perpendicular to the plane of motions), the expression of the potential energy V may be written as:

      2 2 2 1 2 2 2 2 2 1 1 1 ) ( ) ( 5 . 0         k q q q q k V a e a e
. [START_REF] Hunt | Special Configurations of Robot-Arms via Screw Theory[END_REF] The expression of the kinetic energy is

          2 2 2 1 2 2 5 . 0 a a a q q I y x m T         , ( 10
)
where m is the mass of the payload and x  and y  are the velocities of point C along x and y axes, respectively.

Thus the dynamic model can be obtained from ( 3) and (4):

0 W A C p         T T k 2 1 2   ,        y x m     p W (11a) 0 W A B p           T T a e a e n q q n q q k / / 2 2 1 1 1 (11b)                 n q q n q q n k q q I a e a e a a a / / 2 2 1 1 1 2 1      (11c)
Matrices A, B and C may be found from the loop closure equations: 0

) ( ) sin ( ) cos ( 2 1 2 2 1 1 2 1 1 1          L q L y q L a x f e e (12a) 0 ) ( ) sin ( ) cos ( 2 2 3 2 2 4 2 2 4 2          L q L y q L a x f e e
(12b) from which it comes:

                             e e e e i q L y a q L x q L y a q L x a a a a f
L q a q a L f q B (14)                           3 2 2 3 1 2 0 0 2 0 0 2 L L L L f i   ε C (15) with   T y x,  x ,   T e e e q q 2 1 ,  q and   T 2 1 ,   ε
. As the trajectory x(t) is known, from (11a), we can obtain the values of . Introducing it into (12), the parameters e q are computed. Finally, from (11b), the values of the actuated variables a q are obtained. Differentiating it two times with respect to the time allows obtaining the values of the input torques.

Analyzing these expressions, it could be observed that the dynamic model depend on the position, velocity and acceleration of point C, but also of the jerk and its first derivative. Therefore, a ninth-degree polynomial has to be applied as a control law when the end-effector is not in singular configuration.

In order to avoid infinite values of the input torques when crossing a Type 2 singularity, Eq. ( 8) has to be satisfied. From matrix A, one can find that the twist of the infinitesimal displacement in the singularity can be written under the form:

T ] cos , sin [ 1 1     s t (16)
Thus, the examined manipulator can pass through the given singular positions if the wrench W p determined by ( 11) is orthogonal to the direction of the uncontrollable motion t s described by [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF].

Let us now consider the motion planning, which makes it possible to satisfy this condition. For this purpose the following parameters of manipulator's links are specified: With regard to the prescribed trajectory generation, the point C should reproduce a motion along a straight line between the initial position C 0 (x 0 , y 0 ) = C 0 (0.1, 0.345) and the final point C f (x f , y f ) = C f (-0.05, 0.195) in t f = 1.5 s (Fig. 2). Thus, the given trajectory can be expressed as follows: However, the manipulator will pass by a Type 2 singular position at point C s (x s , y s ) = C s (0, 0.245) (Fig. 2). Developing the condition (8) for passing through the singular position for the planar 5R parallel manipulator at point C s , we obtain: 0  y  (18) Then, taking into account that the velocity and the acceleration of the end-effector in initial and final positions are equal to zero, the following thirteen boundary conditions are found:

a = 0.2 m, L 1 = L 2 = L 3 = L 4 = 0.
                  ) ( ) ( ) ( ) ( ) ( ) ( 0 0 0 0 y y t s y x x t s x t y t x f f x (17)
s (t 0 ) = 0, ( 19
)
s (t f ) = 1, (20)
s (t s = 1 s) = 2/3, (21) 0 ) ( 0  t s  , (22) 0
) (  f t s  , (23) 0 )
(  s t s  (in our case 5 . 1 ) (  s t s  ), (24) 0 )
( 0  t s  , (25) 0
) (  f t s  , (26) 0
) (  s t s  . ( 27
) 0 ) ( 0  t s   , (28) 0
) (  f t s   , ( 29
) 0 / )) ( ( 0  dt t s d    , (30) 0
/ )) ( (  dt t s d f    , (31)
Thus the generation of the motion by the obtained twelfth order polynomial makes it possible to pass through the singularity without perturbation and the input torques remain in the limits of finite values.

In order to compare the different cases of trajectory planning, in Figs. 3 and4 are given the values of the input torques obtained using the software ADAMS for the following numerical simulations: A: a trajectory between points C 0 and C' f (x' f , y' f ) = C' f (-0.05, 0.295) (Fig. 2) without meeting any singularity. For such a case, the following ninth order polynomial law is used s(t)= 16.59 t 5 -36.87 t 6 + 31.6 t 7 -12.29 t 8 + 1.82 t 9 for the trajectory planning out of the singular zone of the manipulator. In this case the values of the input torques are finite. Thus, the numerical simulations show that the obtained optimal dynamic conditions assume the passing of the flexible manipulator through the singular position.

CONCLUSION

In our previous work, we have shown that any parallel manipulator can pass through the singular positions without perturbation of motion if the wrench applied on the endeffector by the legs and external efforts is orthogonal to the twist along the direction of the uncontrollable motion [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF]. This condition was applied to the rigid-link manipulators. The obtained results showed that the planning of motion for assuming the optimal force generation can be carried out by a eight order polynomial law. In our other study [START_REF] Briot | On the Dynamic Properties of Rigid-Link Flexible-Joint Parallel Manipulators in the Presence of Type 2 Singularities[END_REF] the rigidlink flexible-joint manipulators have been studied. It was shown that the degree of the polynomial law should be different, when the flexibility of actuated joints is introduced. The obtained results disclosed that the planning of motion for assuming the optimal force generation in the rigid-link flexible-joint manipulators must be carried out by a twelfth order polynomial law.

In this paper, we have expanded the information about the dynamic properties of parallel manipulators in the presence of Type 2 singularity by including in the studied problem the link flexibility and the payload. The obtained results have shown that the planning of motion for assuming the optimal force generation in the structurally flexible manipulators with payload must be also carried out by a twelfth order polynomial law.

The suggested technique was illustrated by a 5R planar parallel manipulator. The obtained dynamic properties have been validated by numerical simulations carried out using the software ADAMS.

  the presence of a Type 2 singularity, it is possible to satisfy conditions[START_REF] Zlatanov | Con-straint Singularities of Parallel Mechanisms[END_REF] if the wrench applied on the platform by the legs and external efforts

Figure 1 .

 1 Figure 1. Kinematic chain of the planar 5R parallel manipulator.

  25 m; m = 10 kg; k 1 = 25000 Nm/rad; k 2 = 350000 N/m. Moreover, we use the physical parameters of a real Harmonic Drive® system (HDUC-1U-11) coupled with an actuator Gamdrive (11-50-MB-SP-024-CR), that are I a = 0.064.10 -4 kg.m² and n = 50.

Figure 2 .

 2 Figure 2. Initial, singular and final positions of the planar 5R parallel manipulator.

Figure 3 .

 3 Figure 3. Torques values for the actuator 1.

Figure 4 .

 4 Figure 4. Torques values for the actuator 2.


Sébastien Briot
email: sbriot@yahoo.fr

G C G M Ea3913

Vigen Arakelian
email: vigen.arakelian@insa-rennes.fr

On the Dynamic Properties of Rigid-Link Flexible-Joint Parallel Manipulators in the Presence of Type 2 Singularities

, the dynamic properties of rigid-link parallel manipulators, in the presence of Type 2 singularities, have been studied. It was shown that any parallel manipulator can pass through the singular positions without perturbation of motion if the wrench applied on the end-effector by the legs and external efforts of the manipulator are orthogonal to the twist along the direction of the uncontrollable motion. This condition was obtained using symbolic approach based on the inverse dynamics and the study of the Lagrangian of a general rigid-link parallel manipulator. It was validated by experimental tests carried out on the prototype of a four-degrees-of-freedom parallel manipulator. However, it is known that flexibility of the mechanism may not always been neglected. Indeed, joint flexibility is the main source contributing to overall manipulator flexibility and it leads to the trajectory distortion. Therefore, in this paper, the condition for passing through a Type 2 singularity of parallel manipulators with flexible joints is studied. The suggested technique is illustrated by the example of a 5R parallel manipulator with flexible joints. It is shown that passing through singularity is possible if the twelfth order polynomial trajectory planning is applied. The obtained results are validated by numerical simulations carried out using the ADAMS software.

Introduction

Over the past decades, with the large development of parallel manipulators, more attention has been paid to their kinematic, kinetostatic and dynamic properties, and in particular, to their singularities. Several papers deal with singularity analysis of parallel manipulators [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF][START_REF] Zlatanov | Singularity Analysis of Mechanisms and Robots via a Velocity-Equation Model of the Instantaneous Kinematics[END_REF][START_REF] Zlatanov | Constraint Singularities of Parallel Mechanisms[END_REF][START_REF] Hunt | Special Configurations of Robot-Arms via Screw Theory[END_REF][START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF][START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF][START_REF] Briot | Singularity Analysis of PAMINSA Manipulator[END_REF][START_REF] Gosselin | The Optimum Design of Robotic Manipulators using Dexterity Indices[END_REF][START_REF] Merlet | Jacobian, Manipulability, Condition Number, and Accuracy of Parallel Robots[END_REF][START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF][START_REF] Alba-Gomez | Consistent Kinetostatic Indices for Planar 3-DOF Parallel Manipulators, Application to the Optimal Kinematic Inversion[END_REF][START_REF] Arakelian | Increase of Singularity-Free Zones in the Workspace of Parallel Manipulators using Mechanisms of Variable Structure[END_REF][START_REF] Hubert | Singularity Analysis through Static Analysis[END_REF]. Most of them present the analysis of singular configurations from a kinematic point of view [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF][START_REF] Zlatanov | Singularity Analysis of Mechanisms and Robots via a Velocity-Equation Model of the Instantaneous Kinematics[END_REF][START_REF] Zlatanov | Constraint Singularities of Parallel Mechanisms[END_REF][START_REF] Hunt | Special Configurations of Robot-Arms via Screw Theory[END_REF][START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF][START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF][START_REF] Briot | Singularity Analysis of PAMINSA Manipulator[END_REF]. Algebraically, a singularity analysis is based on the degeneracy of Jacobian matrices of the mechanical structure, or of the system of reciprocal screws (wrenches) applied to the platform by the legs. However, it is also known that, when parallel manipulators have Type 2 singularities [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF], they lose their stiffness and their quality of motion transmission, and as a result, their payload capability. Therefore, the singularity zones in the workspace of manipulators may be analyzed not only in terms of kinematic criterions, from the theoretically perfect model of manipulators, but also in terms of kinetostatic approaches [START_REF] Gosselin | The Optimum Design of Robotic Manipulators using Dexterity Indices[END_REF][START_REF] Merlet | Jacobian, Manipulability, Condition Number, and Accuracy of Parallel Robots[END_REF][START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF][START_REF] Alba-Gomez | Consistent Kinetostatic Indices for Planar 3-DOF Parallel Manipulators, Application to the Optimal Kinematic Inversion[END_REF][START_REF] Arakelian | Increase of Singularity-Free Zones in the Workspace of Parallel Manipulators using Mechanisms of Variable Structure[END_REF][START_REF] Hubert | Singularity Analysis through Static Analysis[END_REF].

Moreover, while it is demonstrated using the kinetostatic approach that, when subjected to Type 2 singularities parallel manipulator lock up, it has been shown experimentally that, via optimal dynamic control of manipulators, it is possible to pass through these singular zones. Thus, it is evident that singular configurations should also be examined in terms of dynamic aspects.

The further study of singularity in parallel manipulators has revealed an interesting problem that concerns the path planning of parallel manipulators under the presence of singular positions, i.e. the motion feasibility in the neighborhood of singularities. In this case the dynamic conditions can be considered in the design process. One of the most evident solutions for the stable motion generation in the neighborhood of singularities is to use redundant sensors and actuators [START_REF] Alvan | On the Control of the Spatial Parallel Manipulators with Several Degrees of Freedom[END_REF]- [START_REF] Kotlarski | Singularity Avoidance of a Kinematically Redundant Parallel Robot by a Constrained Optimization of the Actuation Forces[END_REF]. However, it is an expensive solution to the problem because of the additional actuators and the complicated control of the manipulator caused by actuation redundancy. Another approach concerns with motion planning to pass through singularity [START_REF] Bhattacharya | Comparison of an Exact and Approximate Method of Singularity Avoidance in Platform Type Parallel Manipulators[END_REF]- [START_REF] Perng | Inverse Kinematic Solutions for a Fully Parallel Robot with Singularity Robustness[END_REF], i.e. a parallel manipulator may track a path through singular poses if its velocity and acceleration are properly constrained. This is a promising way for the solution of this problem. However, only a few research papers on this approach have addressed the path planning for obtaining a good tracking performance. But they have not adequately addressed the physical interpretation of dynamic aspects.

In our recent work [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF], optimal force generation in parallel manipulators for passing through the singular positions has been studied. It was shown that any parallel manipulator can pass through the singular positions without perturbation of motion if the wrench applied on the end-effector by the legs and external efforts of the manipulator are orthogonal to the twist along the direction of the uncontrollable motion. This paper was concerned with the study of rigid-link parallel manipulators without any flexibility. However, it should be noted that several factors may bring a loss of rigidity in parallel manipulators (elasticity of links, clearance in joints, etc.). But their contribution can be considerably reduced for a properly designed and constrained mechanical system. However, even for the most optimum design of manipu-lator one main source of flexibility remains and it cannot be easily reduced: it is the flexibility in actuated joints, due to the use of Harmonic Drive® systems.

Therefore, in the present work the dynamic condition for passing through the singular positions is defined in general for parallel manipulators taking into account the elasticity in the actuated joints.

The paper is organized as follows. The next section presents theoretical aspects of the examined problem. As in our previous work, using the Lagrangian formulation, the condition of force distribution is defined, that allows the passing of any parallel manipulator through the Type 2 singular positions. In section 3, the suggested solution is illustrated via 5R planar parallel manipulator. In section 4, the conclusions are given.

Optimal dynamic conditions for passing through Type 2 singularity

Let us consider a parallel manipulator of m links, n degrees of freedom and driven by n actuators.

The general Lagrangian dynamic formulation for a manipulator with elasticity in actuated joints can be expressed as [START_REF] Spong | An Integral Manifold Approach to the Feedback Control of Flexible Joint Robots[END_REF]:

d d L L q q 0                dt d , ( 1a
) a a L L q q τ                dt d , (1b)
where, -L is the Lagrangian of the manipulator; L = T-V, where T is the kinetic energy and V the potential energy due to gravitational forces, friction and elasticity; -

T a n a a a q q q ] ,..., , [ 2 1  q and T a n a a a q q q ] ,..., , [ 2 1      q
represent the vectors of position and velocity of the actuators, respectively; -

T d n d d d q q q ] ,..., , [ 2 1
 q and

T d n d d d q q q ] ,..., , [ 2 1      q
represent the vectors of position and velocity of controlled links, respectively, i.e. the position and velocity of the links that are controlled by the displacement of the actuators in which there are elasticity; -  is the vector of the actuators efforts. However, for a parallel mechanism, the position (velocity, resp.) of the end-effector is a non-trivial function of the position (velocity, resp.) of the controlled links, therefore it is preferable to rewrite Eq. (1) using the Lagrange multipliers, as follows:

λ B W 0 b T   , d d L L q q W b                dt d (2a) a a L L q q τ                dt d , ( 2b
)
where  is the Lagrange multipliers vector, which is related to the wrench applied on the platform by:

p W A T    ,                     x v W p L L dt d ( 3
)
where, - A and B are two matrices relating the vectors v and q  according to

d q B Av  
. They can be found by differentiating the closure equations with respect to time;

- T z y x ] , , , , , [     x and T z y x ] , , , , , [           v
are trajectory parameters and their derivatives, respectively; x, y, z represent the position of the controlled point and  and  the rotation of the platform about three axes a   a  and a  ; - W p is the wrench applied on the platform by the legs and external forces [START_REF] Khalil | A Novel Solution for the Dynamic Modeling of Gough-Stewart Manipulators[END_REF] expressed along axes a   a  and a  .

Expressing W p in the base frame, one can obtain:

p R b W J W 0 0 T   , ( 4a
) a a L L q q τ                dt d , (4b)
where

  B A J R 1 0  
is the Jacobian matrix between twist t of the platform (expressed in the base frame) and d q  , D A A R  0 is the expression of matrix A in the base frame, where D is a transformation matrix, of which expression is given in [START_REF] Merlet | Parallel Robots[END_REF].

For any prescribed trajectory x(t), the values of vectors d q   , d q  and q d can be found using the inverse kinematics. Thus, taking into account that the manipulator is not in a Type 1 singularity [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF], the terms W b and p R W 0 can be computed. However, for a trajectory passing through a Type 2 singularity, the determinant of matrix J is indefinite. Numerically, the values of the efforts applied by the actuators become infinite. In practice, the manipulator either is locked in such a position of the end-effector or it can not follow the prescribed trajectory.

It is known that a Type 2 singularity appears when the determinant of matrix A R 0 vanishes, in other words, when at least two of its columns are linearly dependant [START_REF] Merlet | Parallel Robots[END_REF]. So, one may obtain such a relationship:

   6 1 j j j 0 A  , (5)
where A j represents the j-th column of matrix A R 0 and  j are coefficients, which in general can be functions of d p q (p = 1, …, n). It should be noted that the vector t s =      2 , …,  6 ] T represents the direction of the uncontrollable motion of the platform in a Type 2 singularity. By substituting ( 5) into (3), we obtain

j T j W  λ A , j = 1, …, 6 (6)
where W j is the j-th row of vector p R W 0 .

Then, from ( 5) and ( 6) the following conditions are derived:

          6 1 6 1 0 j j j j T j j W   λ A . ( 7
)
The right term of Eq. ( 7) corresponds to the scalar product of vectors t s and p R W 0 . Thus, in the presence of a Type 2 singularity, it is possible to satisfy conditions [START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF] if the wrench applied on the platform by the legs and external efforts p R W 0 are orthogonal to the direction of the uncontrollable motion t s . Otherwise, the dynamic model is not consistent. Obviously, in the presence of a Type 2 singularity, the displacement of the end-effector of the manipulator has to be planned to satisfy [START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF]. Therefore, our task will be to achieve a trajectory which will allow the manipulator to pass trough the Type 2 singularities, i.e. which will allow the manipulator respecting condition [START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF].

In the dynamic model of the rigid-link flexible-joint manipulator [START_REF] Spong | An Integral Manifold Approach to the Feedback Control of Flexible Joint Robots[END_REF], the efforts  applied on the actuators may be expressed as follows:

) ( ) , ( ) ( dt d d d d d d d d d L L q D q q q C q q M q q τ                       , (8)
This equation only depends on the acceleration, velocity and position of the actuators of the manipulator. Therefore, in order to avoid some discontinuity on the efforts , the polynomial used for the trajectory planning should be at least of degree 5 (because initial and final positions are known, and the velocities and accelerations at the beginning and the end of the trajectory should be equal to 0). In our previous work [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF], it was shown that the condition for passing through the Type 2 singular configurations added 3 supplementary conditions, and therefore the polynomial used for the trajectory planning should be at least of degree 8 (the position, the velocity and the acceleration when passing through the singularity should be constrained).

In the present study, it will be shown that the degree of the polynomial law should be different, when taking into account the flexibility on the actuating system. Indeed, it will be presented further that the efforts  applied on the actuators depend as previously on the position, the velocity and the acceleration of the actuators, but also on the jerk and its first derivative. So, due to the addition of elasticity in the actuated joints, the polynomial used for the trajectory planning should be at least of degree 12.

In the next section, an example illustrates the obtained results discussed above. This example presents a planar 5R parallel manipulator, which allows obtaining relatively simple symbolic model for demonstrating the expected results by numerical simulations. The results are validated using the ADAMS software.

Illustrative example

In the planar 5R parallel manipulator, as shown in Fig. 1, the output point is connected to the base by two legs, each of which consists of three revolute joints and two links. In each of the two legs, the revolute joint connected to the base is actuated. Thus, such a manipulator is able to position its output point in a plane. As shown in Fig. 1, the input joints are denoted as A and E with input parameters d q 1 and d q 2 . The common joint of the two legs is denoted as C, which is also the output point with con-trolled parameters x and y. A fixed global reference system xOy is located at the center of AE with the y-axis normal to AE and the x-axis directed along AE. The lengths of the links AB, BC, CD, DE are respectively denoted as L 1 , L 2 , L 3 and L 4 . The positions of the centers of masses S i of links from joint centers A, B, D and E are respectively denoted by dimensionless lengths r 1 , r 2 , r 3 and r 4 , i.e. . Actuators 1 and 2 are connected to links 1 and 4, respectively, via harmonic drive systems which are presented by a model similar to that given in [START_REF] Spong | An Integral Manifold Approach to the Feedback Control of Flexible Joint Robots[END_REF]. The position of actuator i is denoted as a i q . It is assumed that the actuator i is capable to deliver a couple  i to the motor shaft, which is elastically coupled to the link j of the robot (i = 1 or 2, j = 1 or 4). The flexibility is modeled by a torsion spring with stiffness k i . The gear ratio is denoted n i . In the following of this paper, n 1 = n 2 = n and k 1 = k 2 = k. a i I is the axial moment of inertia of the motor i plus the Harmonic drive system.

The singularity analysis of this manipulator [START_REF] Liu | Kinematics, Singularity and Workspace of Planar 5R Symmetrical Parallel Mechanism[END_REF] shows that the Type 2 singularities appear when legs 2 and 3 are parallel (see also Fig. 2 in [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF]).

In both cases, the gained degree of freedom is an infinitesimal translation perpendicular to the legs 2 and 3. However, if L 2 = L 3 , the gained degree of freedom may become a finite rotary motion.

In order to simplify the analytic expressions, we consider that the gravity effects are along the z-axis and consequently the input torques are only due to inertia effects. We also admit that there is no friction in the system. To simplify the computation, it is also preferable to replace the masses of moving links by concentrated masses [START_REF] Seyferth | Massenersatz duch Punktmassen in Räumlichen Getrieben[END_REF][START_REF] Wu | On the Dynamic Balancing of Multi-DOF Parallel Mechanisms with Multiple Legs[END_REF]. For a link i with mass m i and its axial moment of inertia I i , we have:

                                 i i i i i i i i i i i I m m m m L r L r r r 0 ) 1 ( 0 1 0 1 1 1 3 2 1 2 2 2 2
, (i = 1, 2, 3, 4) [START_REF] Gosselin | The Optimum Design of Robotic Manipulators using Dexterity Indices[END_REF] where m ij (j = 1, 2, 3) are the values of the three point masses placed at the centers of the revolute joints and at the center of masses of the link i.

In this case, the potential energy V can be written as:

   n n k V a d T a d / / 2 q q q q    , (10)
where

  T d d d q q 2 1 ,  q ,   T a a a q q 2 1 ,  q
, and the kinetic energy T as:

a T a a D D C C B B S S S S S S S S I m m m m m m m T q q V V V V V V V   2 1 ) ( 2 1 2 2 2 2 4 4 2 3 3 2 2 2 2 1 1         (11)
where, . The terms m ij (i = 1, 2, 3, 4) are deduced from the relation (9), V Si is the vector of the linear velocities of the center of masses S i ; V B , V C and V D are the vectors of the linear velocities of the corresponding axes.

Thus the dynamic model can be obtained from (2):

p b W J W 0 T R 5   (12a)   n n k I a d a a / q q q       (12b)
taking into account that for examined manipulator:

  n k a d / * q q W W b b    , D D B B b F J F J W T T   * (13) where,        0 cos 0 sin 1 1 1 1 q L q L B J ,         2 4 2 4 cos 0 sin 0 q L q L D J , ( 14
) C B B F   1 1 C B m m   , C D D F   3 2 C D m m   , (15)
                      1 1 2 1 1 1 1 1 sin cos cos sin q q q q q q L    B  , ( 16a
)                       2 2 2 2 2 2 2 4 sin cos cos sin q q q q q q L    D  ,        y x     C  (16b) 2 2 2 2 1 1 1 ) 1 ( r m m r m m S B S B     , ) 1 ( 2 2 2 1 r r m m S C   , ( 17
) ) 1 ( 3 3 3 3 r r m m S C   , 2 3 3 2 4 4 2 ) 1 ( r m m r m m S D S D    
. ( 18) The term W p is given by:

D C B p W    3 2 1 C C C m m m    , ( 19
)
                    2 4 2 4 1 1 1 1 22 21 12 11 5 sin cos sin cos 2 q L y a q L x q L y a q L x a a a a R A , ( 22
)           ) cos sin ( 0 0 ) cos sin ( 2 22 2 21 4 1 12 1 11 1 5 q a q a L q a q a L R B . ( 23
)
For a given trajectory q d , the vector of the positions of the actuators q a can be deduced from (12a), as well as the acceleration of the actuators a q   :

  d T R a n k n q W J W q p b     5 * (24a)   d T R a n dt d k n q W J W q p b         5 * 2 2 (24b)
Introduction Eq. ( 24) into Eq. (12b), one can deduce the vector of actuator torques :

    p b p b W J W q W J W T R d T R a n n dt d k n I 5 * 5 * 2 2 1                  . (25)
Analyzing this expression, it could be observed that, as terms * b W and p W depend on the position, velocity and acceleration of the input links, the input torques depend not only on these parameters, but also on the jerk and its first derivative. Therefore, on the contrary of rigid manipulators for which, in order to avoid discontinuities in the input torques, a fifth-degree polynomial is sufficient as a control law when the end-effector is not in singular configuration, for non rigid robots, the degree of the polynomial should be increased (indeed, it should be at least a ninth-degree polynomial).

In order to avoid infinite values of the input torques when crossing a Type 2 singularity, Eq. ( 7) has to be satisfied. From matrix A 5R , one can find that the twist of the infinitesimal displacement in the singularity can be written under the form:

T ] cos , sin [ 1 1     s t (26)
Thus, the examined manipulator can pass through the given singular positions if the wrench W p determined by ( 19) is orthogonal to the direction of the uncontrollable motion t s described by [START_REF] Spong | An Integral Manifold Approach to the Feedback Control of Flexible Joint Robots[END_REF].

Let us now consider the motion planning, which makes it possible to satisfy this condition. For this purpose the following parameters of manipulator's links are specified: L 1 = L 2 = L 3 = L 4 = 0.25 m; r 1 = r 2 = r 3 = r 4 = 0.5; a = 0.2 m; m 1 = m 4 = 2.81 kg; I 1 = I 4 = 0.02 kg.m 2 ; m 2 = m 3 = 1.41 kg; I 2 = I 3 = 0.01 kg.m 2 ; I a = 0.067 kg/m 2 ; k = 250 Nm/rad; n = 50.

With regard to the prescribed trajectory generation, the point C should reproduce a motion along a straight line between the initial position C 0 (x 0 , y 0 ) = C 0 (0.1, 0.345) and the final point C f (x f , y f ) = C f (-0.1, 0.145) in t f = 2 s (Fig. 2). Thus, the given trajectory can be expressed as follows:

                    ) ( ) ( ) ( ) ( ) ( ) ( 0 0 0 0 y y t s y x x t s x t y t x f f x (27)
However, the manipulator will pass by a Type 2 singular position at point C s (x s , y s ) = C s (0, 0.245) (Fig. 2).

Developing the condition for passing through the singular position [START_REF] Bonev | Singularity Analysis of 3-DOF Planar Parallel Mechanisms via Screw Theory[END_REF] for the planar 5R parallel manipulator at point C s , we obtain:

0 6 3 ) 48 248 ( 2 2 2 1 1    y m y x L m C C    
(28) Then, taking into account that the velocity and the acceleration of the end-effector in initial and final positions are equal to zero, the following thirteen boundary conditions are found:

s (t 0 ) = 0, (29) s (t f ) = 1, (30) s (t s = 1 s) = 0.5, (31) 0
) ( 0  t s  , (32) 0
) (  f t s  , (33) 1
) /( ) /( ) ( 0 0      x x x y y y t s f s f s s    , ( 34
) 0 ) ( 0 0   s t s     , (35) 0
) (   f f s t s     , ( 36
) ) 6 ) ( 3 /( ) 48 248 ( ) ( 2 0 2 2 1 1 C f s s C s s m y y y x L m s t s           , (37) 0 ) ( 0 0   s t s       , (38) 0
) (   f f s t s       , ( 39
) 0 ) ( ) 4 ( 0 0   s t s dt d    , ( 40
) 0 ) ( ) 4 (   f f s t s dt d    . ( 41
)
Table 1. Variation of the input torques as a function of the polynomial law used for the trajectory. From ( 28)-(41), the following twelfth order polynomial trajectory planning is found:   Thus the generation of the motion by the obtained twelfth order polynomial makes it possible to pass through the singularity without perturbation and the input torques remain in the limits of finite values.

In order to compare the different cases of trajectory planning, in table 1 are given the values of the input torques obtained using the software ADAMS for the following numerical simulations: A1: a trajectory between points C 0 and C' f (x' f , y' f ) = C' f (-0.1, 0.345) (Fig. 2) without meeting any singularity. For such a case, the following fifth order polynomial law is used   for the trajectory planning of the rigid-link flexible-joint manipulator inside the singular zone. The numerical simulation shows that the values of the input torques close to the singular positions tend to infinity. D:

the twelfth order polynomial law (42) for the trajectory planning of the rigid-link flexible-joint manipulator inside of the singular zone. The values of the input torques are finite and there are no discontinuities. Thus, the numerical simulations show that the obtained optimal dynamic conditions assume the passing of the rigid-link flexiblejoint manipulator through the singular position.

Conclusion

At a singular configuration, in the case of an arbitrary generation of forces, a manipulator may not reproduce stable motion with prescribed trajectory. Nevertheless it is approved that there are several motion planning techniques, which allow passing through these singular zones. These approaches are simulated by numerical examples and illustrated on several parallel structures. However, in these studies much more attention was focused only on control aspects of this problem and little attention has been paid to the dynamic interpretation, which is a crucial factor for governing the behavior of parallel manipulators at the singular zones.

In our previous work [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF], the dynamic properties of parallel manipulators in the presence of Type 2 singularity have been studied. It was shown that any parallel manipulator can pass through the singular positions without perturbation of motion if the wrench applied on the end-effector by the legs and external efforts of the manipulator are orthogonal to the twist along the direction of the uncontrollable motion. This condition was applied to the rigid-link manipulators without clearance or flexibility in the joints. The obtained results showed that the planning of motion for assuming the optimal force generation can be carried out by a eight order polynomial law.

In the present paper the rigid-link flexible-joint manipulators have been studied. It was shown that the degree of the polynomial law should be different, when the flexibility of actuated joints is introduced into conditions of the optimal force generation in the presence of singularity. The obtained results disclosed that the planning of motion for assuming the optimal force generation in the rigid-link flexible-joint manipulators must be carried out by a twelfth order polynomial law. The suggested technique was illustrated by an example, which presents 5R planar parallel manipulator with flexible joints. The numerical simulations carried out using the software ADAMS validated the obtained theoretical results.

Fig. 1 .

 1 Fig. 1. Kinematic chain of the planar 5R parallel manipulator.

  and the Jacobian matrix J 5R by:

Fig. 2 .

 2 Fig. 2. Initial, singular and final positions of the planar 5R parallel manipulator.

  0 and C f inside the singular zone for the rigid-link manipulator without taking into account the flexibility in the actuated joints. In this case the values of the input torques close to the singular positions tend to infinity. B: the eight order polynomial law s(t) = -0.25851 t 3 + 3.84228 t 4 -5.72792 t 5 + 3.58909 t 6 -1.07101 t 7 + 0.12606 t 8 for the trajectory planning of the rigid-link manipulator without flexibility in the actuated joints inside the singular zone. The obtained results show that the values of the input torques are finite.

			s	t		. 1	25	t	3		9375 . 0	t	4		1875 . 0	t	5	for the trajec-
		tory planning out of the singular zone of the rigid-link
		manipulator without taking into account the flexibility
		in the actuated joints. In this case the values of the input
		torques are finite.
	A2:	s	the   . 1 t 	25	t	3	fifth 9375 . 0 	t	4	order 1875 . 0 	t	5	polynomial for the trajectory law
	planning between C C: the ninth   5 56 . 6 94 . 3 t t t s  	6	order 22 . 4 	t	7		polynomial 8 14 . 0 23 . 1 t 	t	9	law


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

Wisama Khalil
email: wisama.khalil@irccyn.ec-nantes.fr

Recursive and Symbolic Calculation of the Stiffness and Mass Matrices of Parallel Robots

Keywords: Parallel robots, Flexibilities, Natural frequencies

This paper presents a symbolic and recursive calculation of the stiffness and mass matrices of parallel robots. In order to reduce the computational time required for simulating the elastodynamic behavior of robots, it is necessary to minimize the number of operators in the symbolic model expression. Some algorithms have been proposed for the rigid case or for parallel robots with lumped springs. In this paper, we extend the previous works to parallel robots with distributed flexibilities. The proposed algorithm, that takes advantage of recursive calculations for the computation of the Jacobian matrices defining the kinematic constraints, is used to compute the natural frequencies of a robot developed at IRCCyN: the NaVARo.

Introduction

The large computational time required for calculating the natural frequencies of a robot prevents to use them in many applications, such as real-time control, design optimization process, etc. To decrease the computational cost, this paper focuses on the efficient symbolic computation for the stiffness and mass matrices of flexible parallel robots. This approach could be combined with model reduction methods [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF][START_REF] Craig | Structural dynamics[END_REF][START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF], but this is not the main goal of that paper.

For the computation of the robot natural frequencies, two main approaches are generally proposed (see [START_REF] Dwivedy | Dynamic analysis of flexible manipulators, a litterature review[END_REF] for a large literature review): (i) lumped modeling [START_REF] Khalil | Modeling of mechanical systems with lumped elasticity[END_REF][START_REF] Wittbrodt | Dynamics of Flexible Multibody Systems[END_REF] and (ii) modeling using distributed flexibilities [START_REF] Bauchau | Flexible multibody dynamics[END_REF][START_REF] Boyer | An efficient calculation of the flexible manipulator inverse dynamics[END_REF][START_REF] Cammarata | An algorithm to study the elastodynamics of parallel kinematic machines with lower kinematic pairs[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Stachera | Automation and Robotics, chap. 15: Derivation and Calculation of the Dynamics of Elastic Parallel Manipulators[END_REF]. The lumped modeling is generally simpler to use by non-experts in finite element meth-

M ij (M 0ij ) M 0ij B 0ij B ij (B 0ij ) A 0ij =A ij F b0 F ij v ij (A ij ) ω ij (A ij ) u eij =Σ φ dkij (M 0ij )q ekij k=1 Nij undeformed config.
deformed configuration (c) A single flexible link Fig. 1 A general parallel robot composed of flexible elements. ods but, to obtain a correct model accuracy, higher number of elements is required, thus increasing the computational time.

Contrary to lumped modeling, using distributed flexibilities allows the improvement of the model accuracy. However, such methods require highly-skilled users. In [START_REF] Bauchau | Flexible multibody dynamics[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Stachera | Automation and Robotics, chap. 15: Derivation and Calculation of the Dynamics of Elastic Parallel Manipulators[END_REF], some general methodologies are proposed. In the case of closed-loop mechanisms, some Jacobian matrices are computed that allow taking into account the kinematic dependencies. However, such general methodologies are not specifically designed for parallel robots and they do not guarantee the minimization of the number of operators for the symbolic computation of the model. To the best of our knowledge, a systematic procedure to compute the mass and stiffness matrices (using distributed flexibilities) of parallel robots with a minimal numbers of operators has never been proposed.

The present work aims at filling this gap. In order to minimize the number of operations, the Jacobian matrices defined in the principle of virtual powers (PVP) are computed using recursive algorithms. For computing the stiffness and mass matrices of parallel robots, the approach proposed in [START_REF] Briot | Recursive and symbolic calculation of the elastodynamic model of flexible parallel robots[END_REF] is adapted. It proposes to (1) convert the parallel robot into a virtual system defined by a tree-structure robot composed of the kinematic chains of the actual robot for which all joints (passive and active) are considered actuated and a free body (the platform which is considered as rigid) (Fig. 1), (2) compute the elastodynamic model of this new virtual system, and (3) finally, close the loops by using the PVP.

This method is effective, systematic, can be applied to any parallel robot.

Stiffness and Mass Matrices for the Virtual System

Let us consider a parallel robot composed of one rigid fixed base (denoted as element 0), one rigid moving platform and n legs, each leg being a serial kinematic chain composed of (m i -1) elements 1 connected by m i joints of coordinates q ik (revolute, prismatic or fixed joints -i = 1, ..., n) located at points A ik (k = 1, ..., m i -Fig. 1(a)). The j-th element of the i-th leg is denoted by i j and its displacement can be parameterized by the coordinates i j t i j2 which represents the twist of the body i j at the origin of the local frame F i j (Fig. 1(c)) and qe i j that are the generalized velocities characterizing the elastic displacement of the body i j i j t e i j (M i j ) = Φ i j (M 0i j ) qe i j [START_REF] Arakelian | Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF] where i j t e i j (M i j ) is the deformation twist due to the body elasticity that can be parameterized using truncated series of Rayleigh-Ritz shape functions Φ i j [START_REF] Boyer | Robot Manipulators. Modeling, Performance, Analysis and Control, chap. 7: Modeling and Control of Flexible Robots[END_REF].

The vector of generalized coordinates of the tree-structure is given by T , where q t i regroups all joint variables (denoted as q T p i = [q i1 . . . q im i ]) and elastic generalized coordinates q T e i = q T e i1 . . . q T e i,m i for the real i-th leg. The Lagrangian of the tree structure system can be expressed as:

q t = q T t 1 • • • q T t n
L t = 1 2 ∑ i, j i j t T i j qT e i j M i j i j t i j qe i j -q T e i j K i j q e i j ( 2
)
where M i j and K i j are the mass and stiffness matrices of the link i j whose full expressions in the most compact form are given in [START_REF] Boyer | Robot Manipulators. Modeling, Performance, Analysis and Control, chap. 7: Modeling and Control of Flexible Robots[END_REF][START_REF] Shabana | Dynamics of flexible bodies using generalized newton-euler equations[END_REF]. To express the Lagrangian as a function of q t and qt , let us express the displacement of the element i j frame located at A i j using the following equations obtained by a recursive algorithm [START_REF] Boyer | Robot Manipulators. Modeling, Performance, Analysis and Control, chap. 7: Modeling and Control of Flexible Robots[END_REF]

j t i j = J t i j qt with J t i j = i j T i( j-1) J t i( j-1) + Φ q ei j + A i j and (3)
Φ q ei j = 0 • • • i j R i( j-1) Φ i( j-1) (A i j ) • • • 0 , A i j = 0 • • • i j a i j • • • 0 (4)
where a i j is the unit twist describing the joint i j axis [START_REF] Boyer | Robot Manipulators. Modeling, Performance, Analysis and Control, chap. 7: Modeling and Control of Flexible Robots[END_REF] and, in Φ q ei j (A i j , resp.), the term i j R i( j-1) Φ i( j-1) (A i j ) ( i j a i j , resp.) is located at the columns corresponding to the variables qe i( j-1) ( qi j , resp.). Moreover, in the previous expressions,

i j R i( j-1) = i j R i( j-1) 0 0 i j R i( j-1) , i j T i( j-1) = i j R i( j-1) I 3 -i( j-1) ri( j-1) (A i j ) 0 I 3 (5)
where I 3 is the (3 × 3) identity matrix, i j R i( j-1) is the rotation matrix between frames F i j and F i( j-1) , i( j-1) ri( j-1) (A i j ) is the cross-product matrix associated with the vector i( j-1) r i( j-1) (A i j ), i.e the position of point

A i j ≡ B i( j-1) in F i( j-1) (Fig. 1(c)).
Finally, a global Jacobian matrix J i j defined such as i j t T i j qT e i j T = J i j qt can be computed as J T i j = J T t i j O T q ei j where O q ei j is defined such that qe i j = O q ei j qt .

Introducing i j t T i j qT e i j T = J i j qt into (2) leads to L t = 1 2 ∑ i, j qT t J T i j M i j J i j qt -q T t O T q ei j K i j O q ei j q t = 1 2 qT t M t qt -q T t K t q t (6)
where M t and K t are the mass and stiffness matrices of the tree structure.

Adding the contribution of the rigid platform into (6), the Lagrangian of the total system can be written as:

L = 1 2 qt t p T M t 0 0 M p qt t p - q t x p T K t 0 0 0 q t x p = 1 2 qT g M g qg -q T g K g q g (7)
M p is the mass matrix of the rigid platform and x p represents the platform displacement (t p its twist). M g and K g are the total mass and stiffness matrices of the virtual system. q T g = q T t x T p is its vector of generalized coordinates.

Stiffness and Mass Matrices for the Parallel Robot

It is now necessary to determine one possible subset of generalized coordinates for the parallel robot. Using (3) for computing the twist i,m i t i,m i of the tip of leg i:

i,m i t i,m i = J i t i,m i qt i ( 8
)
where J i t i,m i can be obtained from J t i,m i by extracting the columns corresponding to the vector qT t i = qT p i , qT e i , i.e. the vector stacking all variables of the leg i. As the leg extremity is also linked to the rigid platform, its twist can be related to the platform twist t p via the rigid body displacement relation: i,m i t i,m i = J i p t p , where J i p = i,m i R 0 I 3 -0 pi 0 I 3 [START_REF] Craig | Structural dynamics[END_REF] in which 0 pi is the cross product matrix of vector 0 p i that characterizes the position of the attachment point A i,m i wrt the platform center (Fig. 1(a)) and i,m i R 0 is the (6 × 6) rotation matrix between the global frame and the local frame F i,m i . Thus, the following set of equations can be obtained:

   J 1 t 1,m 1 • • • 0 . . . . . . . . . 0 • • • J n t n,mn       qt 1 . . . qt n    -    J 1 p . . . J n p    t p = 0 ⇔ J t qt -J p t p = J t -J p qt t p = J g qg = 0
(10) where J g is a (r n × n q g ) matrix, n q g > r n (r = 6 for a spatial robot, r = 3 for a planar robot). This means that a subset qd of r n variables in vector qg is linked to the others. This subset is not unique. As most of parallel robots have identical legs, an idea is to put in qd the last r components q f t i of each vector qt i which can be decomposed into two parts qT

t i = q0T t i q f T t i : -     J f 1 t 1,m 1 • • • 0 . . . . . . . . . 0 • • • J f n t n,mn        q f t 1 . . . q f t n    =    J 01 t 1,m 1 • • • 0 -J 1 p . . . . . . . . . . . . 0 • • • J 0n t n,mn -J n p         q0 t 1 . . . q0 t n t p      ⇔ J f t q f t = J t p q (11)
This can be rewritten as

q f t = J f t -1 J t p q =    J d 1,1 • • • J d 1,n J d 1,n+1 . . . . . . . . . . . . J d n,1 • • • J d n,n J d n,n+1    q = J d q (12)
If the coordinates q f t i are those of the last elastic element of the leg (which is most often the case), the k-th column of matrix J f i t il corresponds to a unit twist that describes the displacement of the leg extremity due to the k-th coordinate of q f t i , i.e. [12]

J f i t il = i,m i R il -i,m i R il il pil 0 i,m i R il ( 13
)
where i,m i R il is the rotation matrix between the local frame linked at element i, m i and the local frame linked at element il and il pil is the cross product matrix of the vector il p il that characterizes the position of the leg extremity with respect to the frame linked at element il. Thus its inverse is equal to

J f i t il -1 = i,m i R T il pil i,m i R T il 0 i,m i R T il ( 14
)
which does not require much calculation. Finally, from [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] and the definition of qT g = qT t 1 , . . . , qT t n , t T p T , the matrix J defined such that qg = J q can be computed.

Introducing qg = J q into (7) leads to:

L = 1 2 qT J T M g J q -q T J T K g Jq = 1 2 qT M q -q T Kq ( 15
)
from which the natural frequencies f i (i = 1, . . . , n q g -r n) of the parallel robot can be computed as f i = eig (M -1 K)/(2π). It should be noticed that obviously, matrices M and K depend on the robot configuration.

To automatize the calculation of the mass and stiffness matrices of the robot, for each computation, the elements of a vector or a matrix containing at least one mathematical operation are replaced by an intermediate variable. This variable is written in an output file which contains the model. The elements that do not contain any operations are not modified. The obtained vectors and matrices are propagated in the subsequent equations. Consequently, at the end, the model is obtained as a set of intermediate variables. Those that have no effect on the desired output can be eliminated. This algorithm has been successfully implemented with Mathematica.

(a) Prototype of the NaVARo q 11 q 21 q 12 q 22 q 13 q 23 q 41 q 31 q 51 q 42 q 32 q 52 q 43 q 33 The NaVARo is a 3-dof planar parallel manipulator developed at IRCCyN (Fig. 2(a)) and composed of 3 identical legs and one moving platform made up of 3 segments E 1 P, E 2 P and E 3 P rigidly linked at point P. The i-th leg contains four links connected with five revolute joints in such a way that Three double clutches are mounted to the base at points A i , i = 1, 2, 3, in order to actuate either angle q 1i or angle q 2i . As a consequence, the NaVARo has 8 actuation modes [START_REF] Arakelian | Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF][START_REF] Rakotomanga | Kinetostatic performance of a planar parallel mechanism with variable actuation[END_REF]. Therefore, the platform can be moved throughout the workspace without reaching any parallel singularity thanks to a judicious actuation scheme. The kinematics of the i-th leg is described by the modified Denavit-Hartenberg parameters (MDH) [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] given in Table 1, in which γ i = π/2 if i = 1, γ i = -5π/6 if i = 2 and γ i = -π/6 if i = 3. Besides, the circumradius of the moving-platform is equal to 0.2027 m, i.e., l 5i = 0.2027 m. Each link, of rectangular cross-section, is made up of duraluminum alloy. Table 1 gives their cross-section area and inertia.

q 53 A 1 B 1 C 1 D 1 E 1 A 2 B 2 C 2 D 2 E 2 A 3 B 3 C 3 D 3 E 3 P x p
A i B i C i D i is a parallelogram link- age, i = 1,
In the experimental setup, the rotation of links 1i and 2i about point A i , i = 1, 2, 3, is locked thanks to the double clutch mechanisms. A single 3D beam element is used to model links 1i, 2i, 3i and 5i while two 3D beam elements of equal lengths are used to model links 4i, the latter being twice longer than the former.

Thus, the NaVARo is modelled as a spatial mechanism and its model contains 144 generalized coordinates, among which only 90 are independent (see Sect. 3). The model has been calculated using the proposed procedure and compiled into C code to obtain the robot natural frequencies. The computation involve the use of 36183 '+' or '-' and 37341 '*' or '/' operators, while 21383 variables are defined.

Table 1 MDH parameters of the i-th leg and characteristics of the beam cross-sections ji a( ji) σ ji γ ji b ji α ji d ji θ ji r ji 1i 0 0 γ i 0 0 d 1 = 0.4041m q 1i -γ i 0 2i 0 0 γ i 0 0 d 1 = 0.4041m q 2i -γ i 0 3i 2i 0 0 0 0 d 3 = 0.2100m q 3i 0 4i 3i 0 0 0 0 d 4 = 0.2100m q 4i 0 5i 4i 0 0 0 0 d 5 = 0.4200m q 5i 0 link A i j (m 2 ) I y i j (m 4 ) I z i j (m 4 ) I p i j (m 4 ) I 0 i j (m 4 ) 1i, 2i, 3i, 4i 2.4 • 10 -4 1.152 • 10 -8 2.000 • 10 To the best of our knowledge, there exist no works that try to minimize the number of operators in the elastodynamic models of parallel robots. Therefore, the efficiency of the proposed solution may be difficult to analyze. However, for reasons of comparison, the obtained frequencies were validated by means of an equivalent model developed using Cast3M software [START_REF]The castem software webpage[END_REF]. For the simulations, Cast3M gives the result after around 6 sec. of computation while our model send the results in around 0.01 sec. (for a Pentium 4 2.70GHz, 8Go of RAM). Both models give the same values for the first 90 natural frequencies of the NaVARo. Table 2 gives the first 5 natural frequencies of the NaVARo for the 4 robot postures shown in Fig. 3.

Conclusions

The paper has presented a symbolic and recursive calculation of the stiffness and mass matrices of parallel robots. The proposed algorithm, that takes advantage of recursive calculations for the computation of the Jacobian matrices defining the kinematic constraints, is used to compute the natural frequencies of a robot developed at IRCCyN: the NaVARo. Results have shown that the proposed model was able to give the same values as a FEA software for the first 90 natural frequencies of the NaVARo but in a considerably reduces computational time (around 0.01 sec. for our model while FEA results were obtained in several sec.).

  Fig. 2 The NaVARo

2 , 3 .

 23 The base frame F b (O, x 0 , y 0 , z 0 ) (not shown in Fig.2(b)) is defined such as point O is located at the geometric centre of the equilateral triangle A 1 A 2 A 3 . Frame F p (P, x p , y p , z p ) is attached to the moving platform. In the home configuration shown in Fig.2, F b and F p (P, x p , y p , z p ) (attached to the moving platform) coincide. (x p , y p ) are the Cartesian coordinates of point P expressed in frame F b and θ p is the orientation angle of the platform (the angle between x 0 and x p ).

8 --Fig. 3

 83 Fig.3The four test poses {x p , y p , θ p } (positions in meters, orientation in radian)

Table 2

 2 Comparison of the natural frequencies obtained with Cast3M and Matlab. Cast3M Matlab Cast3M Matlab Cast3M Matlab Cast3M Matlab Cast3M Matlab Pose (a) 44.10 44.10 44.10 44.10 53.98 53.98 60.63 60.63 95.62 95.62 Pose (b) 45.71 45.71 45.71 45.71 54.58 54.58 65.35 65.35 97.92 97.92 Pose (c) 36.98 36.98 49.31 49.31 53.37 53.37 67.28 67.28 91.80 91.80 Pose (d) 40.17 40.17 50.32 50.32 52.99 52.99 67.36 67.36 91.52 91.52

	(Hz)	f 1	f 2	f 3	f 4	f 5

Note that each robot link can be composed of one element or several elements.

In what follows, the preceeding superscript "i j" denotes that the vector is given in the frame F i j .

This work was supported by the French ANR project ARROW (ANR 2011BS3 006 01).