Julian Erskine

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

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

DETC2018-85378 WRENCH CAPABILITY ANALYSIS OF AERIAL CABLE TOWED SYSTEMS

Keywords: Quadrotor, Cable Robot, Wrench Feasibility, Capacity Margin

Aerial cable towed systems (ACTSs) can be created by joining unmanned aerial vehicles (UAVs) to a payload to extend the capabilities of the system beyond those of an individual UAV. This paper describes a systematic method of evaluating the available wrench set and the robustness of equilibrium of ACTSs by adapting wrench analysis techniques used in cable-driven parallel robots to account for the constraints of quadrotor actuation. Case studies are provided to demonstrate the analysis of different classes of ACTSs, as a means of evaluating the design and operating configurations.

INTRODUCTION

Aerial cable towed systems (ACTSs) are a class of aerial manipulator that is gaining interest as a means of payload transportation and manipulation. Currently most ACTSs are slung loads from helicopters for the purpose of military or civilian logistics, but the recent increase in the commercialization of unmanned aerial vehicles (UAVs) due to improved onboard sensing and computation, and higher energy density batteries has opened the commercial possibility of using many small UAVs to collaborate towards single tasks [START_REF] Kumar | Opportunities and challenges with autonomous micro aerial vehicles[END_REF][START_REF] Tagliabue | Collaborative transportation using mavs via passive force control[END_REF][START_REF] Tognon | Dynamics, control, and estimation for aerial robots tethered by cables or bars[END_REF]. ACTSs are composed of aerial vehicles connected to either a rigid body payload or to a point mass through cables. Systems of multiple UAVs give greater * Address all correspondence to this author. versatility in the control of the payload and expand the available task set of the system, as well as delivering an attractive degree of modularity. Many types of UAVs exist, however much of the current research is focused on quadrotors due to their high payload to weight ratio, versatile control, and simple mechanical design [START_REF] Lee | Geometric tracking control of a quadrotor uav on se(3)[END_REF].

Much of the current research into ACTSs is focused on designing trajectory tracking controllers for the suspended body. In [START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF], differential flatness is used to generate feasible trajectories for the system. The authors of [START_REF] Masone | Cooperative transportation of a payload using quadrotors: a reconfigurable cable-driven parallel robot[END_REF] use control methods adapted from cable-driven parallel robots (CDPRs) to ensure stability and trajectory convergence, while in [START_REF] Alothman | Using iterative lqr to control two quadrotors transporting a cablesuspended load[END_REF], a linear quadratic regulator is used to control the robot. In [START_REF] Manubens | Motion planning for 6-d manipulation with aerial towedcable systems[END_REF], a novel ACTS is presented and uses RRT graphs to determine feasible payload trajectories. In all these works the wrench feasibility and stability of the ACTS is considered, but there is no systematic quantification of the degree of wrench feasibility in task space.

There are strong similarities between ACTSs and CDPRs, indeed an ACTS can be considered as a CDPR with moving pulley anchor points. Wrench capabilities have been an important topic of research for CDPRs, playing an important role in their design, workspace analysis, and trajectory generation [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF][START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF][START_REF] Azizian | The dimensional synthesis of spatial cable-driven parallel mechanisms[END_REF][START_REF] Gosselin | Cable-driven parallel mechanisms: state of the art and perspectives[END_REF]. The link between CDPRs and ACTSs has been investigated by using tension distribution algorithms for the purposes of control and stability [START_REF] Masone | Cooperative transportation of a payload using quadrotors: a reconfigurable cable-driven parallel robot[END_REF][START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF]. In [START_REF] Ruiz | Arachnis: Analysis of robots actuated by cables with handy and neat interface software[END_REF], the degree of wrench inclusion is used as a performance index for CDPRs. Despite the similitudes of the two types of systems, there are differences in the modelling of the available wrench set due to the actuation limits of quadrotors. The purpose of this paper is therefore to develop a broadly applicable means of evaluating the available wrench set of ACTSs, with possible applications in design, motion planning, and for evaluating the robustness of equilibrium. In this paper, only quasi-static operations are considered, however the methods shown here may later be adapted to dynamic flight.

F 0 O 3 x 0 y 0 z 0 F p p R p F 1 F 2 F 3 u 1 u 2 u 3 t 1 t 2 t 3 B 1 B 2 B 3 A 1 A 2 A 3 v 1 v 2 v 3
Section 2 of this paper describes the parameterization of ACTSs. In section 3, the equilibrium conditions of the platform and quadrotors are formalized. The methodology used for calculating the available wrench set of the ACTS is shown in section 4 for a simple planar system. Several case studies are presented in section 5, showing the application for more complex planar and spatial systems. The paper finishes with the conclusion, list of nomenclature, and an appendix of mathematical proofs.

GEOMETRIC PARAMETERIZATION

There are different types of ACTSs so a parameterization allowing for design variability is important. From a practical perspective, it is assumed that the pose of the bodies are known with respect to an inertial frame F 0 (many experiments use multiview visual motion capture). The platform frame F p is placed at the center of mass and inertia of the payload. Its position vector p is expressed in F 0 , and in the case of rigid body loads, the orientation is expressed by the rotation matrix R p .

The ACTS is taken to have n quadrotors where the i th quadrotor has frame F i at the center of mass x i of the quadrotor, and has a rotation matrix R i defined with respect to F 0 . The point where the j th cable connects to the moving platform is called B j , and is expressed in F 0 . The cartesian coordinate vector B j is defined in F p as b j . The j th cable connects the payload to the i th quadrotor and is attached to that quadrotor at point (here assumed congruent with x i ), and the cable has a known length l j . The configuration of the ACTS is described by the set of cable vectors u 1...,m , with u i being the unit vector pointing from point B i to point A i , i = 1, . . . , m. As a function of measured states, the cable vector u i is given by Eq. [START_REF] Kumar | Opportunities and challenges with autonomous micro aerial vehicles[END_REF].

A i in space b 1 b 2 b 3 m p g t 1 u 1 t 2 u 2 t 3 u 3 m i g -t i u i f i v i
l i u i = x i -p -R p b i (1)

STATIC EQUILIBRIUM CONDITIONS OF ACTS

This paper describes the available wrench set of ACTSs in a quasi-static equilibrium situation. It is considered that an external wrench w e = [f T e , m T e ] T expressed in F 0 acts on the platform. The equilibrium equations acting on the platform and on the quadrotors are expressed separately in determining the system equilibrium, effectively dividing the systems into three components; the quadrotors, the platform, and the cables linking the two together.

Platform

Assumption 1: Cables are considered to be massless and inelastic. This means that the cables can be considered as straight lines, and as such the i th cable tension vector can be approximated as t i = t i u i .

The wrenches on the payload are generated through the sum of m tension vectors and their corresponding moments, along with gravity acting at the center of mass. There are therefore m + 1 constraining wrenches on the payload. In cases where the payload is considered a point mass, the moments generated by the cables are disregarded. As a result, the static equilibrium equations of the platform expressed in frame F 0 are expressed as follows:

m p g + m ∑ j=1 (t j u j ) + f e = 0 3 (2) m ∑ j=1 (R p b j × t j u j ) + m e = 0 3 (3)
where m p is the mass of the platform, g is the gravitational vec-tor [0, 0, -9.81] T ms -2 , and m is the total number of cables connected to the platform.

Quadrotors

The actuation wrench of the i th quadrotor is expressed in

F i as i w i = [ i f T i , i m T i ] T
, where i f i is the thrust vector [0, 0, f i ] and i m i is the moment vector [m i,x , m i,y , m i,z ] [START_REF] Lee | Geometric tracking control of a quadrotor uav on se(3)[END_REF]. The components of i w i take the form:

f i = i=4 ∑ i=1 C F ω 2 i (4a) m x,i = C F r i ω 2 1 -ω 2 3 (4b) m y,i = C F r i ω 2 2 -ω 2 4 (4c) m z,i = C M ω 2 1 -ω 2 2 + ω 2 3 -ω 2 4 ( 4d
)
where ω i is the anguar velocity of the i th propellor, r i is the distance from the center of the i th propellor to the quadrotor's center of mass, and C F and C M are the coefficients of thrust and drag of the propellors, respectively. It can be seen that a quadrotor exerts a thrust normal to the plane containing its rotors, and three independent moments about its center of inertia. This results in four degrees of actuation controlling the orientation of the unit vector v i = R i [0, 0, 1] T , and the magnitude of the thrust f i along v i . Assumption 2: Cables pass through the center of mass of the quadrotor to which they are attached. This allows the assumption that in static equilibrium, f i is the only non-null component of i w i .

The static equilibrium equations of the i th quadrotor expressed in F 0 are written as:

f i v i + m i g - k ∑ j=1 t j u j = 0 3 ( 5
) i m i = 0 3 ( 6
)
where k is the number of cables connected to the i th quadrotor and m i is the mass of the i th quadrotor.

AVAILABLE WRENCH SET

With the ultimate goal of determining the available wrench set of ACTSs, the mapping of cable tension space used with CD-PRs [START_REF] Ruiz | Arachnis: Analysis of robots actuated by cables with handy and neat interface software[END_REF][START_REF] Bouchard | On the ability of a cable-driven robot to generate a prescibed set of wrenches[END_REF] has been adapted to take into account quadrotor constraints. The main difference between the CDPR and the ACTS available wrench sets is that CDPRs have a constant available tension set regardless of the state of the robot, so for any configuration of the CDPR the minimum and maximum available cable tensions are constant. Because a component of the actuation of the quadrotor is required to keep itself in the air, different

x 0 z 0 u 1 u 2 u 3 θ 3 θ 2 φ 3 φ 2 φ 1 θ 1 FIGURE 3: ACTS WITH 3 QUADROTORS, 3

CABLES, AND

A POINT-MASS PLATFORM configurations will result in variable tension spans achievable by the ACTS, in turn affecting the set of available wrenches in task space. To this end, an extra mapping step is proposed to go from the thrust space of the quadrotors to the tension space of the cable system. During the explanation of the available wrench set calculations, an example of a planar ACTS (Figure 3) composed of three quadrotors with a single cable each attached to a point mass will be used. It has the configuration

[θ 1 , θ 2 , θ 3 ] = [0, π 4 , - 7π

Thrust Space

As the i th cable of this ACTS connects with the i th quadrotor at x i and can only transmit internal tensions, the thrust f i is the only component of i w i that gets transmitted to the payload. Quadrotor i's thrust f i is bounded by the minimum thrust f i required for hovering, and the maximum thrust fi that the quadrotor can generate. In this paper, an assumption is made that the tension force exerted by the cable on the quadrotor always has a positive projection along g (each quadrotor is above its attachment point(s)), allowing the minimum thrust to be determined by Eq. ( 7), as a function of the minimum allowable cable tension t, where t i > 0.

f i = k ∑ j=1 t j u j -m i g (7)
The upper bounds of the thrust space are simply the maximum thrust the quadrotors are able to generate. In the planar case study with three quadrotors, the thrust space is a 3-orthotope (Figure ). The thrust space H can be generalized for an arbitrary ACTS as the n-dimensional set of independently achievable quadrotor thrusts, forming an n-orthotope:

A) THRUST SPACE B) TENSION SPACE C) WRENCH SPACE H f 1 (N) f1 f 1 f2 ( N ) f2 f 2 f 3 (N) f3 f 3 T t 1 (N) t1 t 1 t2 ( N ) t2 t 2 t 3 (N) t3 t 3 W -m p g f x (N) f z (N) ¢ ¢ γ
H = f ∈ R n : f ≤ f ≤ f (8)
where

f = [ f 1 , . . . , f n ] T .

Tension Space Assumption 3: Cables tensions remain always positive

Tension space is the set of tensions that the quadrotors are able to exert on the cables. Each cable j has tension t j such that 0 < t j ≤ t j ≤ t j . The minimum tension is a chosen such that there will not be excessive sag in the cable, and the maximum tension is the maximum force the quadrotor can exert on the cable (assumed to be less than the cable's safe operating tension). Accordingly, the tension space can be generalized as:

T = {t ∈ R m : t ≤ t ≤ t} ( 9
)
where t = [t 1 , . . . ,t m ] T . As each cable is supported by a quadrotor of finite thrust and non-zero mass, a component of the thrust of the quadrotor must counteract gravity, with the remaining being available to support the cable tension. The balance of forces for the i th quadrotor in this example is given by Eq. 10.

f i v i + m i g -t i u i = 0 ( 10
)
The vector form can be decomposed into a system of three equa-tions and four variables, namely t i , f i , v ix , v iz .

f i v ix -t i u ix = 0 (11a) f i v iz + m i g -t i u iz = 0 (11b) v 2 ix + v 2 iz = 1 (11c)
with g being the g T [0, 1] T . Eliminating v ix and v iz reduces the problem to a single equation linking thrust and tension.

t 2 i -2u iz m i gt i + m 2 i g 2 -f 2 i = 0 ( 12
)
It is apparent that the maximum achievable tension must occur at the maximum quadrotor thrust, therefore the maximum tension can be mapped from thrust space to tension space by Eq. 13.

ti

= m i gu iz + f 2 i + m 2 i g 2 u 2 iz -1 (13)

Available Wrench Set

The wrench that each cable exerts on the platform can be mapped from tension space to wrench space with the wrench matrix W. The platform is in static equilibrium when Wt + w e = 0, where t ∈ T [START_REF] Azizian | The dimensional synthesis of spatial cable-driven parallel mechanisms[END_REF][START_REF] Ruiz | Arachnis: Analysis of robots actuated by cables with handy and neat interface software[END_REF].

W = u 1 . . . u m R p b 1 × u 1 . . . R p b m × u m (14)
The wrench space is the set of wrenches that the cables are able to exert on the payload expressed in F 0 as described in Eq. ( 15), where w i is the i th column of the wrench matrix W and ∆t i = tit i . A wrench w is inside the available wrench set if and only if it satisfies Eq. ( 15).

W = w ∈ R d |w = m ∑ j=1 α j ∆t j w j + Wt, 0 ≤ α j ≤ 1 ( 15
)
As the tension space in this example is an orthotope which is necessarily a convex set, the linear mapping of the points bounding the tension space is also a convex set. Therefore the mapping of the vertices of the tension space by the wrench matrix fully and exclusively contains all feasible wrenches. Using the convex hull method described in [START_REF] Bouchard | On the ability of a cable-driven robot to generate a prescibed set of wrenches[END_REF], the achievable wrench set of the system can be found as the convex hull of the set of d-dimensional points W a .

W a = Wdiag(∆t)α α α + Wt 1 (1,2 n ) ( 16
)
where 1 (1,2 n ) is a row vector of ones with 2 n columns. In the planar example considered here, the matrix α α α is the set of all combinations of tension limits bounding the tension space

α α α =   0 0 0 0 1 1 1 1 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 1   (17)
It is apparent that the mapping from tension space vertices to wrench space generates eight wrenches, however the wrench space in Fig. 4 has six vertices only. That is because, due to the redundant nature of the cable system (m > d), two mapped tension vertices are internal to the convex set formed by the other six wrenches.

Capacity Margin

The capacity margin γ is an index used to evaluate the robustness of static equilibrium [START_REF] Ruiz | Arachnis: Analysis of robots actuated by cables with handy and neat interface software[END_REF][START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF]. It is defined as the shortest distance from the expected load wrench to the nearest face of the available wrench set. Throughout the entirety of the wrench set determination, the orientation of the quadrotor thrust was held to be arbitrary, constrained solely by its maximum magnitude. This leads to the corollary that the capacity margin can be used to determine the minimum set of external wrenches on the platform that the ACTS can compensate for solely with modifications to the attitude of the quadrotor. As the dynamics of the attitude control of quadrotors are much faster than those of position control, this means that compensating for disturbances with attitude only will provide a more robust system. In this planar case, the attitude of the quadrotor is given by v i = [sin(φ ), cos(φ )] T , such that Eq. ( 10) is satisfied.

u 1 v 1 u 2 u 3 v 2 φ 1 φ 2

CASE STUDIES

Several case studies are provided to demonstrate the computation of the available wrench set. In the first case, a planar ACTS composed of a rigid platforms with two quadrotors and three cables is analyzed, differing from the previous example through a coupling between two dimensions of the tension space, as well as a non-homogenous wrench space. The study is then extended to two 6-DOF spatial systems, the first of which is a redundantly actuated platform with eight quadrotors and eight cables, the other being a system of three quadrotors and six cable such as is used in [START_REF] Manubens | Motion planning for 6-d manipulation with aerial towedcable systems[END_REF]. All examples use quadrotors with a mass of 0.8kg and a maximum thrust of 32N.

Planar System with 2 Quadrotors and 3 Cables

If the object being transported is large, the moments that can be applied to the platform become of greater importance. Multiple cables attached to a single quadrotor may improve the orientation control of the platform without increasing the number of quadrotors. For the ACTS shown in Fig. [START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF], the thrust space to tension space mapping of quadrotor 1 to cable 1 can be done using the same methodology as the previous example. The second quadrotor is connected to two cables, and has the equilibrium equation shown in Eq. ( 18)

f 2 v 2 + m 2 g -t 2 u 2 -t 3 u 3 = 0 (18)
To determine the tensions in the two cables, there are 5 variables ( f 2 ,t 2 ,t 3 , v 2x , v 2z ). An additional equation is found as v 2 is a unit

A) THRUST SPACE B) TENSION SPACE C) WRENCH SPACE H f 2 f 2 f 1 f1 f 2 (N) f 1 (N) T t 2 t 3 t 1 t2 t3 t1 t 2 (N) t 3 (N) t 1 (N) W f x (N) m y (Nm) f z (N) -m p g FIGURE 6:
THRUST, TENSION, AND WRENCH SPACES OF PLANAR ACTS WITH 2 QUADROTORS AND 3 CABLES vector, resulting in as system of 3 equations.

f 2 v 2x -t 2 u 2x -t 3 u 3x = 0 (19a) f 2 v 2z -t 2 u 2z -t 3 u 3z + m 2 g = 0 (19b) v 2 2x + v 2 2z = 1 (19c)
This set of equations can be reduced to an equation relating the thrust of the quadrotor to the magnitude of the two cable tensions.

t 2 2 + t 2 3 + 2t 2 t 3 u T 2 u 3 -2m 2 g(t 2 u 2z + t 3 u 3z ) + m 2 2 g 2 -f 2 2 = 0 (20)
By writing Eq. (20) to solve for t 2 as t 2 = h 2 (t 3 , f 2 ) and likewise t 3 = h 3 (t 2 , f 2 ) the cable tension space can be formulated as

T = t ∈ R 3      t 1 ≤ t 1 ≤ t1 t 2 ≤ t 2 ≤ h 2 (t 3 , f2 ) t 3 ≤ t 3 ≤ h 3 (t 2 , f2 ) (21)
The main difference in computing the available wrench set for this system is the α α α matrix used in Eq. ( 16). Appendix A proves that, given the sufficient conditions of all cables projecting along the positive z 0 axis, and the angle between cables j and k being less than π 2 rad, the function h j (t k , fi ) is convex over the domain t k ≤ t k ≤ tk . This property guarantees that the mapping of any set of points in the tension space produces a convex set exclusively containing feasible wrenches. In Fig. (6b), six vertices are clearly identifiable, however mapping only these points will result in an over conservative set of available wrenches. To approximate the full wrench set, the arc between t2 and t3 is divided into p + 1 segments. The combination matrix of the tension span of the second and third cables is denoted as α α α 23 .

α α α = 0 (1,3+p) 1 (1,3+p) α α α 23 α α α 23 ( 22
)
The α α α 23 matrix includes the three vertices of the t 2 ,t 3 plane, and a set of vertices evaluated at interpolation points

t 3,k = h 3 (t 2,k , f2 ) for all interpolation points k ∈ Z | 1 ≤ k ≤ p.
The full α α α 23 matrix is shown in Eq. ( 23).

α α α 23 =    0 1 0 k p+1 . . . p p+1 0 0 1 h 3 ( k∆t 2 p+1 + t 2 ) -t 3 ∆t 3 . . . h 3 ( p∆t 2 p+1 + t 2 ) -t 3 ∆t 3    (23)
As the available wrench set of this robot is w ∈ [ f x , f z , m y ], it must be homogenized to allow for a calculation of the capacity margin. In [START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF], the capacity margin is homogenized by dividing the moment by the radius of gyration r g = I/A of the platform, where A is the area of the platform and I is the second area moment of inertia. If W a is the matrix whose columns form a set of points bounding the convex wrench space, the homogenized wrench space boundary is given by Eq. (24).

W h =   I 2 0 0 1 r g   W a ( 24
)
where I 2 is an identity matrix of rank 2. In this case study the system parameters are defined as: pose p = [0, 0, pi The capacity margin of the ACTS in this configuration is 0.3N, so this configuration is on the border of equilibrium. It is noteworthy that solely reducing m p here does little to change γ, but that γ increases significantly with a combination of w e = [0, 0, -m y ] and decreased m p . The attitudes φ 1 and φ 2 of the quadrotors are chosen such that for the point in tension space t = W -1 (m p gw e ), Eq. ( 10) is satisfied for the first quadrotor, and Eq. ( 18) is satisfied for the second quadrotor.

Spatial System with 8 Quadrotors and 8 Cables

In this case, the platform is desired to have a controllable orientation without holonomic constraints on the twist. The ACTS, shown in Fig. [START_REF] Alothman | Using iterative lqr to control two quadrotors transporting a cablesuspended load[END_REF] has a cubic moving platform with side lengths of 0.5m, mass m p = 4kg, and orientation of [ π 12 , π 12 , 0]rad. It is manipulated by eight quadrotors so that all twists are achievable with two degrees of actuation redundancy in non-singular states. A configuration has been chose to avoid singularities for the evaluation of the available wrench set. The vector u i of the cables connected to the lower vertices of the platform are rotated by -π 6 rad from the direction of b i about z 0 , and the given an elevation of π 6 rad from the xy plane. The top vertices cables are rotated in the opposite direction by π 6 rad from b i about z 0 , and then elevated by π 3 rad. In this case, the thrust space is computed for each quadrotor as in the previous section, resulting in a 8-orthotope. The tension space is likewise calculated to result in an 8-orthotope by extending Eq. (11a) to a 3rd dimension, adding the equation f i v iyt i u iy = 0. The result is the same as in the planar system, with the thrust space to tension space mapping being described by Eq. ( 13). The vertices of the thrust and tension spaces are combinations of maximum and minimum thrust and tension values in Table [START_REF] Kumar | Opportunities and challenges with autonomous micro aerial vehicles[END_REF].

Mapping the vertices of the tension space orthotope to task space results in a convex 6-zonotope [START_REF] Gosselin | Cable-driven parallel mechanisms: state of the art and perspectives[END_REF][START_REF] Bouchard | On the ability of a cable-driven robot to generate a prescibed set of wrenches[END_REF] bounding the available wrench set of the ACTS. Fig. [START_REF] Manubens | Motion planning for 6-d manipulation with aerial towedcable systems[END_REF] shows the evaluation of the available wrench set at zero moment (8a) and at f = -m p g (8b). From this, it is seen that the system can exert a force wrench between 10N and 120N along z 0 , and forces between -20N and 20N laterally, with the largest lateral wrench range occuring when the vertical force is around 60N. This ACTS has little ability to withstand positive moments about z 0 limiting the capacity margin to γ = 6.5N. This ACTS configuration would be useful in situations that required (primarily vertical) movements, as well as roll and pitch capabilities, but without strict yaw requirements. While the wrench feasibility can be quantified using the capacity margin after normalizing W a by Eq. ( 25), the configuration of this particular ACTS does not lend itself well toward numerical optimization, as the set of variables (the altitude and azimuth of each cable) is 16 dimensional. Intuitive configurations avoiding common singularities can however, be evaluated for their robustness.

W h =   I 3 0 0 I 3 r g   W a (25)
The attitude of the quadrotors must be controlled so as to chose v i satisfying f i v i + m i gt i u i = 0. As the system is redundant, tension distribution algorithms may be used to specify the desired tensions, and therefor attitudes. In this case, using the Moore-Pensrose pseudo-inverse of W finds the least norm cable tensions, and vector λ λ λ is chosen so that t is within T . The vector v i may be found for all quadrotors using Eq. [START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF].

t = W † (w e -m p g) + Nλ λ λ , where N = null(W) (26)
The rotation matrix of each quadrotor is under-constrained by 1-DOF, as the rotation of the i th quadrotor around v i has no affect on the system equilibrium [START_REF] Lee | Geometric control of cooperating multiple quadrotor uavs with a suspended payload[END_REF], and any rotation matrix of the under-constrained quadrotor may be chosen so long as the equation R i [0, 0, 1] T = v i is satisfied.

Spatial System with 3 Quadrotors and 6 Cables

The Flycrane ACTS has been studied for motion planning in [START_REF] Manubens | Motion planning for 6-d manipulation with aerial towedcable systems[END_REF], and is composed of three quadrotors and six cables. A benefit of this design is the ability to generate six independent wrenches using only three quadrotors. It also has the advantage of being significantly easier to optimize, with a three dimensional configuration space compared to the 16 dimensional configuration space in the previous case.

The ACTS in this case is displayed in Fig. [START_REF] Bosscher | Wrench-feasible workspace generation for cable-driven robots[END_REF]. platform orientation is [ π 12 , π 12 , 0]rad and m p = 4kg. The platform is parameterized as b i = 0.5 cos 2πi 6 , sin 2πi 6 , 0 m, i ∈ Z, 1 ≤ i ≤ 6. The cable pairs have lengths of L 1,2 = 0.35m, L 3,4 = 0.70m, and L 5,6 = 1.40m. Each quadrotor is at an inclination of π 4 rad from the z p axis.

The 3-orthotope thrust space is mapped to the 6-dimensional cable tension space, with the tension in each pair of cables ( j and k) being related to the thrust of a single quadrotor (i) through Eq. (27).

f i v i + m i g -t j u j -t k u k = 0 (27)
To determine the tensions in the two cables, there are 6 variables (t j ,t k , f i , v ix , v iy , v iz ). An additional equation is found as v i is a unit vector, resulting in 4 equations with 6 unknowns. This reduces to the same result as seen in Eq. (20), generalized in Eq. (28) for cables j and k.

t 2 j + t 2 k + 2t j t k u T j u k -2m i g(t j u jz + t k u kz ) + m 2 i g 2 -f 2 i = 0 (28)
Each of the three orthogonal 2-dimensional tension space components, which are all composed of straight lines along the lower bound tensions, and the upper bound tension is described by an analytic function h, where t j = h j (t k , fi ), creating a 6 dimensional solid composed of a combination of flat and curved It is seen that for longer cables (u T j u k ≈ 1), the upper bound of the curve approaches a straight line joining the vertices t j , h k (t j , fi ) and h j (t k , fi ), t k . As the number of interpolation points p along the curve t j = h j (t k , fi ) increases, the size of the convex set of points increases by the cube of the number of points (3 + p) 3 . Short cables relative to the baseline separating the cables increases the curvature of the maximum tension curve, increasing the importance of interpolations. The vertices of the thrust space and the cable tension space are shown in Table (2). In the same manner as is done for the planar system with two quadrotors and three cables, the available wrench space of this ACTS can be mapped through interpolation along the convex curved surfaces. The combination matrix α α α is built as all combinations of the columns of α α α jk , derived as in Eq. ( 23), and arranged as shown in Eq. ( 29) The available wrench set of the ACTS is shown in Fig. [START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF]. The ACTS in this configuration can apply a peak force of 50N along z 0 , and a peak force of 30N in the positive x 0 direction, however it has can only exert a force of 2N along the negative x 0 axis. It is able to exert moments from -10 Nm to 15 Nm about z 0 , but much less around x 0 and y 0 .

α α α = comb     α α α 12 α α α 34 α α α 56     (29)

CONCLUSIONS AND FUTURE WORK

In this paper, it is shown how the available wrench set of the payload of an aerial cable towed system (ACTS) can be calculated. This method is adaptable to systems where two cables are coupled due to sharing the actuation of a single quadrotor. It is proven that the relationship between tensions in such a coupled system is necessarily convex given certain operational and design conditions, showing that the available wrench set obtained from a mapping of the tension space is guaranteed to be conservative throughout the wrench space.

From a practical point of view, this has applications in the optimal design and motion planning of ACTS with the objective of increasing the ability to resist external wrenches acting on the moving platform. As the available wrench set is determined through a mapping of the full achievable tension space of the system, it can be seen that the position of the platform and the configuration of the system can be maintained while the external wrench is within the bounds of the available wrench set through changes in attitude of the quadrotor, requiring no change in the quadrotors positions. The higher dynamic response of attitude control compared to position control may lead to improvements in response time to system disturbances.

Future work will include experimental validation of the available wrench set calculations by testing the reaction of ACTSs to wrenches lying inside and outside of the available wrench space. The analysis will then be extended to dynamic systems, and the implementation of the available wrench set calculation for generating trajectories optimized for the robustness of equilibrium of the system will be performed.

NOMENCLATURE

FIGURE 1 :

 1 FIGURE 1: PARAMETERIZATION OF AN AERIAL CABLE TOWED SYSTEM COMPOSED OF THREE QUADROTORS, THREE CABLES, AND A RIGID PLATFORM

FIGURE 2 :

 2 FIGURE 2: FREE BODY DIAGRAMS OF THE MOVING PLATFORM AND THE i th QUADROTOR

  18 ]rad, and a payload mass of m p = 4kg. The quadrotors used have a mass of m i = 0.8kg and maximum thrust of fi = 32N. In a general case, there are n quadrotors, m cables, and a task space of dimension d, satisfying the following requirements: 1 ≤ d ≤ 6 and d ≤ m and n ≤ m

FIGURE 4 :

 4 FIGURE 4: THRUST, TENSION, AND WRENCH SPACES OF PLANAR ACTS WITH 3 QUADROTORS AND 3 CABLES

FIGURE 5 :

 5 FIGURE 5: DIAGRAM OF PLANAR SYSTEM WITH 2 QUADROTORS AND 3 CABLES

/ 6 ]

 6 , platform mass m p = 4kg, platform dimensions b 1 = [-0.75, -0.25] T , b 2 = [0.5, 0.5] T , b 3 = [1.0, -0.25] T , cables of 0.75m in length, and u 1 = [0, 0, 1] T . Fig. (6c) shows the available wrench set, with the dotted line indicating the line along f z passing through f x = 0, m y = 0.

FIGURE 7 :

 7 FIGURE 7: ACTS WITH 8 QUADROTORS AND 8 CABLES

FIGURE 8 :

 8 FIGURE 8: ACHIEVABLE WRENCH SPACE OF 8 QUADRO-TOR, 8 CABLE ACTS EVALUATED AT MOMENT = [0, 0, 0]Nm (LEFT) AND FORCE = -m p g N (RIGHT)

FIGURE 9 :

 9 FIGURE 9: OF FLYCRANE ACTS

FFIGURE 10 :

 10 FIGURE 10: ACHIEVABLE WRENCH SPACE OF FLY-CRANE ACTS EVALUATED AT MOMENT = [0, 0, 0]Nm (LEFT) AND FORCE = -m p gN (RIGHT)

TABLE 1 :

 1 THRUST SPACE AND CABLE TENSION SPACE LIMITS FOR 8 QUADROTOR, 8 CABLE CASE STUDY.

	Quadrotor f(N) f (N) Cable t(N) t(N)
	1	8.7	32.0	1	1.0	25.0
	2	8.4	32.0	2	1.0	27.4
	3	8.7	32.0	3	1.0	25.0
	4	8.4	32.0	4	1.0	27.4
	5	8.7	32.0	5	1.0	25.0
	6	8.4	32.0	6	1.0	27.4
	7	8.7	32.0	7	1.0	25.0
	8	8.4	32.0	8	1.0	27.4

TABLE 3 :

 3 LIST OF NOMENCLATURE, DIVIDED INTO GENERAL, QUADROTOR, AND CDPR CLASSES Number of quadrotors x 0 , y 0 , z 0 x, y, and z unit vectors of F 0 g Gravitational vector [0, 0, -9.81] T ms -2 g Scalar gravity g T z 0 F i Reference frame of quadrotor i A i Attachment point of cables to quadrotor i (in F 0 ) C F ,C M Coefficients of thrust and drag of propellors f i , fi Min. and max. allowable thrust of quadrotor i m i Mass of quadrotor i R i Rotation matrix of F i with respect to F 0 v i Unit vector normal to propellors of quadrotor i x i Position of quadrotor i expressed in F 0 ω i Angular rotation of propellor i Vector from p to cable i connection point (in F p ) B i Point where cable i connect to platform (in F 0 ) l i Length of cable i m p Mass of platform p Position of platform center of mass expressed in F 0 0 R p Rotation matrix of F p with respect to F 0 t i Force vector exerted by cable i at point B i t Vector of cable tensions [t 1 , . . . ,t m ] T t i , ti Min. and max. allowable tensions in cable i u i Unit vector from B i to A i W Wrench matrix of platform W a Set of available wrench vectors of the ACTS w e External wrench acting on platform

	Symbol Description
	F 0	Inertial reference frame
	H	Thrust space of ACTS
	T	Tension space of ACTS
	W	Wrench space of ACTS
	d	Dimension of task space
	m	Number of cables
	n
	F p	Reference frame of platform
	α α α	Combination matrix
	b i

Appendix A: Proof of Convex Tension Space

This appendix proves that the function relating the maximum tensions of two cables i, j connected to quadrotor Q forms a convex tension space. Equation 20 can be rewritten as

where

The maximum tension of cable j, t j relates to ti through the function t j = h j (t i , f Q ), with t j = h j t i , fQ . The function is convex if, over the span of allowable tensions, it can be shown that dt j dt i ≤ 0 and that

Using implicit differentiation, it can be found that the derivative of t j with respect to t i is

Assumption number 4 requires that u z be positive for all cables, therefore k 2 and k 3 are both positive. Applying a simple design rule l 2 i + l 2 j ≥ B j -B i forces 0 ≤ k 1 ≤ 1. With strictly positive tensions, and k 1 , k 2 , k 3 ≥ 0, dt j dt i is strictly negative. By setting dt j dt i = -g h , the quotient rule and the law of implicit differentiation can be used to solve the second derivative.

It is therefore shown that within the domain of t i , given that all cable vectors have a positive projection along the z axis, and that the angle between two cables connected a quadrotor is acute, the tension space is a convex shape.


Coralie Germain
email: germain@irccyn.ec-nantes.fr

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

Stéphane Caro
email: caro@irccyn.ec-nantes.fr

Philippe Wenger
email: wenger@irccyn.ec-nantes.fr

Natural Frequency Computation of Parallel Robots

come   L'archive ouverte pluridisciplinaire

Introduction

Parallel robots have been increasingly used in industry in the last few years for pick-and-place applications and high-speed machining [START_REF] Brogardh | Present and future robot control development an industrial perspective[END_REF][START_REF] Chanal | A study of the impact of machine tool structure on machining processes[END_REF]. This interest is mainly due to their high stiffness and good dynamic performance compared with their serial manipulator counterparts.

Having a good knowledge of the elastodynamic behaviour of a manipulator (especially its natural frequencies) is a crucial point. As a consequence, accurate elastodynamic models of parallel robots are useful at both their design [START_REF] Voglewede | Overatching framework for measuring the closeness to singularities of parallel manipulators[END_REF][START_REF] Bouzgarrou | New approach for dynamic modelling of flexible manipulators[END_REF][START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF] and control [START_REF] Singer | Preshaping Command Inputs to Reduce System Vibration[END_REF][START_REF] Singh | Tutorial on input shaping/time delay control of maneuvering flexible structures[END_REF][START_REF] Pelaez | Input shaping reference commands for trajectory following cartesian machines[END_REF] stages in order to define their optimal dimensions and shapes while improving their vibratory behaviour.

Several models have been proposed and used in order to compute the natural frequencies of a mechanism. For instance, the following three methods are commonly used in the literature:

• Finite Element Analysis (FEA); the FEA method is proved to be the most accurate and reliable, since the links and joints are modeled with their true dimensions and shape [START_REF] Bouzgarrou | Rigidity analysis of t3r1 parallel robot uncoupled kinematics[END_REF][START_REF] Bettaieb | Modeling of a high-speed machining center with a multibody approach: the dynamic modeling of flexible manipulators[END_REF]. However, its accuracy is limited by the mesh size and is usually used at the final design stage of the robot because the meshing is time-consuming.

• Matrix Structural Analysis (MSA) method is a common technique in mechanical engineering [START_REF] Bouzgarrou | New approach for dynamic modelling of flexible manipulators[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Rognant | A systematic procedure for the elastodynamic modeling and identification of robot manipulators[END_REF][START_REF] Bauchau | Flexible multibody dynamics[END_REF][START_REF] Cammarata | An algorithm to study the elastodynamics of parallel kinematic machines with lower kinematic pairs[END_REF]. It incorporates the main ideas of the FEA method but operates with rather large flexible elements (beams, arcs, cables, and so on). As a result, the MSA method is less time consuming than the FEA method, but requires good knowledge in FEA.

• Virtual joint method (VJM ) [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF][START_REF] Wittbrodt | Dynamics of Flexible Multibody Systems[END_REF], which is also referred to as "lumped modeling", is based on the expansion of the traditional rigid model by adding virtual joints (localized springs), which describe the elastic deformations of the manipulator components (links, joints and actuators). Lumped modeling is simpler to use but is less accurate than MSA.

Some general methodologies dealing with MSA or VJM are presented in [START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Rognant | A systematic procedure for the elastodynamic modeling and identification of robot manipulators[END_REF][START_REF] Bauchau | Flexible multibody dynamics[END_REF][START_REF] Wittbrodt | Dynamics of Flexible Multibody Systems[END_REF]. These methodologies can differ from the type of elements or flexible models used. Nevertheless, all of them require the determination of some Jacobian matrices in order to characterize kinematic dependencies in closed-loop mechanisms. The main drawback of such methodologies is that they are not specifically designed for parallel robots and they do not propose a systematic and straightforward way for computing the Jacobian matrices required for closing the loops.

An approach for the systematic computation of these matrices was proposed in [START_REF] Bouzgarrou | New approach for dynamic modelling of flexible manipulators[END_REF]. However, this approach is complicated and the obtained results have not been verified either with a commercial software or experimentally. In [START_REF] Cammarata | An algorithm to study the elastodynamics of parallel kinematic machines with lower kinematic pairs[END_REF], the authors proposed an interesting algorithm to study the elastodynamic behaviour of parallel robots that considers both the joint and link flexibilities. The way the authors close the robot loops and choose the independant coordinates is not straightforward, resulting in a potentially asymmetrical description of the leg variables.

This paper aims at developing a simple and straightforward procedure for the computation of mass and stiffness matrices of parallel robots in any stationary configuration. Indeed,

• a simple way to compute the Jacobian matrices required for closing the robot loops is presented, and

• the set of independent coordinates for parallel robots with a symmetrical arrangement are chosen cleverly.

Moreover, contrary to most of the existing methods, the proposed approach does not contain any numerical matrix inversion, which is better to avoid numerical issues that may lead to a loss in the result accuracy. The natural frequencies calculated using this procedure will be compared with simulation data obtained using a FEA/MSA software and with experiments carried out on a parallel robot. A methodology similar to the one presented in [START_REF] Ibrahim | Inverse and direct dynamic models of hybrid robots[END_REF] to obtain the rigid dynamic model of a parallel robot is used. For a parallel robot composed of rigid or flexible links connected by passive or active joints, such as the one described in Fig. 1, the methodology is decomposed into two steps:

1. All closed-loops are virtually opened in such a way that the platform is virtually disassembled from the robot architecture (Fig. 1(b)). Each joint is virtually considered actuated so that the robot becomes a tree structure with a free body: the platform. The mass and stiffness matrices of the tree structure and the ones of the free platform are then computed;

2. The loops are closed using Jacobian matrices so that the mass and stiffness matrices of the actual parallel robot can be obtained.

The paper is organized as follows. In Section 2, a method used for the computation of the mass and stiffness matrices of one single flexible body using MSA techniques is recalled. The mass and stiffness matrices of the virtual system composed of the tree structure and the free platform are introduced in Section 3. A straightforward way for the computation of the Jacobian matrices linking the Cartesian coordinates of each body to the generalized coordinates of the virtual system is proposed. Section 4 describes the proposed method to close the loops in order to obtain the mass and stiffness matrices of the actual parallel robot, that keeps a symmetrical description of the leg variables. The natural frequencies of the NaVARo [START_REF] Rakotomanga | Kinetostatic performance of a planar parallel mechanism with variable actuation[END_REF][START_REF] Caro | Kinematic and dynamic modeling of a parallel manipulator with eight actuation modes[END_REF] are then evaluated in Section 5. NaVARo is a three-degree-of-freedom (dof) planar parallel manipulator with multiple actuation modes developed at IRCCyN. Moreover, the simulation results obtained with our approach are compared with both Cast3M software [START_REF] Castem | The castem software webpage[END_REF] and experimentations. Finally, some conclusions and future works are proposed in Section 6.

Fixed base

Moving platform passive joints actuated joints

links C 1,1 C 1,2 C 1,3 C 1,m 1 C 1,m 1 -1 C n,1 C n,2 C n,3 C n,m n C n,m n -1 C k,1 C k,2 C k,3 C k,m k C k,m k -1 C 2,1 C 2,2 C 2,3 C 2,m 2 C 2,m 2 -1 (a) Kinematic chain
M j (M 0j ) M 0j B 0j B j (B 0j ) A 0j =A j F b0 F j v j (A j ) ω j (A j ) u e j =Σ φ dk j (M 0j )q ek j

Kinematics of a Flexible Free Body

Let us consider body j described in Fig. 2(a) to which a local frame F j is attached at point A j . The flexible body is considered as the superposition of the undeformed configuration referred to as 0, which corresponds to a rigid transformation of body j described by variable q j and the deformed configuration related to the elastic deformation of the body. The translational velocity of any material point M j of this flexible body can be expressed as [START_REF] Boyer | Dynamics of a 3-dof spatial parallel manipulator with flexible links[END_REF]:

v j (M j ) = v j (A j ) + ω j (A j ) × r j (M j ) + v e j (M j ) (1)
where v j (A j ) and ω j (A j ) are the translational and rotational velocities of body j expressed in frame F j , respectively, r j (M j ) is the position vector of point M j (of the deformed body) expressed in frame F j , v e j (M j ) is the translational velocity due to the body elasticity that can be parameterized as truncated series of Rayleigh-Ritz shape functions:

v e j (M j ) = Φ d j (M 0j ) qe j (2)
with

Φ d j = φ d1 j • • • φ dN j j of size (3 × N j ), φ dk j (M 0j
) being a three-dimensional vector containing the k-th shape functions for the displacement of the flexible body at point M 0j , and

qe j = qe1 j • • • qeN j j T
, qek j being the k-th elastic generalized velocity of body j and N j the number of considered shape functions. It should be noted that vector r j (M j ) in (1) can be expressed as:

r j (M j ) = r j (M 0j ) + u e j (M 0j ) (3)
where r j (M 0 j ) is the position vector of point M 0j expressed in frame F j ,

u e j (M 0j ) = Φ d j (M 0j )q e j (4)
is the elastic displacement that transforms M 0j to M j and q e j = q e1 j • • • q eN j j T is the vector of elastic generalized coordinates of body j.

Equation (1) can thus be written the following matrix form:

v j (M j ) = I 3 [r j (M j )] T × Φ d j (M 0j )   v j (A j ) ω j (A j ) qe j   (5)
where I 3 is the (3 × 3) identity matrix and [r j (M j )] × is the cross-product matrix of vector r j (M j ).

Equations ( 1) to ( 5) define the kinematic model of flexible free body j. This model is thus parameterized by the following set of variables:

• v j (A j ) and ω j (A j ): Cartesian velocities (Euler variables) characterizing the rigid displacement of body j expressed in frame F j ;

• q e j : generalized coordinates (Lagrange variables) characterizing the elastic displacements of body j.

It should be mentioned that this description can be applied to both robot links and joints, as long as all the shape functions can be defined.

Kinetic Energy of a Flexible Free Body

The kinetic energy of body j (denoted as Σ j in the following integral) is given by:

T j = 1 2 Σ j v j (M j ) T v j (M j )dm = ρ j 2 V j v j (M j ) T v j (M j )dV (6)
where ρ j is the material density and V j is the body volume. Introducing (5) into (6) leads to:

T j = ρ j 2 V j   v j (A j ) ω j (A j ) qe j   T   I 3 [r j (M j )] × Φ d j (M 0j ) T   I 3 [r j (M j )] T × Φ d j (M 0j )   v j (A j ) ω j (A j ) qe j   dV = 1 2   v j (A j ) ω j (A j ) qe j   T M j   v j (A j ) ω j (A j ) qe j   (7)
where

M j = ρ j V j    I 3 [r j (M j )] T × Φ d j (M 0j ) [r j (M j )] × [r j (M j )] × [r j (M j )] T × [r j (M j )] × Φ d j (M 0j ) Φ d j (M 0j ) T Φ d j (M 0j ) T [r j (M j )] T × Φ d j (M 0j ) T Φ d j (M 0j )    dV (8)
M j is the mass matrix of body j expressed in its frame F j .

Elastic Potential Energy and Stiffness Matrix of a Flexible Free Body

The elastic potential energy of any body is given by [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]:

V e j = 1 2 V j σ T j I t ǫ j dV (9)
where σ j and ǫ j are the six-dimensional stress and strain vectors due to the small elastic displacement u e j (M 0j ) in body j. I t is a (6 × 6) diagonal matrix. The first three diagonal terms are equal to one, whereas the last three diagonal terms are equal to two, because of the two multipliers associated with the shear strains [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]. The strain vector is defined as ǫ j = ǫ j 11 ǫ j 22 ǫ j 33 ǫ j 12 ǫ j 13 ǫ j 23 T , where:

  ǫ j 11 ǫ j 12 ǫ j 13 ǫ j 12 ǫ j 22 ǫ j 23 ǫ j 13 ǫ j 23 ǫ j 33   = 1 2 ∇ u e j (M 0j ) + ∇ u e j (M 0j ) T (10)
with

∇ u e j (M 0j ) = ∂u e j ∂x (M 0j ) ∂u e j ∂y (M 0j ) ∂u e j ∂z (M 0j ) =          ∂Φ 1 d j ∂x (M 0j )q e j ∂Φ 1 d j ∂y (M 0j )q e j ∂Φ 1 d j ∂z (M 0j )q e j ∂Φ 2 d j ∂x (M 0j )q e j ∂Φ 2 d j ∂y (M 0j )q e j ∂Φ 2 d j ∂z (M 0j )q e j ∂Φ 3 d j ∂x (M 0j )q e j ∂Φ 3
d j ∂y (M 0j )q e j ∂Φ 3
d j ∂z (M 0j )q e j          (11)
where Φ k d j corresponds to the k-th line of matrix Φ d j , k = 1, 2, 3. As a result,

ǫ j =                        ∂Φ 1 d j ∂x (M 0j ) ∂Φ 2 d j ∂y (M 0j ) ∂Φ 3 d j ∂z (M 0j ) 1 2 ∂Φ 1 d j ∂y (M 0j ) + ∂Φ 2 d j ∂x (M 0j ) 1 2 ∂Φ 1 d j ∂z (M 0j ) + ∂Φ 3 d j ∂x (M 0j ) 1 2 ∂Φ 2 d j ∂z (M 0j ) + ∂Φ 3 d j ∂y (M 0j )                        q e j = Φ ǫ j q e j (12)
The stress vector is expressed as:

σ j = σ j 11 σ j 22 σ j 33 σ j 12 σ j 13 σ j 23 T = E j ǫ j (13)
where matrix E j is given by the Hooke's law [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]. Thus, introducing ( 12) and ( 13) into (9) leads to:

V e j = 1 2 q T e j K j q e j ( 14
)
where K j is the stiffness matrix of body j and takes the form:

K j = V j Φ T ǫ j E T j I t Φ ǫ j dV (15)

Case of the Elastic Beam: the Bernoulli Model

The computation of the mass and stiffness matrices of 3D beams is useful for the elastodynamic modeling of parallel manipulators. The Bernoulli model describes the beam deformation under the assumption that the shear effect is negligible, that the cross-sections remain perpendicular to the neutral axis and that the rotational inertia of sections is assumed to be zero [START_REF] Blevins | Formulas for natural frequency and mode shape[END_REF]. With such a model, the 3D beam deformation u e j (M 0j ) (see Section 2.1) can be characterized with the six shape functions φ dx j , φ dy j , φ dz j , φ rx j , φ ry j and φ rz j , i.e. N j = 6, defined as:

φ dx j = ξ 0 0 0 0 0 (16a) φ dy j = 0 3ξ 2 -2ξ 3 0 0 0 l j (ξ 3 -ξ 2 ) (16b) φ dz j = 0 0 3ξ 2 -2ξ 3 0 -l j (ξ 3 -ξ 2 ) 0 (16c) φ rx j = 0 0 0 ξ 0 0 (16d) φ ry j = 0 0 -6 (ξ -ξ 2 ) /l j 0 3ξ 2 -2ξ 0 (16e) φ rz j = 0 6 (ξ -ξ 2 ) /l j 0 0 0 3ξ 2 -2ξ (16f)
where ξ = x/l j and l j is the beam length.

x, y and z denote the Cartesian coordinates of point M 0j expressed in the local frame F j and Φ d j (M 0j ) is a (3 × 6) matrix that takes the form:

Φ d j (M 0j ) =   φ dx j -yφ rz j + zφ ry j φ dy j -zφ rx j φ dz j + yφ rx j   (17)
Moreover, in the beam model, it is assumed that

σ j 22 = σ j 33 = σ j 23 = 0 (18) ǫ j 22 = ǫ j 33 = ǫ j 23 = 0 (19) σ j 11 = E j ǫ j 11 (20) σ j 12 = G j ǫ j 12 (21) σ j 13 = G j ǫ j 13 ( 22
)
where E j is the Young modulus of body j and G j is its shear modulus.

Introducing (16a) to ( 22) into ( 8) and [START_REF] Wittbrodt | Dynamics of Flexible Multibody Systems[END_REF] for q e j = 0, the stiffness matrix of body j takes the form:

K j = 1 l 3 j         E j A j l 2 j 0 0 0 0 0 0 12E j I z j 0 0 0 -6E j I z j l j 0 0 12E j I y j 0 6E j I y j l j 0 0 0 0 I 0 j G j l 2 j 0 0 0 0 6E j I y j l j 0 4E j I y j l 2 j 0 0 -6E j I z j l j 0 0 0 4E j I z j l 2 j         (23)
where A j is the beam cross-section area, I y j and I z j are the second moments of area around axes y and z of the local frame, I 0 j is the torsion constant.

Similarly, the mass matrix of body j is expressed as:

M j = M 11 j M 12 j M T 12 j M 22 j (24)
where

M 11 j =         m j 0 0 0 m j z G j -m j y G j 0 m j 0 -m j z G j 0 m j x G j 0 0 m j m j y G j -m j x G j 0 0 -m j z G j m j y G j J xx j J xy j J xz j m j z G j 0 -m j x G j J xy j J yy j J yz j -m j y G j m j x G j 0 J xz j J yz j J zz j         (25)
M 12 j =                  m j 2 0 0 0 0 0 0 m j 2 0 0 0 m j l j 12 0 0 m j 2 0 m j l j 12 0 0 0 0 ρ j l j I p j 2 0 0 0 0 -ρ j I y j - 7m j l j 20 0 - m j l 2 j 20 0 0 ρ j I z j + 7m j l j 20 0 0 0 - m j l 2 j 20                  (26)
M 22 j =            m j 3 0 0 0 0 0 0 13m j 35 + 6ρ j Iz j 5l j 0 0 0 - 11m j l j +21ρ j Iz j 210 0 0 13m j 35 + 6ρ j Iy j 5l j 0 11m j l j +21ρ j Iy j 210 0 0 0 0 ρ j l j Ip j 3 0 0 0 0 11m j l j +21ρ j Iy j 210 0 m j l 2 j +14ρ j Iy j l j 105 0 0 - 11m j l j +21ρ j Iz j 210 0 0 0 m j l 2 j +14ρ j Iz j l j 105            (27) and
• x G j , y G j and z G j are the Cartesian coordinates of the center of mass of body j expressed in frame F j ;

• J xx j , J yy j , J zz j , J xy j , J xz j and J yz j are the terms of the rigid inertia matrix expressed at the origin of frame F j ;

• I p j = I y j + I z j is the polar moment of inertia.

Computation of the Stiffness and Mass Matrices of a Tree Structure Robot

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 elements1 connected by m i joints (revolute, prismatic or fixed joints -i = 1, ..., n) located at points C ik (k = 1, ..., m i -Fig. 1(a)). The j-th element of the i-th leg is denoted by ij. The nominal values of the actuated variables (of the passive variables, resp.) are denoted by q a (q p , resp.). The nominal values of the Cartesian coordinates of the platform are denoted by x p . The dimension n a of vector q a must be greater than or equal to the number of degrees of freedom of the parallel robot. Under the assumption of an elastic deformation, the variations in those variables are denoted by δq a , δq p and δx p , respectively. The number of shape functions for the element ij is denoted by N ij (j = 1, ..., m i -1). As a result, the dimension n e of the vector of elastic variables q e is equal to n i=1 m i -1 j=1 N ij . Thus, the vector of generalized coordinates of the tree-structure is given by q

t = q T t 1 • • • q T tn T ,
where

q t i = δq T a i δq T p i q T e i
T . δq a i , δq p i and q e i are the vectors of the actuated, passive and elastic generalized coordinates for the i-th leg.

Relationships between the Generalized Coordinates of a Flexible Body and the Generalized Coordinates of the Tree Structure

The generalized velocities v ij (A ij ), ω ij (A ij ) and qe ij of body ij are linked to the generalized velocities2 qt by the relations:

v ij (A ij ) ω ij (A ij ) = J i v ij qt i = J v ij qt (28)
with

J v ij = 0 • • • J i v ij • • • 0 (29a) qt =        qt 1 . . . qt i . . . qtn        (29b)
and

qe ij = J e ij qt (30)
The foregoing equations can be rewritten in the following compact form:

  v ij (A ij ) ω ij (A ij ) qe ij   = J ij qt (31)
with

J ij = J v ij J e ij ( 32
)
J ij is the Jacobian matrix mapping the generalized velocities qt of the tree structure into the generalized velocities of the element ij. J v ij is the Jacobian matrix mapping qt into the twist of the local frame attached to element ij. J e ij is a matrix composed of 0 and 1 terms that is used to extract vector qe ij from vector qt . Matrix J v ij is expressed thanks to a generic and systematic approach described thereafter and derived from [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF].

Let us first compute the homogeneous transformation matrix T ij that defines the location and orientation of the local frame attached to the element ij. By definition, this element is located in the leg i and is preceded by j -1 elements, each element ik (k = 1, ..., j) being linked to another by a joint (revolute, prismatic or fixed) described by its nominal coordinate q ik (Fig. 1(b)). Thus, matrix T ij is defined by:

T ij = T i base j-1 k=1 V a (q ik )T ik elt V e (q e ik ) V a (q ij )T ij end ( 33
)
where

• T i base denotes the rigid transformation between the global base frame and the frame attached to i-th leg;

• The matrix function V a (.) is a transformation matrix corresponding to an elementary rotation or an elementary translation;

• T ik elt denotes the rigid transformation matrix between the frame attached to element ik and the frame attached to element i, k + 1 in the undeformed case;

• The matrix function V e (.) represents the translations and rotations due to the deformations of the flexible link (if the element is rigid, this matrix will be the identity matrix q e ik = 0 in this case);

• T ij end is a matrix that allows the rotational part of T ij to be equal to the identity matrix. Then, let us gather all variables introduced in (33) into vector

q t ij = q i1 q T e i1 • • • q T e i,j-1 q ij T ( 34
)
and assume that for the l-th variable q t ijl of vector q t ij , (33) is rewritten as:

T ij = H R ijl V ijl (q t ijl )H L ijl (35)
where the first and the third multipliers are constant homogeneous matrices while the second multiplier is either an elementary translation matrix or an elementary rotation matrix. Then the partial derivative of the homogeneous matrix T ij with respect to q t ijl at the configuration q nom t ijl (q nom t ijl = 0 for an elastic variable and may not vanish for a joint variable) can be computed from a similar product where the internal term is replaced by the matrix V d ijl (.) = ∂V ijl (q t ijl )/∂q t ijl that takes a simple analytical form. For instance, for elementary translations and rotations about the x axis, these derivatives take the form:

V d ijl (q t ijl ) =     0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0     for a translation along x axis V d ijl (q t ijl ) =     0 0 0 0 0 -sin(q t ijl ) -cos(q t ijl ) 0 0 cos(q t ijl ) -sin(q t ijl ) 0 0 0 0 0     for a rotation about x axis (36)
For the elementary translations and rotations about the other axes, the matrices V d ijl (.) can be obtained similarly.

Only small displacements occur in the determination of the natural frequencies. Therefore, the derivative of the homogeneous matrix

T d ij = ∂T ij /∂q t ijl = H R ijl V d ijl (q t ijl )H L
ijl can be expressed as:

T d ij =     0 -γ d z ijl γ d y ijl λ d x ijl γ d z ijl 0 -γ d x ijl λ d y ijl -γ d y ijl γ d x ijl 0 λ d z ijl 0 0 0 0     (37)
It is noteworthy that vector 0 j

qt ijl v ij = λ d x ijl λ d y ijl λ d z ijl γ d x ijl γ d y ijl γ d z ijl T
, which can be obtained by extracting the terms of matrix T d ij defined in (37) is the column of the Jacobian matrix, corresponding to variable qt ijl . The latter maps vector q t into the velocity of the local frame attached to body ij expressed in the global frame [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. Let ij R 0 be the (6 × 6) extended rotation matrix:

ij R 0 = ij R 0 0 3 0 3 ij R 0 ( 38
)
where ij R 0 is the rotation matrix between the global frame and the local frame attached to element ij, which is evaluated in the robot undeformed configuration and can be extracted using (33). Thus, multiplying ij R 0 by vector 0 j qt ijl v ij yields the column j

qt ijl v ij = ij R 0 0 j qt ijl
v ij of Jacobian matrix J v ij defined in (28) corresponding to variable qt ijl .

Stiffness and Mass Matrices of the Tree Structure

From ( 6) and ( 9), the Lagrangian of the tree structure system can be expressed as:

L t = i,j T ij -V e ij = 1 2 i,j      v ij (A ij ) ω ij (A ij ) qe ij   T M ij   v ij (A ij ) ω ij (A ij ) qe ij   -q T e ij K ij q e ij    (39)
Introducing ( 28) into (39) leads to

L t = 1 2 i,j qT t J T ij M ij J ij qt -q T t J T e ij K ij J e ij q t = 1 2 qT t M t qt -q T t K t q t ( 40
)
where M t and K t are the mass and stiffness matrices of the tree structure.

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

L = 1 2 qT t M t qt + ẋT p M p ẋp -q T t K t q t = 1 2 qT t ẋT p M t 0 0 M p qt ẋp -q T t δx T p K t 0 0 0 q t δx p = 1 2 qT tot M tot qtot -q T tot K tot q tot (41)
where M p is the mass matrix of the rigid platform. M tot and K tot are the total mass and stiffness matrices of the virtual system. q tot = q T t δx T p T is the vector of all generalized coordinates of the virtual system.

Computation of the Stiffness and Mass Matrices of the Parallel Robot

The model of the virtual tree structure and of the free moving platform does not consider the closed-loop kinematic chains. As a matter of fact, the n qtot components of vector q tot are dependent. The independent components are gathered into vector q and their determination is described thereafter.

Determination of the Generalized Coordinates of the Parallel Robot

For determining one possible subset of generalized coordinates for the parallel robot, let us express the relations between the vector of generalized velocities qt i and the twist of the last element m i for each leg i. Using (28) for computing the twist

3 t i,m i = v T i,m i (A i,m i ) ω T i,m i (A i,m i )
T of the extremity of each leg, it comes:

t i,m i = J i v i,m i qt i (42)
As the leg extremity is also linked to the rigid platform, its twist can be related to the platform twist ẋp via the rigid body displacement relation:

t i,m i = J i p ẋp (43)
where

J i p = i,m i R 0 I 3 [p i ] × 0 I 3 (44)
in which J i p is a (6×6) matrix, [p i ] × is the cross product matrix of vector p i that characterizes the position of the attachement point C i,m i with respect to the platform centre position (Fig. 1(a)) and i,m i R 0 is the (6 × 6) rotation matrix between the global frame and the local frame attached to element i, m i , evaluated in the robot undeformed configuration.

Thus, expressing the twist t i,m i for each leg as a function of the platform twist ẋp and generalized coordinates qt i , the following set of equations is obtained:

   J 1 v 1,m 1 • • • 0 . . . . . . . . . 0 • • • J n vn,m n       qt 1 . . . qtn    -    J 1 p . . . J n p    ẋp = 0 (45)
which can be equivalently written as

J v qt -J p ẋp = J v -J p qt ẋp = J tot qtot = 0 (46)
where J tot is a (r n × n qtot ) matrix in the case of a spatial robot, n qtot > r n (r = 6 for a spatial robot, r = 3 for a planar robot). This means that a subset q d of r n variables in vector q tot is linked to the others. This subset is not unique. An idea could be to put in this subset all passive joints and platform variables, i.e., q * d = δq T p δx p . However, for over-constrained parallel robots, dim (q * d ) < r n. As a result, this vector should be completed using some other elastic variables that could be chosen arbitrarily. Meanwhile, it must be mentioned that most of parallel robots have identical legs and that such a methodology will lead to an asymmetrical description of the leg variables, which is not ideal. In order to avoid this problem, we had better put in q d the last r components q f t i of each vector q t i that is now decomposed into two parts:

q t i = q 0T t i q f T t i T .
Thus, variables q f t i are related to the others using (45):

-    J f 1 v 1,m 1 • • • 0 . . . . . . . . . 0 • • • J f n vn,m n       qf t 1 . . . qf tn    =    J 01 v 1,m 1 • • • 0 -J 1 p . . . . . . . . . . . . 0 • • • J 0n vn,m n -J n p         q0 t 1 . . . q0 tn ẋp      (47)
which can be rewritten as

-J f v    qf t 1 . . . qf tn    = J 0 v -J p      q0 t 1 . . . q0 tn ẋp      (48)
or also

   qf t 1 . . . qf tn    = -J f v -1 J 0 v -J p      q0 t 1 . . . q0 tn ẋp      =    J d 1,1 • • • J d 1,n J d 1,n+1 . . . . . . . . . . . . J d n,1 • • • J dn,n J d n,n+1    q (49)
where

• J 0i v i,m i (J f i v i,m i
, resp.) collects the columns of matrix J i v i,m i corresponding to variables q0 t i ( qf t i , resp.);

• J d ij is the matrix that maps q0 t j into qf t i , j = 1, . . . , n;

• J d i,n+1 is the matrix that maps ẋp into qf t i . It is noteworthy that the inversion of matrix J f v involves only the inversion of the (r × r) matrices

J f i v i,m i
, which is more efficient in terms of computational time. Moreover, when 3D beam elements are used for leg i, if the coordinates q f t i are the elastic coordinates of l-th element of this leg (previesly denoted as qe i,l ), it can be proven that, as the k-th column of matrix J f i v il corresponds to a unit twist that describes the displacement of the leg extremity due to the k-th coordinate of vector q f t i , J f i v il is equal to [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

J f i v il = i,m i R il i,m i R il [p il ] × 0 i,m i R il (50)
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 to element il and [p il ] × is the cross product matrix of the vector p il that characterizes the position of the leg extremity with respect to the frame linked to element il. Thus, its matrix inverse is equal to

J f i v il -1 = i,m i R T il -[p il ] × i,m i R T il 0 i,m i R T il ( 51
)
which requires few calculations and, before all, avoid any numerical inversion that could lead to numerical issues. If 2D beam elements are used, some similar relations can be obtained. Finally, the Jacobian matrix relating all variables qtot to the independent variables q = [ q0T

t 1 • • • q0T tn ẋp ] T can be obtained as: qtot =          q0 t 1 qf t 1 . . . q0 tn qf tn ẋp          =          I c 1 • • • 0 0 J d 1,1 • • • J d 1,n J d 1,n+1 . . . • • • . . . . . . 0 • • • I cn 0 J d n,1 • • • J dn,n J d n,n+1 0 • • • 0 I 6               q0 t 1 . . . q0 tn ẋp      = J q (52)
where I c i is the (c i ×c i ) identity matrix, c i being the dimension of vector q0 t i . Under the assumption of small displacements, the following relation holds true:

q tot = Jq (53)

Computation of the Natural Frequencies of the Parallel Robot

Introducing ( 52) and ( 53) into (41) leads to:

L = 1 2 qT J T M tot J q -q T J T K tot Jq = 1 2 qT M q -q T Kq (54)
Since in the natural frequency problem, matrices M and K are evaluated in the robot undeformed configuration, namely, for q e = 0, and as a result for δq a = 0 and δq p = 0, it turns out that the Lagrange equations yield:

d dt ∂L ∂ q - ∂L ∂q = Mq + Kq = 0 (55)
A solution q l of this equation satisfies:

ω 2 l M -K q l = 0 (56)
where ω l = 2πf l , f l is the natural frequency associated with the l-th natural mode of vibrations and q l is its associated eigenvector. Therefore, the natural frequencies of the parallel robot are evaluated by solving the following eigenvalue problem:

det

ω 2 l M -K = 0 (57)
In the next section, the natural frequencies of the NaVARo, a parallel robot developed at IRCCyN [START_REF] Rakotomanga | Kinetostatic performance of a planar parallel mechanism with variable actuation[END_REF], are computed using the proposed method and compared with the results obtained with the Cast3M software [START_REF] Castem | The castem software webpage[END_REF] and experimentally.

5 Case Study: Computation of the Natural Frequencies of the NaVaRo

Description of the NaVARo

The NaVARo (acronym for Nantes Variable Actuation Robot) was developed at IRCCyN and is shown in Fig. 3(a). The NaVARo is a 3-dof planar parallel manipulator composed of three identical legs and one moving platform made up of three segments E 1 P , E 2 P and E 3 P rigidly linked at point P . The i-th leg contains four links

A i B i , B i C i , C i E i , A i D i (named link 2i
, link 3i, link 4i and link 1i, respectively) connected with five revolute joints in such a way that

A i B i C i D i is a parallelogram linkage, i = 1, 2, 3.
The base frame F b (O, x 0 , y 0 , z 0 ) (not shown in Fig. 3(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. 3, F b and F p 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 moving platform, namely, the angle between x 0 and x p . q 1i denotes the angle between axis x 0 and link 1i. q 2i denotes the angle between link 1i and link 2i. Three double clutches are mounted to the base and located 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 eight actuation modes as described in [START_REF] Rakotomanga | Kinetostatic performance of a planar parallel mechanism with variable actuation[END_REF][START_REF] Arakelian | Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF]. Therefore, the moving platform can be moved throughout the (a) Prototype of the NaVARo located at IRCCyN, Nantes, France 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 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

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
γ 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.

Table 1

: MDH parameters of the i-th leg 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
Each link, of rectangular cross-section, is made up of duraluminum alloy (E = 74000 MPa, G = 28900 MPa, ρ = 2800 kg/m 3 ). Table 2 gives the cross-section area and the moments of inertia of the robot links. 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. The elasto-dynamic modeling of the NaVARo is complex due to the closed-loop chain in each leg and is obtained by following those three steps:

Table 2: Characteristics of the beam cross-sections link

A ij (m 2 ) I y ij (m 4 ) I z ij (m 4 ) I p ij (m 4 ) I 0 ij (m 4 ) 1i, 2i, 3i,
1. Computation of the mass and stiffness matrices of the virtual system assuming that the moving platform is cut at point P and the parallel linkages are opened at points D i , i = 1, 2, 3;

2. Computation of the mass and stiffness matrices of the legs including the closed-loop chains;

3. Computation of the mass and stiffness matrices of the NaVARo.

A single 3D beam element is used to model links 1i, 2i, 3i and 5i (see Section 2.4) while two 3D beam elements of equal lengths l (l = l C i D i = l D i E i ) are used to model links 4i. Decomposing links 4i into two beam elements is mandatory in order to be able to close the leg loops as mentioned previously in the step 2. Thus, the NaVARo is modelled as a spatial mechanism and its elastodynamic model contains 144 generalized coordinates: (i) 108 elastic coordinates; (ii) 12 passive joint coordinates, i.e., four passive joint angles per leg; (iii) 18 intermediary coordinates for the assembly of the legs; (iv) 6 coordinates for the moving-platform pose. From Sec. 4, it turns out that there are only 90 independent coordinates amongst those 144 coordinates.

Numerical analysis

A Matlab code was written to compute the robot mass and stiffness matrices using the modeling procedure presented in Secs. 2 to 4. The obtained robot mass and stiffness matrices were validated by means of an equivalent model developed using Cast3M software [START_REF] Castem | The castem software webpage[END_REF]. Cast3M aims to determine the elastodynamic model of structures modeled with beams. Both models give the same values for the first 90 natural frequencies of the NaVARo. Table 3 gives the first five natural frequencies of the NaVARo for the eight robot postures shown in Fig. 4.

The natural frequencies of the NaVARo are the same for poses 3, 5 and 7 (4, 6 and 8, resp.) as they correspond to a rotation of the robot base frame of ±120deg with respect to pose 3 (pose 4, resp.).

Experiments

Some experimental tests were carried out using the setup presented in Fig. 5. The application of experimental modal testing to the NaVARo was done through impact hammer excitation, a 3-D accelerometer response and data postprocessing, conducted using the DataBox software developed at IRCCyN and sold by MITIS company. The points and directions of excitation were chosen on points B i and E i of each leg along all axes in order to get the maximal number of resonance frequencies. Piezoelectric triaxial accelerometers with a sensitivity of 1000 mV/g were used to pick up the three acceleration responses. The acquisitions were performed for the eight robot postures shown in Fig. 4. Each measurement resolution is equal to 1 Hz as the acquisition time and the sampling time are equal to 1 s and 40 µs, respectively. The resonance frequencies were obtained with a fast Fourier transform of the signals given by the triaxial accelerometer. As a result, the measured resonance frequencies between 0 and 80 Hz for poses 1 to 4 are given in Table 4. As the results for poses 3, 5 and 7 (poses 4, 6 and 8, resp.) are similar due to the manipulator symmetry, only the results for poses 3 and 4 are given in Table 4 and the redundant poses were used to highlight some resonance frequencies with low energy level.

It is noteworthy that the resonance frequencies of the NaVARo amount to its natural frequencies as the damping is supposed to be negligible. 4 do not match with those shown in Table 3. As a matter of fact, the elasticity of the clutches has not been modelled and the joint masses have Thus, the natural frequencies of the NaVARo computed with this refined model and the measured frequencies are gathered in Table 4 by comparing the computed mode shapes with the hammer impact direction and the direction of the vibration signals, the latter being measured by the triaxial accelerometer.

We can notice that there is a good correlation between the measured frequencies and the computed natural frequencies.

Nevertheless, few predicted frequencies do not match with the measurements and vice-versa. Indeed, the theoretical and experimental results may differ due to the following reasons:

• The NaVARo has not been calibrated yet and there are some errors in the estimated moving platform pose;

• The passive joint elasticity has not been considered;

• The robot links are supposed to be coplanar in the theoretical model, whereas they are not in the prototype for collision avoidance;

• The robot links are not perfect beams as both ends are widened to insert ball bearings;

• The theoretical elastodynamic model does not consider any damping effect.

However, from those experiments, we can claim that the theoretical model is satisfactory and the proposed modeling procedure is efficient for reproducing the real behavior of any parallel robot.

Conclusion

Parallel robots have been increasingly used in the industry in the last few years and the characterization of their elastodynamic behaviour is still an issue. Accurate elastodynamic models of parallel robots are useful at both their design and control stages in order to define their optimal dimensions and shapes while improving their vibratory behaviour. Several models have been proposed in the literature. However, even if they can be adapted to any type of mechanism, they are not directly devoted to parallel manipulators and they do not provide a systematic and straightforward way for computing the Jacobian matrices associated with the kinematic constraints. Moreover, they do not take into account the symmetry in the robot leg description for choosing the independent coordinates describing the robot motion.

Therefore, a systematic method for the natural frequency computation of parallel robots has been developed in this paper. Indeed, the Jacobian matrices related to the kinematic constraints of the parallel robots are obtained in a straightforward way. Moreover, a way of choosing a symmetrical set of leg variables has been proposed. Contrary to most of the existing methods, the proposed approach does not contain any numerical matrix inversion, which is better to avoid numerical issues that may lead to a loss in the result accuracy.

This proposed approach has been used to compute the natural frequencies of the NaVARo, which is a planar parallel manipulator with multiple actuation modes developed at IRCCyN. The foregoing computed natural frequencies and those obtained with Cast3M software by using an equivalent robot model turned out to be identical. Then, some experiments have been carried out through impact hammer excitation and measurements of the platform displacements with a 3D accelerometer. The resonance frequencies obtained with a fast Fourier transform of the signals given by the triaxial accelerometer have been compared with the frequencies computed from a refined model of the robot. It appeared that there is a good correlation between the natural frequencies of the NaVARo computed with this refined model and the measured excitation frequencies.

Figure 1 :

 1 Figure 1: Schematic of a parallel robot for its dynamic modeling

Figure 2 :

 2 Figure 2: Schematics of the flexible elements into consideration

Figure 3 :

 3 Figure 3: The NaVARo

  4i 2.4 • 10 -4 1.152 • 10 -8 2.000 • 10 -9 1.352 • 10 -8 5.902 • 10 -9 5i 4 • 10 -4 3.333 • 10 -8 5.333 • 10 -8 8.666 • 10 -8 1.123 • 10 -8

  Pose 2 x = 0 m, y = 0 m, θ = -π/3 rad Pose 3 x = 0.117 m, y = 0.068 m, θ = -π/3 rad Pose 4 x = 0.182 m, y = 0.105 m, θ = -π/3 rad Pose 5 x = -0.117 m, y = 0.068 m, θ = -π/3 rad Pose 6 x = -0.182 m, y = 0.105 m, θ = -π/3 rad Pose 7 x = 0 m, y = -0.135 m, θ = -π/3 rad Pose 8 x = 0 m, y = -0.21 m, θ = -π/3 rad

Figure 4 :

 4 Figure 4: The eight poses used for the experiments

Figure 5 :

 5 Figure 5: Experimental setup: DataBox

Table 3 :

 3 Comparison of the natural frequencies obtained with Cast3M and the Matlab model.

	(Hz)	Pose 1 Pose 2 Pose 3 Pose 4 Pose 5 Pose 6 Pose 7 Pose 8
	f 1 (Cast3M)	44.10	45.71	36.98	40.17	36.98	40.17	36.98	40.17
	f 1 (Matlab model) 44.10	45.71	36.98	40.17	36.98	40.17	36.98	40.17
	f 2 (Cast3M)	44.10	45.71	49.31	50.32	49.31	50.32	49.31	50.32
	f 2 (Matlab model) 44.10	45.71	49.31	50.32	49.31	50.32	49.31	50.32
	f 3 (Cast3M)	53.98	54.58	53.37	52.99	53.37	52.99	53.37	52.99
	f 3 (Matlab model) 53.98	54.58	53.37	52.99	53.37	52.99	53.37	52.99
	f 4 (Cast3M)	60.63	65.35	67.28	67.36	67.28	67.36	67.28	67.36
	f 4 (Matlab model) 60.63	65.35	67.28	67.36	67.28	67.36	67.28	67.36
	f 5 (Cast3M)	95.62	97.92	91.80	91.52	91.80	91.52	91.80	91.52
	f 5 (Matlab model) 95.62	97.92	91.80	91.52	91.80	91.52	91.80	91.52

Table 4 :

 4 NaVARo natural frequencies (measured and computed using refined model) between 0 and 80 Hz.

	(Hz)	Pose 1 Pose 2 Pose 3 Pose 4
	f 1 meas.	22	19	17	18
	f 1 calc.	19.25	19.46	17.91	18.44
	f 2 meas.	24	21	19	20
	f 2 calc.	20.43	20.49	19.71	19.26
	f 3 meas.	32	-	23	22
	f 3 calc.	40.25	41.88	20.91	21.28
	f 4 meas.	-	44	27	33
	f 4 calc.	43.16	45.55	-	36.88
	f 5 meas.	42	45	32	43
	f 5 calc.	44.10	47.05	36.88	40.60
	f 6 meas.	50	53	43	44
	f 6 calc.	-	-	41.86	46.13
	f 7 meas.	52	54	46	50
	f 7 calc.	-	56.37	45.61	55.29
	f 8 meas.	62	56	48	56
	f 8 calc.	67.94	-	50.52	57.81
	f 9 meas.	66	60	57	58
	f 9 calc.	68.81	63.10	55.45	62.27
	f 10 meas.	77	-	60	66
	f 10 calc.	79.79	-	61.04	-
	f 11 meas.	-	-	61	-
	f 11 calc.	-	-	-	-
	f 12 meas.	-	-	65	-
	f 12 calc.	-	-	65.00	-

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

It is assumed that the generalized velocities are equal to d dt (δq t ) = qt

Note that index ij is written i, j in this section for a better understanding of the equations.

Acknowledgements

This work has been partially funded by the French ANR project ARROW (ANR 2011 BS3 006 01). The authors would like to thank the MO2P team of IRCCyN for the impact hammer and accelerometers loan. Dr. Mathieu Ritou and Mr. Jean-Pierre Regoin are also acknowledged for their great help and good advice during the experiments on the NaVARo.