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

Maxime Gautier

Global Identification of Joint Drive Gains and Dynamic Parameters of Parallel Robots

Keywords: Parallel robot, Force calibration, Drive gains, Inertial parameters identification, Total Least Squares

Off-line robot dynamic identification methods are based on the use of the Inverse Dynamic Identification Model (IDIM ), which calculates the joint forces/torques (estimated as the product of the known control signal -the input reference of the motor current loop -with the joint drive gains) that are linear in relation to the dynamic parameters, and on the use of linear least squares technique to calculate the parameters (IDIM-LS technique). Most of the papers dealing with the dynamic parameters identification of parallel robots are based on simple models, which take only the dynamics of the moving platform into account. However, for advanced applications such as output force control in which the robot interaction force with the environment are estimated from the values of the input reference, both identifications of the full robot model and joint drive gains are required to obtain the best results. In this paper a systematic way to derive the full dynamic identification model of parallel robots is proposed in combination with a method that allows the identification of both robot inertial parameters and drive gains. The method is based on the total least squares solution of an over-determined linear system obtained with the inverse dynamic model. This model is calculated with available input reference of the motor current loop and joint position sampled data while the robot is tracking some reference trajectories without load on the robot and some trajectories with a known payload fixed on the robot. The method is experimentally validated on a prototype of parallel robot, the Orthoglide.

Introduction

Parallel robots have increasingly been used for a few decades. This is due to their main advantages compared to their serial counterparts that are: (i) a higher intrinsic rigidity, (ii) a larger payload-to-weight ratio and (iii) higher velocity and acceleration capacities [START_REF] Merlet | Parallel Robots[END_REF]. In order to obtain these interesting properties, a good controller should be implemented. Several approaches could be envisaged [START_REF] Amiral | Design and control of a new six-dof parallel robot: application to equestrian gait simulation[END_REF][START_REF] Vivas | Predictive functional control of a parallel robot[END_REF], but it appears that, for high-speed robots or when varying loads have to be compensated for (e.g. in pick-and-place operations or machining), computed torque control is generally used [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]. This approach requires an accurate identification of the dynamic model of the robot with the load [START_REF] Khalil | Identification of the payload inertial parameters of industrial manipulators[END_REF], which can be obtained if two main conditions are satisfied:

1. a well-tuned derivative band-pass filtering of actuated joints position is used to calculate the actuated joints velocities and accelerations, and 2. the values of actuator drive gains g τ are accurately known to calculate the actuator force/torque as the product of the known control signal (computed by the numerical controller of the robot), i.e. the current references, by the drive gains.

However, it is often difficult or impossible to obtain robot manufacturers' data on joint drive gain values. And if available, data are usually given with an uncertainty greater than 10%, thus leading to identification and force calibration errors. Thus, drive gains must be calibrated. Some methods have been developed in the past [START_REF] Restrepo | Calibration of drive chain of robot joints[END_REF][START_REF] Corke | In situ measurement of robot motor electrical constants[END_REF] to calibrate the drive train constituted by a current controlled voltage source amplifier with gain G i which supplies a permanent magnet DC or a brushless motor with torque constant K t coupled to the link directly or through gear train with gear ratio N . Each parameter was identified separately to then compute the drive gain g τ = N G i K t , but these procedures were very sensitive to measurement errors and time consuming, requiring heavy tests on the drive chain. This sensitivity to errors directly affects the accuracy of the output force estimation.

More recent works [START_REF] Gautier | New method for global identification of the joint drive gains of robots using a known inertial payload[END_REF][START_REF] Gautier | New method for global identification of the joint drive gains of robots using a known payload mass[END_REF] have proposed to identify the global joint drive gains g τ for each actuated joint separately by using a known payload fixed on the endeffector. By carrying out with the robot two types of trajectories (trajectories with a known payload and trajectories without it), it appeared that it was possible to identify the payload terms using each joint j equations. However, as the drive gain g τ j of joint j is unknown, the estimated payload terms are indeed equal to the payload value divided by g τ j . As the payload was known, it was thus possible to estimate g τ j .

In [START_REF] Gautier | New method for global identification of the joint drive gains of robots using a known inertial payload[END_REF], the approach requires the use of the ten payload inertial parameters estimated with CAD software which implies some errors on the estimated values. Another drawback is to estimate the drive gain of one joint using only data coming from the corresponding joint equation which implies the loss of information about the coupled data on the other joints. As a result, we only succeed to identify the first four robot drive gains over six actuated joints. [START_REF] Gautier | New method for global identification of the joint drive gains of robots using a known payload mass[END_REF] is an improvement using only the accurate value of the weighed payload mass but with the drawback of using a sequential procedure for the computation of the payload inertial parameters that involves the accumulation of errors in the estimation of the drive gain parameters. Moreover, as previously, the estimation of the drive gain of one joint was done using only data coming from the corresponding joint equation which implies the loss of information about the coupled data on the other joints. With this approach, we only succeed in obtaining good results for the first four robot joints.

In this paper it is proposed a new method for the global identification of all robot dynamic parameters, including joint drive gains, using the input reference of the motor current loop and the joint position sampled data while the robot is tracking some reference trajectories without load fixed on the robot and some trajectories with a known payload fixed on the robot whose inertial parameters are measured or calculated with a CAD software. Contrary to the previous works, all dynamic parameters and drive gains are calculated altogether as the total Least Square (LS ) solution of an over-determined system that takes into account the coupling between the robot axes.

This method is combined with the use of an Inverse Dynamic Identification Model (IDIM ) that takes into account all the robot link dynamic parameters. Indeed, most of papers dealing with dynamic parameters identification for parallel robots propose a simplified dynamic model that takes into account the platform and drive chain dynamics only [START_REF] Vivas | Experimental dynamic identification of a fully parallel robot[END_REF][START_REF] Guégan | Identification of the dynamic parameters of the orthoglide[END_REF][START_REF] Renaud | Kinematic and dynamic identification of parallel mechanisms[END_REF]. Only few papers propose a systematic computation of the full IDIM. The most interesting and comprehensive works are presented in [START_REF] Grotjahn | Identification of friction and rigid-body dynamics of parallel kinematic structures for model-based control[END_REF][START_REF] Diaz-Rodriguez | A methodology for dynamic parameters identification of 3-dof parallel robots in terms of relevant parameters[END_REF], where the authors propose methods for computing the IDIM based on Jourdain's principle or Lagrange multipliers. The authors of [START_REF] Grotjahn | Identification of friction and rigid-body dynamics of parallel kinematic structures for model-based control[END_REF] focus on the identification of the friction terms. They have proposed a way to decouple the identification procedure: first the inertial parameters are identified, then the friction terms. This decoupled identification is obtained through the use of special trajectories with constant velocities. This work however does not take into account the asymmetry in the friction coefficients. [START_REF] Diaz-Rodriguez | A methodology for dynamic parameters identification of 3-dof parallel robots in terms of relevant parameters[END_REF] is a pedagogical work on the identification of parallel robots which presents the identification procedure in detail. However, the method proposed to compute the base parameters may lead, for a parallel robot with identical legs (which is the case for most of parallel robots), to a set of base parameters which does not conserve the symmetry properties of the robot legs. Thus, the identification results are sometimes difficult to analyze. Furthermore, for all these works, the way to identify the drive gains is not treated and some Jacobian matrices, whose computation is not straightforward, are not clearly derived.

The present paper has two main goals:

to propose a systematic and straightforward procedure for the computation of the IDIM for parallel robots which is combined with a way to obtain a set of symmetric base parameters when robot legs are identical, -to present a new method for the global identification of all robot dynamic parameters, including joint drive gains.

The work will be decomposed as follows. First, in Section 2, the way to compute the IDIM for parallel robots is developed. A straightforward procedure is proposed for the computation of the Jacobian matrices required for closing the loops in the dynamic model. Section 3 presents the usual identification procedure. Moreover, some brief recalls are made for the computation of the base parameters of robots and the method is adapted to parallel robots so that, for a parallel robot with identical legs, it leads to a set of base parameters which conserves the symmetry properties of the robot legs. Then, Section 4 details the new procedure for the global identification of all robot dynamic parameters, including joint drive gains. In Section 5, the method is experimentally validated on a prototype of parallel robot developed at the IRCCyN: the Orthoglide [START_REF] Chablat | Architecture optimization of a 3-dof parallel mechanism for machining applications, the Orthoglide[END_REF].

A condensed version of this work has been presented in [START_REF] Briot | Global identification of drive gains and dynamic parameters of parallel robots -part 1: Theory[END_REF][START_REF] Briot | Global identification of drive gains and dynamic parameters of parallel robots -part 2: Case study[END_REF]. The present paper contains detailed proofs to enlighten the theoretical understanding of the method and gives additional experimental results to show its practical efficiency.

Inverse Dynamic Identification Model of Parallel Robots

A Systematic Procedure for the IDIM Computation

This work focuses only on parallel robots without actuation redundancy. A parallel robot is a complex multi-body system having several closed loops (Fig. 1(a)). It is composed of a moving platform connected to a fixed base by n legs, each composed of m i elements. It is considered here that there is one actuator per leg, but the method can be easily extended to robots with several actuators per legs.

For the computation of IDIM of parallel robots, a method similar to the one presented in [START_REF] Ibrahim | Inverse and direct dynamic models of hybrid robots[END_REF] is applied. The proposed method is decomposed into two steps:

1. first, all closed loops are virtually opened to make the platform virtually disassembled from the rest of the structure (Fig. 1(b)); each leg joint is virtually considered actuated (even for unactuated actual joints) so that the robot becomes a tree structure with a free body: the platform; the dynamic model of the tree structure and of the free platform is then computed using a systematic procedure based on the Newton-Euler principle, 2. then, the loops are closed using the loop-closure equations and the Lagrange multipliers, which involve the computation of robot Jacobian matrices.

In the following, the computation of the IDIM of the virtual tree structure and of the platform is recalled, and then a straightforward way to compute the Jacobian matrices for closing the loops is detailed.

IDIM of Tree Open Loop Robots

It is known that the complete rigid dynamic model of any open-loop tree structure can be linearly written in terms of a (n t × 1) vector with respect to the standard parameters χ st t [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] (n t denotes the total number of joints for the virtual tree structure), τ idm t (q t , qt , qt ) = IDM st t (q t , qt , qt )

χ st t (1)
where τ idm t is the (n t × 1) vector of the input efforts of the virtual tree structure, IDM st t is the (n t × n st t ) Jacobian matrix of τ idm t , with respect to the (n st t × 1) vector χ st t of the standard parameters given by χ T st t = [χ 1T st , χ 2T st , . . . χ n t T st ] and q t , qt , qt are the vectors of the joint positions, velocities and accelerations, respectively.

For rigid robots, the vector χ j st of link j is composed of 14 standard parameters described as:

C 2,1 C 3,1 C m 1 ,1 C m 1 -1,1 C 1,n C 2,n C 3,n C m n ,n C m n -1,n C 1,k C 2,k C 3,k C m k ,k C m k -1,k C 1,2 C 2,2 C 3,2 C m 2 ,2 C m 2 -1,2
(a) Kinematic chain (C j,k is the joint j of the leg k, and m k is the total number of joints for the leg k) xx j , xy j , xz j , yy j , yz j , zz j are the 6 components of the inertia matrix of link j at the origin of frame j, -mx j , my j , mz j are the 3 components of the first moment of link j, -m j is its mass, -ia j is the total inertia moment for rotor and gears of the drive train, -f v j , f s j are the viscous and Coulomb friction coefficients in the joint, respectively, and τ of f j = τ of f f s j + τ of f τ j is an offset parameter which regroups the current amplifier offset τ of f τ j and the asymmetrical Coulomb friction coefficient τ of f f s j .

In the same vein, the IDIM of the platform can be obtained as:

τ p x, t, ṫ = IDM p x, t, ṫ χ p (2)
where τ p is the (6 × 1) vector of platform reaction wrench, IDM p is the (6 × 10) Jacobian matrix of τ p , with respect to the (10 × 1) vector χ p of the platform inertial standard parameters1 and x, t, ṫ are the platform position, twist and acceleration screw, respectively. Several methods can be used to systematically derive these equations. Here, an algorithm based on the use of the modified Denavit-Hartenberg robot geometric description and the Newton-Euler principle is applied. This modeling is known to give the dynamic model equations in the most compact form [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF].

IDIM of Parallel Robots

The IDIM of the virtual tree structure and of the free moving platform does not take into account the closed loop characteristics of parallel robots: among all joint and platform coordinates q t and x of the virtual robot (Fig. 1(b)), resp., only a subset denoted as q is independent in the real robot (the actual actuated joints positions, that are indeed a subset of q t ). All these variables are linked through the loop-closure equations that can be obtained by expressing the (translational and rotational) displacement x k of the last joints of each leg located at C m k ,k (that belong to both the platform and each leg k -Fig. 1) in two different ways: (i) as a function of the platform coordinates x and (ii) as a function of all joint coordinates q t , such that:

f (x, q t ) =    x 1 (x) -x 1 (q t )
. . .

x n (x) -x n (q t )    = 0 (3)
The main problem with (3) is that, for most parallel robots, it is difficult to solve the forward kinematic problem (fkp) using these equations in a straightforward way. Therefore, it is better to express the reduced loop-closure equations that directly relate the displacements q of the actuated joints to the moving platform coordinates x which is almost simpler to obtain [START_REF] Merlet | Parallel Robots[END_REF]:

f p (x, q) = 0 (4)
and to solve then the reduced fkp that gives x as a function of q. Obviously, this problem can be also tedious, but:

the equations are simpler to solve than when using the equations (3), -if the problem cannot be solved because the obtained polynomial equations have a too high degree, a numeric procedure can be used [START_REF] Merlet | Parallel Robots[END_REF].

Then, once the values of x are found as a function of q, it is possible to introduce them into (3) in order to express all joint coordinates q t as a function of x, and thus of q. This problem is generally easy for usual parallel robots [START_REF] Merlet | Parallel Robots[END_REF] and, even for more complicated cases, can now be solved using advanced mathematical methods [START_REF] Pfurner | New Trends in Mechanism Science, chapter Implementation of a new and efficient algorithm for the inverse kinematics of serial 6R chains[END_REF].

Differentiating (3) and ( 4) with respect to time, the following expressions can be obtained:

A p v + B p q = 0 ⇒ v = -A -1 p B p q = J p q, (5)
A p v + Ȧp v + B p q + Ḃp q = 0 ⇒ v = -A -1 p ( Ȧp v + B p q + Ḃp q) (6)
and

J tk v -J k qt = 0 ⇒ qt = J -1 k J tk v, (7)
J tk v + Jtk v -J k qt -Jk qt = 0 ⇒ qt = J -1 k (J tk v + Jtk v -Jk qt ) (8)
where

A p = ∂f p ∂x ind T, B p = ∂f p ∂q J tk = ∂f ∂x ind , J k = - ∂f ∂q t (9)
with x ind a vector of the independent platform coordinates, that are a subset of x, T a transformation matrix between the platform twist and the time derivatives of the terms x ind [START_REF] Merlet | Parallel Robots[END_REF], and v a vector of the independent coordinates in the platform twist (dim v = dim x ind ≤ 6), i.e.

t = Dv (10)
In the case of robots with 6 dof, D is the identity matrix. In these expressions, it is to be noted that the matrix J k stacks all Jacobian matrices corresponding to the independent motions of the last joints of each leg due to the joint displacements of each serial leg and is thus a square matrix of dimension ((n (m i -1)) × (n (m i -1))), -the matrix J tk is a matrix of dimension ((n (m i -1)) × n) that can be obtained by considering the rigid body displacement of any point of the robot platform as a function of the platform twist, and in the case of parallel robots without actuation redundancy, the matrices A p and B p are square of dimension (n × n).

Finally, by introducing ( 5) and ( 6) into ( 7) and ( 8), the expressions of qt and qt as a function of q, q and q can be obtained.

It should be mentioned that all the previous expressions are valuable as long as the robot does not meet any singularity and as long as the number of actuators is the same as the number of platform dof to control. Singularity avoidance or crossing is not the main topic of this paper, and the reader should refer to [START_REF] Briot | Optimal force generation of parallel manipulators for passing through the singular positions[END_REF][START_REF] Nenchev | Dynamic analysis of parallel manipulators under the singularity-consistent parameterization[END_REF] for further developments. In the following of the Section 2.3, it is considered that all these matrices are regular.

To take into account the loop-closure constraints in the dynamic model of the parallel robot, Lagrange multipliers λ T = λ T 1 λ T 2 can be used [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] to compute the (n × 1) vector of the actuated joint force/torque τ idm of the closed-loop structure. τ idm can be obtained in relation of the Lagrange multipliers λ by

τ idm = [0, -B T p ] λ 1 λ 2 = -B T λ, ( 11
)
where λ is calculated from the relation:

J T k 0 -J T tk A T p λ 1 λ 2 = A T λ = τ idm t τ pr (12)
with A a square matrix of dimension ((n m i ) × (n m i )) and

τ pr = D T τ p (13)
where τ p is defined in (2) and τ pr is a subset of forces/moments in τ p that can be found through the use of the principle of virtual powers, which states that:

v * T τ pr = t * T τ p = v * T D T τ p (14)
In this equation, the superscript " * " stands for a virtual velocity. It should also be mentioned that, in ( 11) and ( 12):

λ 1 stacks the wrenches applied by the virtual tree structure on the platform at points C m k ,k , -λ 2 stacks the values of the norms of the wrenches due to the platform dynamics in the platform joints located at C m k ,k .

Thus, the second equation of the system (12) represents the platform equilibrium so that the loops of the parallel robots can be closed. Solving [START_REF] Guégan | Identification of the dynamic parameters of the orthoglide[END_REF] and using the right part of ( 5) and ( 7), it can be demonstrated that:

τ idm = J T t τ idm t + J T p D T τ p = J T t IDM st t χ st t + J T p D T IDM p χ p = J T t IDM st t J T p D T IDM p χ T st t χ T p T = IDM st (q, q, q) χ st (15)
where J t = J -1 k J tk J p . Finally, because of perturbations due to noise measurement and modeling errors, the actual force/torque τ differs from τ idm by an error e, such that:

τ = τ idm + e = IDM st χ st + e ( 16
)
where τ is calculated with the drive chain relation:

τ = v τ g τ =    v τ 1 0 0 0 . . . 0 0 0 v τ n       g τ 1 . . . g τ n    ( 17
)
v τ is the (n × n) matrix of the actual motor current references of the current amplifiers (v τ j corresponds to actuator j) and g τ is the (n × 1) vector of the joint drive gains (g τ j corresponds to actuator j) that is given a priori by the manufacturer's data or measured with special time-consuming and heavy tests [START_REF] Restrepo | Calibration of drive chain of robot joints[END_REF][START_REF] Corke | In situ measurement of robot motor electrical constants[END_REF]. Equation ( 16) represents the Inverse Dynamic Identification Model (IDIM ).

IDIM of Parallel Robots Including the Payload

The payload is considered as an additional link (denoted as link l) fixed to the robot platform [START_REF] Khalil | Identification of the payload inertial parameters of industrial manipulators[END_REF]. Only n kl among its ten parameters are considered to be known. The model [START_REF] Chablat | Architecture optimization of a 3-dof parallel mechanism for machining applications, the Orthoglide[END_REF] becomes:

τ idm = IDM st IDM ul IDM kl   χ st χ ul χ kl   (18)
where:

χ kl is the (n kl × 1) vector of the known inertial parameters of the payload; -χ ul is the (n ul × 1) vector of the unknown inertial parameters of the payload, where n ul = 10 -n kl -IDM kl is the (n × n kl ) Jacobian matrix of τ idm , with respect to the vector χ kl , -IDM ul is the (n × n ul ) Jacobian matrix of τ idm , with respect to the vector χ ul .

Usual Identification Procedure

This part proposes a procedure for the computation of the base parameters for parallel robots so that, for a parallel robot with identical legs, it leads to a set of base parameters which conserves the symmetry properties of the robot legs, and then presents some necessary recalls on the identification procedure.

Computation of the Base Parameters

The off-line identification of the dynamic parameters is considered, given measured or estimated off-line data for τ and (q, q, q), collected while the robot is tracking some planned trajectories. The model ( 16) is sampled at frequency f m in order to get an over-determined linear system of r f m equations and n st unknowns:

Y f m (τ ) = W st f m (q, q, q)χ st + ρ f m (19)
where (q, q, q) is an estimation of (q, q, q), respectively, obtained by sampling and band-pass filtering the measure of q [START_REF] Gautier | Dynamic identification of robots with power model[END_REF], ρ f m is the (r f m × 1) vector of errors, Y f m is the (r f m × 1) vector of the input torques/force, sampled at frequency f m and W st f m (q, q, q) is the (r f m × n st ) observation matrix. The force/torque τ is perturbed by high frequency unmodelled friction and flexibility force/torque of the joint drive chain which is rejected by the closed loop control. These force/torque ripples are eliminated with a parallel decimation procedure which low pass filters in parallel Y f m and each column of W st f m and resamples them at a lower rate, keeping one sample over n d . This parallel decimation can be carried out with the MATLAB decimate function, where the low pass filter cut-off frequency, ω f p = 2π0.8f m / (2n d ), is chosen in order to keep Y f m and W st f m in the same frequency range of the model dynamics. After the data acquisition procedure and the parallel decimation of [START_REF] Pfurner | New Trends in Mechanism Science, chapter Implementation of a new and efficient algorithm for the inverse kinematics of serial 6R chains[END_REF], we obtain the over-determined linear system Y(τ ) = W st (q, q, q)χ st + ρ (20

)
where ρ is the (r × 1) vector of errors, Y is the (r × 1) vector of the input torques/force and W st (q, q, q) is the (r × n st ) observation matrix.

In Y and W st , the equations of each joint j are sorted in order to regroup the equations of each joint altogether such that:

Y T = (Y 1 ) T , • • • , (Y n ) T , W st = (W 1 ) T , • • • , (W n ) T T
, where Y j and W j represent the r/n equations of joint j. It is to be noted that no error is introduced by the parallel filtering process in the linear relation [START_REF] Briot | Optimal force generation of parallel manipulators for passing through the singular positions[END_REF] compared with [START_REF] Pfurner | New Trends in Mechanism Science, chapter Implementation of a new and efficient algorithm for the inverse kinematics of serial 6R chains[END_REF]. In [START_REF] Gautier | Dynamic identification of robots with power model[END_REF], practical rules for tuning this filter are given.

The identifiable parameters are the base parameters which are the minimum number of dynamic parameters from which the dynamic model can be calculated [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]. The minimal dynamic model can be written using the n b base dynamic parameters χ as follows:

Y = W (q, q, q) χ + ρ ( 21
)
where W is a subset of independent columns in W st which defines the identifiable parameters. Several methods exist for the computation of these subsets (analytical [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] or numerical [START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF]). In this work, it is preferred to use a numerical method based on QR factorization.

There is an infinity of possible subsets of base parameters, as presented in [START_REF] Diaz-Rodriguez | A methodology for dynamic parameters identification of 3-dof parallel robots in terms of relevant parameters[END_REF][START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF]. In [START_REF] Diaz-Rodriguez | A methodology for dynamic parameters identification of 3-dof parallel robots in terms of relevant parameters[END_REF], the authors test different subsets (obtained via the SVD of the observation matrix [START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF]) and keep the one which leads to the best conditioning index of the observation matrix. Even if it is computationally efficient, this method has a drawback: for a parallel robot with identical legs, it can lead to a set of base parameters which does not conserve the symmetry properties of the robot legs. Obviously, for avoiding this problem, it can be set a priori that some parameters are equivalent, which involves to sum their corresponding columns in the observation matrix. However, as there can be some small variations in the parameters values due to the manufacturing process, it is worth to avoid this a priori regrouping and check it a posteriori on the identified values.

Here, a method that avoids these drawbacks is described. For presenting it, let us make some brief recalls on the computation of the base parameters via QR factorization.

The QR factorization of the matrix W st of (20) takes the form:

Q T W st = R 0 ( 22
)
where Q is a r × r orthogonal matrix and R is upper triangular.

If the absolute value |R kk | of the k-th component located on the diagonal of R is inferior to α (α is the numerical rank -different from 0 because of round-off errors -and can be chosen such that α = ǫ max |R jj |, where ǫ is a small coefficient depending on the level of perturbations in W st (due to noise measurement and error modelling) and max |R jj | is the largest diagonal absolute value of R [START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF]), the k-th column W st k of W st can be deleted. At the end of the procedure, (n stn b ) columns of W st have been deleted that corresponds to (n st -n b ) standard parameters removed from vector χ st to keep a set of n b base parameters χ.

Because the QR algorithm starts from the last columns to the first of W st , the (n st -n b ) standard parameters to delete are dependent on the ordering of the columns of that matrix. For serial robots, the matrix W st is build such that the columns with the smaller indices are those corresponding to the links closest from the base. Thus, using the previous algorithm, the parameters with the smallest influence (those of the wrist) are eliminated from the base parameters.

For parallel robots, to take into account the symmetry in the leg dynamic parameters, it is preferable to order the columns of W st such that:

W st r = W st p W st χ 1 ,1:n W st χ 2 ,1:n • • • W st χ n st leg ,1:n (23)
where n st leg is the number of standard parameters for one leg, matrix W st p is the observation matrix corresponding to the platform inertial parameters and matrices W st χ k ,1:n concatenates the columns of matrix W st corresponding to the parameters χ k that are a priori identical for the n legs. Then, (n st -n b ) columns of W st r can be deleted using the previous approach based on the QR factorization to obtain a new observation matrix W associated with a set of symmetrical base parameters denoted as χ.

Recalls on Least Squares Identification of the Base Dynamic Parameters (IDIM-LS )

Using the base parameters and tracking "exciting" reference trajectories, i.e. optimized trajectories that can be computed by nonlinear minimization of a criterion function of the condition number of the W st matrix [START_REF] Gautier | Exciting trajectories for the identification of the inertial parameters of robots[END_REF][START_REF] Swevers | Optimal robot excitation and identification[END_REF], a well-conditioned matrix W st can be obtained. Examples of such trajectories will be shown in Section 5. The LS solution χ of ( 21) is given by:

χ = W + Y, where W + = W T W -1 W T (24)
It is computed using the QR factorization of W.

Standard deviations σ χi can be estimated assuming that W is a deterministic matrix and ρ is a zero mean additive independent noise [START_REF] Gautier | Dynamic identification of robots with power model[END_REF], with a covariance matrix C ρρ such that

C ρρ = E ρρ T = σ 2 ρ I r (25)
E is the expectation operator and I r , the r × r identity matrix. An unbiased estimation of the standard deviation σ ρ is:

σ 2 ρ = Y -W χ 2 / (r -n b ) (26)
The covariance matrix of the estimation error is given by:

C χ χ = E (χ -χ) (χ -χ) T = σ 2 ρ W T W -1 (27)
σ 2 χi = C χ χ(i, i) is the i-th diagonal coefficient of C χ χ (27).
The ordinary LS can be improved by taking into account different standard deviations on actuated joint j equations errors [START_REF] Gautier | Dynamic identification of robots with power model[END_REF]. Data in Y and W of ( 20) are weighted with the inverse of the standard deviation of the error calculated from ordinary LS solution of the equations of joint j [22]

Y j = W j χ + ρ j (28)
This weighting operation normalizes the errors in [START_REF] Briot | Optimal force generation of parallel manipulators for passing through the singular positions[END_REF] and gives the weighted LS estimation of the parameters (IDIM-WLS ).

Payload Identification

In order to identify the dynamic parameters of both the robot and the payload, using the model [START_REF] Ibrahim | Inverse and direct dynamic models of hybrid robots[END_REF], it is necessary that the robot carries out two types of trajectories [START_REF] Khalil | General solution for the dynamic modeling of parallel robots[END_REF]:

1. trajectories without the payload, and 2. trajectories with the payload fixed to the end-effector.

The sampling and filtering of the model IDIM (18) can be then written as:

Y = W a 0 0 W b W ul W kl   χ χ ul χ kl   + ρ (29)
where -W a is the observation matrix of the robot in the unloaded case, -W b is the observation matrix of the robot in the loaded case, -W ul is the observation matrix of the robot corresponding to the unknown payload inertial parameters, -W kl is the observation matrix of the robot corresponding to the known payload inertial parameters.

Thus, these two types of trajectories avoid the regrouping of the payload parameters with those of the platform and allow their independent identification. The next section presents the way to identify both link inertial parameters and joint drive gains.

Global Identification of the Robot Dynamic Parameters and the Drive Gains

In the usual IDIM-LS or IDIM-WLS, accurate values of the drive gains are necessary to compute vector Y. However, it is often difficult or impossible to obtain manufacturers' data. And if available, data are given with an uncertainty greater than 10%, thus leading to identification and force calibration errors. Therefore, it is preferable to introduce the drive gains into the base parameters and to use the Total Least Squares Identification (IDIM-TLS ) procedure.

Details on the TLS identification method can be found in [START_REF] Van Huffel | The Total Least Squares Problem: Computational Aspects and Analysis[END_REF] and many papers of the same authors. This method has been applied in [START_REF] Gautier | Identification of inertial and drive gain parameters of robots[END_REF] for the identification of the drive gains and the dynamic parameters on a two dof serial robot but gave arguable results due to the lack of an accurate scale factor.

In this paper a major improvement is proposed: the scaling of parameters using the accurate value of an additional payload mass.

Total Least Square Identification of the Robot Dynamic Parameters and the Drive Gain (IDIM-TLS )

Introducing ( 17) into [START_REF] Ibrahim | Inverse and direct dynamic models of hybrid robots[END_REF], and keeping only the base parameters, we obtain:

τ = v τ g τ = IDM IDM ul IDM kl   χ χ ul χ kl   + e ( 30
)
where IDM is the Jacobian matrix that relates the base parameters χ to the input torques τ . Then, using the method presented in the Section 3.3, the sampling and parallel decimation of the model IDIM (30) can be written as:

Y = V τ a V τ b g τ = W a 0 0 W b W ul W kl   χ χ ul χ kl   + ρ (31)
where V τ a is the matrix of v τ samples in the unloaded case, V τ b is the matrix of v τ samples in the loaded case.

V τ i =    V 1 τ i • • • 0 . . . . . . . . . 0 • • • V n τ i    , V j τ i =     v j τ i,1
. . .

v j τ i,r/n     , with i = a, b (32)
in which v j τ i,k is the k -th sample of current reference for actuator j, V j τ i regroups all the samples of the current reference for actuator j.

Rearranging the terms inside, Eq. ( 31) becomes:

W tot χ tot = -W a V τ a 0 0 -W b V τ b -W ul -W kl χ kl     χ g τ χ ul 1     = ρ (33)
where W tot is a (r × c) matrix (c = (n b + n + n ul + 1)) and χ tot is a (c × 1) vector.

Without perturbation, ρ = 0 and W tot should be rank deficient to get the solutions λχ tot = 0 depending on a scale coefficient λ. However because of the measurement perturbations, W tot is a full rank matrix. Therefore, the system (33) is changed to the compatible system closest to (33):

Ŵtot χtot = 0 (34)
where χT tot = χT ĝT τ χT ul 1 and Ŵtot is the rank deficient matrix, closest to W tot with respect to the Frobenius norm, i.e. W tot minimizes the Frobenius norm W tot -Ŵtot F [START_REF] Gautier | Identification of inertial and drive gain parameters of robots[END_REF]. χtot is the solution of the compatible system (34) closest to (33).

Ŵtot can be computed thanks to the "economy size" SVD of W tot [START_REF] Golub | Matrix computation[END_REF]:

W tot = U diag (s i ) 0 V T ( 35
)
where U and V are (r × r) and (c × c) orthonormal matrices, respectively, and diag (s i ) is a (c × c) diagonal matrix with singular values s i of W tot sorted in decreasing order. Ŵtot is given by:

Ŵtot = W tot -s c U c V T c ( 36
)
where s c is the smallest singular value of W tot and U c (V c , resp.) the column of U (V, resp.) corresponding to s c . Then, the normalized optimal solution χn tot ( χn tot = 1) is given by the last column of V, χn tot = V c [START_REF] Gautier | Identification of inertial and drive gain parameters of robots[END_REF]. There is infinity of vectors χtot = λ χn tot which are solution of (34) depending on a scale factor λ. The unique solution χ * tot = λ χn tot for the robot can be found by taking into account that the last value χ * tot c of χ * tot must be equal to 1 according to (33), i.e. λ = 1/ χn tot c , with χn tot c the last value of χn tot .

Statistical Analysis

Standard deviations σ χi on the dynamic and sensor gains parameters, are estimated assuming that all errors in data matrix W tot are independently and identically distributed with zero mean and common covariance matrix C W W such that

C W W = σ2 W I r W ( 37
)
where I r W is the identity matrix of dimension (r × c) × (r × c).

An unbiased estimation of the standard deviation σW [START_REF] Van Huffel | The Total Least Squares Problem: Computational Aspects and Analysis[END_REF]:

σW = s c / √ r -c (38)
The covariance matrix of the estimation error is approximated by [START_REF] Van Huffel | The Total Least Squares Problem: Computational Aspects and Analysis[END_REF]:

C χ χ ≈ σ2 W 1 + χ1:c-1 2 2 ŴT tot 1:c-1 Ŵtot 1:c-1 -1 (39)
with χ1:c-1 the vector containing the (c -1) first coefficients of χ * tot and Ŵtot 1:c-1 a matrix composed of the (c -1) first columns of Ŵtot . Finally, σ 2 χi = C χ χ(i, i) is the i-th diagonal coefficient of C χ χ and the relative standard deviation %σ χri is given by: %σ

χri = 100σ χi / | χi |, for | χi | = 0.
In order to improve the estimation of χ * tot , the rows of W tot are weighted taking into account the confidence on the measurements. As proposed in IDIM-WLS (Section 3.2), to improve the TLS solution, each row corresponding to joint j equation is weighted by the inverse of σW j , i.e. the standard deviation corresponding to the data of the joint j equations. Moreover, to take into account that the confidence on data in V τ a,b is higher than for data in W a,b and W u,kl , the columns of V τ a,b could also be weighted. However, our experiments have shown that the results were not really improved; therefore this last weighting procedure was not used in the next section.

Discussion on the A Priori Knowledge of the Payload Parameters

The accuracy of λ depends on the accuracy of the knowledge χ kl of the payload parameters. The parameters can be obtained by measurements or estimated by CAD software. However, due to the possible errors in CAD values, the most accurate payload parameter is the mass value m l that can be accurately measured using a weighing machine. For the robot under study in the next section, the Orthoglide, only the payload mass appears in the model (the other payload parameters have no effects on the dynamics). However, for other types of robots such as serial six dof robots, our experiments have shown that, due to the error on the parameter estimation using CAD software, it is better to carry out the identification IDIM-TLS using the knowledge of the payload mass only [START_REF] Gautier | Global identification of drive gains parameters of robots using a known payload[END_REF]. These results are not shown here as they are out of the scope of the paper.

5 Case Study: the Orthoglide Robot

Description of the Architecture

The Orthoglide is a parallel robot with three translational dof composed of three identical legs (Fig. 2). Each leg is made of one linearly actuated foot linked at its extremity to a spatial parallelogram (Fig. 3(a)). The parallelogram is also attached to the mobile platform. The directions of the three linear actuators of the Orthoglide are orthogonal (Fig. 2(b)). This aims at creating a mechanism with a workspace shape close to a cube and whose behavior is close to the isotropy wherever it is located in its workspace [START_REF] Chablat | Architecture optimization of a 3-dof parallel mechanism for machining applications, the Orthoglide[END_REF].

For the remainder of the paper, it should be mentioned that the gravity field g is directed along the vector y 0 of Fig. 3(c).

For the dynamic modeling of the Orthoglide, as the bars of the parallelograms are parallel and as the prototype has been designed so that they are identical, it is possible to replace the long bars by only one equivalent bar linked at each extremity by two orthogonal pivots (Fig. 3(b)) without loss of generality. Using the previous assumption, the Orthoglide kinematics necessary for computing the dynamic model are defined using the modified Denavit and Hartenberg notation (MDH ) [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]. The geometric parameters of the virtual open-loop tree structure are given in Tables 1 and2. The platform and payload are considered as supplementary bodies, the payload being fixed on the platform. They are respectively numbered as bodies 6 and 7.

The MDH notation being well known, the parameters of Tables 1 and2 will not be defined here. For more information concerning the MDH parameters, the reader should refer to [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF].

Table 1 MDH parameters for the frames corresponding to the first body of each leg.

j i a(j i ) µ ji σ ji γ ji b ji α ji d ji θ ji r ji 1 1 0 1 1 0 0 0 0 0 q 11 1 2 0 1 1 π/2 a π/2 0 0 q 12 -a 1 3 0 1 1 0 a -π/2 0 -π/2 q 13 + a
Table 2 MDH parameters for the frames corresponding to the leg i.

j i a(j i ) µ ji σ ji γ ji b ji α ji d ji θ ji r ji 2 i 1 i 1 0 0 0 -π/2 0 q 2i 0 3 i 2 i 1 0 0 0 -π/2 0 q 3i 0 4 i 3 i 1 0 0 0 0 d 4 q 4i 0 5 i 4 i 1 0 0 0 π/2 0 q 5i 0
As proposed in Section 2, the loop constraints can be taken into account by using the robot loop-closure equations [START_REF] Merlet | Parallel Robots[END_REF] and the Lagrange multipliers [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF][START_REF] Moon | Applied dynamics[END_REF]. The loop-closure equations of the Orthoglide can be expressed as [START_REF] Chablat | Architecture optimization of a 3-dof parallel mechanism for machining applications, the Orthoglide[END_REF]

f = f T 11 , f T 21 , f T 12 , f T 22 , f T 13 , f T 23 T = 0 with (40)
f 1i =   x y z   -R i 0   d 4 cos q 2i cos q 3i + aδ i -d 4 sin q 3i + aδ i+1 q 1i -a δi+2 -d 4 sin q 2i cos q 3i + d 6   , f 2i = q 2i q 3i + q 5i q 4i , f p =   x 2 + y 2 + (z -q 11 -d 6 ) 2 -d 2 4 (x -q 12 -d 6 + a) 2 + y 2 + (z -a) 2 -d 2 4 x 2 + (y -q 13 -d 6 + a) 2 + (z -a) 2 -d 2 4   = 0, (41)
where R i 0 is the rotation matrix between the leg frame (O i , x 1i , y 1i , z 1i ) and the base frame (O 0 , x 0 , y 0 , z 0 ) (Fig. 3(c)), a, d 4 and d 6 are geometric parameters defined in Fig. 3(b) (a = 0.34m, d 4 = 0.31m and d 6 = 0.03m), δ i is a scalar whose value is equal to 1 if i = 3, 0 if not and δi = 1 ( δi = 0, resp.) if δ i = 0 (δ i = 1, resp.).

The matrices J k , J tk , A p , B p and all positions, velocities and accelerations that are necessary for computing the dynamic model of the real robot can be calculated from ( 40) and (41) using expressions of Section 2 in a straightforward way.

Identification Results

In this part, experimentations are performed and the dynamic model identification is carried out on the Orthoglide. The actuation of each foot is achieved by a rotary DC motor Sanyo Denki linked to a ball-screw. The actuators are powered by current source amplifiers Sanyo Denki. The global drive gains given by the manufacturer are equal to 637π ≈ 2001. Actuator positions and current references are provided by the controller.

The approach presented in Section 3 requires that the robot moves on two types of exciting trajectories:

1. a first type of trajectory without any payload and 2. a second type of trajectory with a payload.

It should be mentioned that those trajectories can be different.

An example of exciting trajectory obtained using the approach recalled in Section 3 is shown in Fig. 4. In our experiments, the payload mass has been measured with an accurate weighing machine (M 7 = 1.983kg ± 0.001kg).

Even paying attention to the choice of the exciting trajectories, some small parameters remain poorly identifiable because they have no significant contribution in the joint torques. These parameters have no significant estimations and can be cancelled in order to simplify the dynamic model. Thus parameters for which the relative standard deviation %σ χri is too high are cancelled to keep a set of essential parameters of a simplified dynamic model with a good accuracy [START_REF] Khalil | General solution for the dynamic modeling of parallel robots[END_REF]. The essential parameters are calculated using an iterative procedure starting from the base parameters estimation. At each step the base parameter which has the largest relative standard deviation is cancelled. A new IDIM-TLS parameter estimation of the simplified model is carried out with new relative error standard deviations %σ χri . The procedure ends when max (%σ χri ) /min (%σ χri ) < r σ , where r σ is a ratio ideally chosen between 10 and 30 depending on the level of perturbation in Y and W. In the following of the paper, this ratio is fixed to 10. Table 3 presents the identification results. Subscript "R" stands for the parameters that have been regrouped using the procedure presented in Section 3. Their expressions are listed below (for i = 1, 2, 3):

m 1 j R = m 1 j + ia 1 j + m 2 j -3.2258mx 3 j + m 3 j , for j = 1, 2 m 1 3 R = m 1 3 -0.10197τ of f 1 3 + m 2 3 -3.2258mx 2 3 + m 2 3 ia 1 3 R = ia 1 3 + 0.10197τ of f 1 3 zz 2 i R = zz 2 i + yy 3 i -0.31mx 3 i + yy 4 i my 2 i R = my 2 i + mz 3 i + mz 4 i f v 2 i R = f v 2 i + f v 5 i f s 2 i R = f s 2 i + f s 5 i τ of f 2 i R = τ of f 2 i -τ of f 5 i xx 3 i R = xx 3 i -yy 3 i + 0.31mx 3 i zz 3 i R = xz 3 i -0.31mz 4 i f v 3 i R = f v 3 i + f v 4 i f s 3 i R = f s 3 i + f s 4 i τ of f 3 i R = τ of f 3 i -τ of f 4 i m 6R = m 6 + 3 i=1 3.2258mx 3 i + m 4 i + m 5 i
As it can be observed, the same regroupings appears for each robot leg. There is a small difference concerning the regrouped parameter m 1iR : as actuated joint 3 is vertical, contrary to joints 1 and 2 that are horizontal, this joint must support the gravity effects applied on the foot of leg 3 (link 1 3 ). This is not the case for joints 1 and 2. Thus, this is the reason why parameters m 13R and ia 13R are identified separately.

In Table 3, the given values of standard deviation and percentage of error for the payload mass m 7 are those estimated by the use of the weighing machine. The results show that the identified drive gains are about 5% off from those given by the manufacturer.

In Fig. 5, the actuator torques calculated with the relation ( 17) τ = v τ ĝτ (where v τ is the measured motor current reference and ĝτ the vector of the identified drive gains) are compared with torques computed using the IDIM (21) τ idm = IDM χ (where χ are the identified dynamic parameters). It should be noted here that the trajectory used for plotting these figures is different from the trajectories used in the identification process, i.e. the model has been crossvalidated. It can be seen that the torques are well calculated using the identified IDIM.

In order to definitively validate our method, a second payload of 1.136 kg is mounted on the platform and a classical IDIM-WLS identification is performed [START_REF] Gautier | Dynamic identification of robots with power model[END_REF]. Two cases are considered:

-(Case 1) the drive gains used are those of the manufacturer and, -(Case 2) the drive gains are those identified in Table 3.

The following results are found:

-(Case 1) m 7 = 1.09 kg ±0.02kg kg, -(Case 2) m 7 = 1.14 kg ±0.02kg kg.
It is thus possible to conclude that the identification of the drive gain using IDIM-TLS leads to better payload estimation. Finally, to test its effectiveness, the identified dynamic model is used in a model based feedforward control (Computed Torque Control (CTC )). Tracking errors are compared with those obtained with four different control laws:

- - f v 21R 5.57 1.36 f v 22R - - f v 23R - - τ of f 21R 0.96 1.10 τ of f 22R - - τ of f 23R -0.47 2.16 f v 31R - - f v 32R - - f v 33R -5
1. control of the robot using a simple PID controller, mean values of the obtained results. It can be observed that CTC with identified gains leads to slightly better results, even if the difference with respect to the CTC with manufacturer's gains is not really relevant. This can be explained by the fact that, for CTC strategies with identified and manufacturer's gains, the computed torques are slightly the same, as the difference in the gains is very small (less than 5%). It can also be observed that when the gains have an error of 25%, the tracking errors are clearly higher. All these result show the necessity of identifying the drive gains. For doing so, the presented effective approach can be used: for calibrating the drive gains, it is only necessary to weigh the payload mass and to carry out exciting trajectories on the robot.

Finally, to conclude this section and to claim that the proposed method can be applied to any kind of robots (serial or parallel, with any types of dof ), we would like to mention that it has been experimentally tested on three 6 dof industrial robots (the Staübli TX-40 (2kg of payload) and RX-90 robots (about 10kg of payload) and the Kuka KR270 robot (270kg of payload)). Experimental results have either confirmed the drive gains values given by the manufacturers or have shown significant improvements of the identification of the drive gains values leading to better payload estimations.

Conclusion

This paper has presented a new method for the global identification of the parallel robot dynamic parameters including the whole gains of the total drive chain. The presented work has:

proposed a systematic and straightforward procedure for the computation of the IDIM for parallel robots combined with a way to obtain a set of symmetric base parameters when robot legs are identical, -shown a new method for the global identification of all robot dynamic parameters, including joint drive gains. Fig. 6 The trajectories for testing the control strategies (blue dotted line: actuator 1; green dash-dotted line: actuator 2; red dashed line: actuator 3). (c) Actuator 3 Fig. 7 Tracking errors for each actuator on trajectory 1 (blue full line: simple PID; red dotted line: advanced PID with velocity feedforward; cyan bold line: CTC with parameters identified using identified drive gains; black thin line: CTC with parameters identified using manufacturer's drive gains; green dash+dotted line: CTC with parameters identified using error on gains).

The identification method is easy to implement and does not need any special test or measurement on the components of the joint drive train. It is based on a IDIM-TLS technique using motor current reference and joint position sampled data while the robot is tracking some reference trajectories without load fixed on the robot and some trajectories with a known payload fixed on the robot endeffector. The ten inertial parameters are measured or calculated by CAD software. The method has been experimentally validated on a prototype of parallel robot with three translational dof, the Orthoglide. The results have shown that the identified gains allows the best payload reconstruction and the best accuracy along a desired trajectory. Fig. 8 Tracking errors for each actuator on trajectory 1 (blue full line: simple PID; red dotted line: advanced PID with velocity feedforward; cyan bold line: CTC with parameters identified using identified drive gains; black thin line: CTC with parameters identified using manufacturer's drive gains; green dash+dotted line: CTC with parameters identified using error on gains).

Fig. 1 A

 1 Fig.1A general parallel robot (the grey circles denote the actuated joints).

Fig. 2 Fig. 3

 23 Fig.2The Orthoglide robot.

Fig. 4

 4 Fig.4Example of exciting trajectory (blue dotted line: actuator 1; green dash-dotted line: actuator 2; red dashed line: actuator 3).

  Actuator accelerations for trajectory 2

Table 3

 3 Essential parameters of the Orthoglide.

	Par.	Val.	%σ χri	Par.	Val.	%σ χri	Par.	Val.	%σ χri
	gτ 1	2170	0.39	gτ 2	2100	0.39	gτ 3	2070	0.24
	m 11R	9.40	0.44	m 12R	8.95	0.44	m 13R	0.87	1.12
	ia 11R	-	-	ia 12R	-	-	ia 13R	8.05	0.25
	f v 11	85.7	0.42	f v 12	86.7	0.50	f v 13	84.0	0.27
	f s 11	35.2	0.43	f s 12	41.7	0.42	f s 13	24.1	0.32
	mx 21	-0.12	3.75	mx 22	-	-	mx 23

Table 4

 4 Tracking errors along the trajectories.

	Trajectory 1	Maximal error (µm)		rms (µm)
	Control strategy	joint 1	joint 2	joint 3	joint 1 joint 2	joint 3
	Simple PID	5165	5202	335	1755	1776	84
	Advanced PID	453	506	126	158	177	18
	CTC with manuf. gains	180	177	96	61	62	10
	CTC with identif. gains	179	177	96	60	61	10
	CTC with error on gains	188	184	93	66	66	10
	Trajectory 2	Maximal error (µm)		rms (µm)
	Control strategy	joint 1	joint 2	joint 3	joint 1 joint 2	joint 3
	Simple PID	9910	12610	10469	5307	6371	5856
	Advanced PID	727	919	842	420	449	429
	CTC with manuf. gains	305	339	383	140	149	156
	CTC with identif. gains	302	339	382	137	147	154
	CTC with error on gains	350	353	409	154	164	169

The number of standard parameters of a free rigid body can be reduced to 10 inertial parameters as it is not necessary to consider the parameters ia j , f v j , f s j and τ of f j that are related to actuated joint drive chains.

Acknowledgements The authors would like to thank Philippe Lemoine for his great help and support during the experimental validations carried out on the Orthoglide in the scope of this work.

This work has been partially funded by the French ANR project ARROW (ANR 2011 BS3 006 01).

2. control of the robot using an advanced PID controller with velocity feedforward, 3. control of the robot using a CTC with a dynamic model identified using the manufacturer's drive gains, 4. control of the robot using a CTC with a dynamic model identified applying 25% of error on the manufacturer's drive gains keeping the same bandwidth (100 rad/s) for the linear part of the five control laws, so that the results can be fairly compared. It should be mentioned that the last case is introduced into the comparison as the identified gains are very close to the manufacturer's gains, so that the effect of drive gains obtained with a larger error can be analyzed. Two different trajectories are tested. They are plotted in Fig. 6. A mass of 1.983 kg is fixed on the platform. The nominal value of this mass is directly introduced into the dynamic models used in the CTC. The tracking errors in each case for all the three actuated axes are plotted in Figs. 7 and8. It is obvious that CTC strategies lead to lower tracking errors and are about 40 times better than the simple PID control. However, it is difficult to analyze which CTC strategy is the best. Therefore, in Table 4 are presented the values of the maximal absolute tracking errors and the root-mean-square (rms) errors along the trajectories for each control strategy. Each trajectory is run twice and the presented values are the


Dinh Quan Nguyen

Sébastien Briot

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

ANALYSIS OF THE DYNAMIC PERFORMANCE OF SERIAL 3R ORTHOGONAL MANIPULATORS

Serial 3R orthogonal manipulators have been studied recently and it has been proved that they can exhibit good performances in term of workspace size and kinematic properties. The aim of this work is to analyze their dynamic performances, and compare them with anthropomorphic manipulators, which are very popular in industry.

Static and dynamic analyses based on the evaluation of the maximal input torques required for moving the manipulator are achieved. It is shown that, as in kinematics, the dynamic performances of the serial 3R orthogonal manipulators are better.

INTRODUCTION

Nowadays most industrial manipulators are based on anthropomorphic architectures (i.e. PUMA type). These manipulators are composed of a vertical revolute joint linked to the ground followed by two parallel joints and a spherical wrist. But other joint arrangements exist such as for the 3R orthogonal manipulators: any two consecutive joint axes are orthogonal. There are few examples of this type of manipulators such as the IRB 6400C by ABB-Robotics launched in 1998 (not longer commercialized) and the DIESTRO manipulator built at McGill University. But it appears that, due to the success of anthropomorphic manipulators, 3R orthogonal manipulators have almost disappeared.

Recently, researchers showed more and more interests in this kind of manipulators. It has been shown that these manipulators may have quite different kinematic properties according to their link lengths and joint offsets. The late studies, [START_REF] Wenger | Design of cuspidal and noncuspidal manipulators[END_REF]- [START_REF] Baili | Classification des manipulateurs 3R orthogonaux[END_REF], exposed deeply the analysis of geometric and kinematic properties of this manipulator family. For example the determination of the workspace boundary, the size and shape of the workspace, the existence of holes and voids, the accessibility inside the workspace, the feasibility of continuous trajectories in the workspace, singularity curves, and other properties like "cuspidality" (possibility of changing posture without crossing a singularity) were well-discussed.

In [START_REF] Zein | Analysis cinématique des manipulateurs sériels 3R orthogonaux et des manipulateurs parallèles plans[END_REF], the orthogonal manipulators and the PUMA-type manipulators were compared in terms of kinematic performances. The aim of this paper is to continue the comparison between these two manipulators based on the analysis on their dynamic performances.

Until now, the fields related to dynamic performances of a manipulator are still opened. The motion of a manipulator, a multi-body mechanical system, involves several effects caused by inertia, centrifugal, Coriolis, gravity and friction forces/ torques. Most of these are nonlinear effects, which makes very difficult to predict and control the dynamic performances of the manipulator. Several measures of dynamic performances of manipulators have been introduced in the literature. The simplest indices to define and understand are probably the maximal input torques or joint acceleration that are required when the manipulator is moving, or also the dissipated energy. But as they depend on the manipulator trajectory, general conclusions may be difficult to give.

In order to avoid these drawbacks, Asada [START_REF] Asada | A geometrical representation of manipulator dynamics and its application to arm design[END_REF] introduced the generalized inertia ellipsoid (GIE) as a tool to measure the capability of changing the end-effectors velocity in different directions for fixed kinetic energy. Yoshikawa [START_REF] Yoshikawa | Dynamic manipulability of robot manipulators[END_REF] proposed to use the dynamic manipulability ellipsoid (DME) for measuring the ease of changing the end-effector configuration by a set of joint torques with fixed magnitude (Euclidean norm). Both the GIE and the DME are based on the relationship between the generalized inertia force of the endeffector and the generalized inertia torques of joints. Tourassis and Neuman [START_REF] Tourassis | The inertial characteristics of dynamic robot models[END_REF] investigated the inertial characteristics of the dynamics equations of manipulators. They introduced the coefficients of coupling to measure the structural coupling of the dynamics equations. Khatib and Burdick [START_REF] Khatib | Optimization of dynamics in manipulator design: the operational space formulation[END_REF] derived a dynamic model in operational space, which is, in fact, the space spanned by the vector representing the pose of the endeffector. Using the operational space approach, they presented an optimization method for finding the largest and uniform bounds of end-effector's acceleration corresponding to the minimum joint torques at both low and high bounds of joint velocities. Graettinger and Krogh [START_REF] Graettinger | The acceleration radius: a global performance measure for robotic manipulators[END_REF] introduced another performance measure termed acceleration radius. For given bounds of joint torques, the corresponding acceleration radius defines the minimum upper bound of the magnitude of endeffector acceleration over the whole workspace. Other measures that attempt to capture manipulator performances as a function of the dynamics can be cited: Bowling and Khatib [START_REF] Bowling | The dynamic capability equations: a new tool for analyzing robotic manipulator performance[END_REF] proposed a general framework for capturing the dynamic capability of a general manipulator that includes the velocity and acceleration characteristics of the end-effector, taking into account factors such as torque and velocity limits of the actuators. Herman [START_REF] Herman | Energy-Based Indices for Manipulator Dynamics Improvement[END_REF] presented several indices for the evaluation of nonlinear effects during the motion of a manipulator based on the kinetic energy, based on the equations of motion expressed in term of inertial quasivelocities. Rao and Bhatti [START_REF] Rao | Probabilistic approach to manipulator kinematics and dynamics[END_REF] proposed a method to construct manipulator kinematic and dynamic model using a probabilistic approach to study the effect of various factors causing uncertainties in the behavior of a manipulator.

All these indices have the same advantage: they do not depend on the manipulator trajectory. However, they may be difficult to understand by engineers. The goal of this paper is to propose a general optimization procedure based on the analysis of the manipulator maximal input torques during motion, for characterizing the dynamic performances of the manipulators under study.

The paper is organized as follows. In the second part, the manipulators under study are presented. The third section shows results from static analysis. Then our optimization procedure is defined and applied to verify the dynamic performances. Last section concludes this paper.

DESCRIPTION OF THE MANIPULATORS UNDER STUDY

This section introduces the geometric properties and kinematic performances of the PUMA and 3R Orthogonal (ORTH) manipulators.

Serial 3R orthogonal manipulators

Geometric description

The orthogonal manipulators are positioning manipulators with three revolute joints in which the two pairs of adjacent joint axes are orthogonal. Figure 1 shows examples of the architecture of the manipulators and their modified Denavit-Hartenberg parameters [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] are presented in table 1.

Recent studies [START_REF] Zein | An exhaustive study of the workspace topologies of all 3R orthogonal manipulators with geometric simplifications[END_REF]- [START_REF] Zein | Analysis cinématique des manipulateurs sériels 3R orthogonaux et des manipulateurs parallèles plans[END_REF] have shown that, by simply setting certain geometric parameters to zero, it is possible to achieve different sub-families for the orthogonal manipulator with different kinematic performances, especially in their workspace topology (Fig. 2).

It is possible to show that the PUMA manipulator exhibits similar properties. Therefore among these many sub-families of each manipulator, it is necessary to choose the most suitable ones in order to perform the comparison between them.

Kinematic performances analysis

The kinematic performances of 3R ORTH manipulator have been studied in [START_REF] Zein | Analysis cinématique des manipulateurs sériels 3R orthogonaux et des manipulateurs parallèles plans[END_REF] based on the analysis of the regular dexterous workspace (RDW).

RDW is a part of the working areafree of singularities and with a regular shape (cube, square…)within which, the dexterity is better than a minimal prescribed value λ. The dexterity can be assessed with the reciprocal of the conditioning index k -1 of the Jacobian matrix J: In order to calculate the RDW, for example for a square shape, one can first find the maximal square inscribed in the workspace cross section (Fig. 3a then find the maximal sub-square where k -1 >λ where λ is a prescribed value that is up to the designer (Fig. 3b).

In [START_REF] Majou | Analysis cinétostatique des machines parallèles à translations[END_REF] a performance index was defined that accounts for workspace compactness and dexterity:

Case G r2=0

Case H d3=0

Case J d3=0

(2) It is shown in [START_REF] Zein | Analysis cinématique des manipulateurs sériels 3R orthogonaux et des manipulateurs parallèles plans[END_REF] that the ORTH manipulator with the maximum value of performance index  has the following proper geometric parameters:

(3) (1) (4)
2 3 3 2 4 0, 0, 0, d d r r d    
Its forward kinematic equations are given by:

        44 44 4 3 2 1 1 3 1 3 2 1 1 3 1 3 2 where cos , sin ii x d C C C d S S y d C C S d S C z d C S Ci Si               Its inverse kinematic equations are given by     2 2 2 3 2 4 2 4 1 asin 1 2 solutions 2 asin 2 solutions 3 3 2 1 3 atan2 3 2 1 3 x y z d z dC C C y S x C C x S y                                  
In this case, because of the simplicity in geometric structure, four distinct Inverse Kinematic Solutions (IKS) of ORTH manipulator result into two groups with equivalent postures (Fig. 5).

The selected optimized ORTH manipulator belongs to type C

with the performance index  =0.56:

Description of the PUMA manipulator

Figure 6 shows examples of the architecture of the manipulators and their modified Denavit-Hartenberg parameters [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] are presented in table 2.

Applying the same procedure for optimizing the kinematic performances of the PUMA manipulator leads to quite interesting results. The PUMA manipulator has very close performance index value and the same geometric parameters (link lengths) as the optimized ORTH one [4]:

2 2 3 3 4 0, 0, 0, d r r d d    
Its forward kinematic equations are given by

           
x d C d C C r r S y d C d C S r r C z d S d S S q q x d C C C y d C C S z d S S                                  and its inverse kinematic equations are     1 2 2 2 3 2 4 2 atan2 2 solutions asin 1 2 solutions 2 1 1 3 3 atan2 1 3 1 3 y x x y z d zC C xS zC S x C                                     
In this case, the PUMA manipulator also has 4 IKS solutions and results into two groups of equivalent postures.

For easiest notation later, we will indicate the two groups of equivalent postures of each manipulator as follows: 1 st +2 nd IGM ORTH and 1 st +2 nd IGM PUMA (here, 1 st group of PUMA manipulator has the elbow-up posture, 2 nd group has the elbow-down posture -Fig. 5).

The workspace of ORTH manipulator is a torus (simply rotate a sphere with radius d 4 around z-axis of the first joint) while PUMA's workspace is a sphere with radius 2×d 4 (Fig. 7). The selected optimized PUMA manipulator has the performance index  =0.6.

Manipulator choice

The mass and link-length parameters for the selected ORTH manipulator and PUMA manipulator are chosen such that they are similar to the parameters of first three links of the real

(5) (6) (8) (7)

Dynamic modelling of ORTH and PUMA manipulators

The inverse dynamic model (IDM) of a rigid manipulator composed of n moving links calculates the motor forces/torques vector, as a function of the generalized coordinates and their derivatives. It can be obtained from the Newton-Euler or the Lagrangian equations [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

 Γ Γ Γ
where:

dyn Γ is the   1 n vector of dynamic forces/torques due to the inertial, centrifugal, Coriolis and gravitational effects:

      , dyn    Γ M q q C q q q Q q ,,
qqqare respectively the   1 n vectors of generalized joint positions, velocities and accelerations,

  Mq is the   nn  manipulator inertia matrix,   , C q q is the   nn  matrix of centrifugal and Coriolis effects,   Qq is the   1 n vector of gravitational forces/torques,
  sign f v c coff     Γ F q F q F
where:

v F is the   nn  diagonal matrix of viscous parameters, c
F is the   nn  diagonal matrix of dry friction parameters, and sign(.) denotes the sign function, F coff is the   1 n vector of asymmetrical Coulomb friction force/torque between positive and negative velocities.

Thus [START_REF] Herman | Energy-Based Indices for Manipulator Dynamics Improvement[END_REF] becomes:

  sign dyn idm v c off     Γ Γ F q F q Γ
where:

off Γ is the   1 n vector of offset force/torque that regroups the amplifier offset and the asymmetric Coulomb friction coefficient, idm Γ defines the motor force/torque which includes the amplifier offset component, Then the inverse dynamic model (IDM) is given by:

    , idm  Γ M q q N q q  
, N q q regroups the   1 n vector of centrifugal, Coriolis, gravitational, friction and offset forces/torques:

        , , sign v c off      N q q C q q q Q q F q F q Γ
Computing these dynamic models is generally tedious, therefore a software, namely SYMORO+ [START_REF] Khalil | SYMORO+: A System for the Symbolic Modelling of Robots[END_REF] has been used. In the remainder, the static and dynamic models of ORTH and PUMA manipulators are computed thanks to SYMORO+.

STATIC ANALYSIS

In this section, the static analysis of the two manipulators is presented. This analysis gives a first insight of the manipulator performance but is also used for the optimization procedure presented in the next part.

The static analysis is a special case where the actuator torques are computed from the IDM with  q0 and .  q0 As a result, the input torques value only depends on the manipulator configuration. Therefore, for having a general insight of the manipulator performances in static, the input torques are good indicators.

Because of symmetry, it is only necessary to consider a cross-section of the workspace (a circle with radius R = d 4 ). Figure 8 shows an example of the static torques profile (absolute values) of the two manipulators (with payload = 5kg) for the first IKS solution.

Figures 9 and10 show the isocontour plots in the workspace of the maximum and minimum static torque between the two IKS of each manipulator. Then the changes of maximum torque absolute values with respect to the changes of the payload (for all IKS) are shown in figure 11. In these figures, T1,T2,T3 are the actuator torques of joint 1, 2 and 3 of each manipulator respectively. The static analysis gives much information:  For the ORTH manipulator, static torque profiles for maximum and minimum of absolute values of the input torques are similar (this was expected since the two postures of the ORTH manipulator are close while those of the PUMA manipulator are quite different)  The joint torques 2 and 3 of ORTH manipulator are distributed quite equally. Meanwhile for the PUMA, the torque values of joint 2 are much larger (more than twice than joint 3) To sum up, in static analysis, the ORTH manipulator exhibits better performance than the PUMA manipulator.

DYNAMIC ANALYSIS

As presented in the introduction, the main drawback of analyzing the dynamic performances with the manipulator input torques is that they depend on the trajectory. In order to avoid this problem, one idea could be to make the manipulator moves on selecting "exciting" trajectories that make it exhibiting the highest values of the input torques. To generate these exciting trajectories, it is necessary to define an optimization procedure that is detailed below. In this section, the following assumptions are made:  The analysis is performed under the same trajectories in the task-space for both manipulators, with the same payloads and under the same maximum operational velocity V = const (magnitude of Tool Control Point's (TCP) velocity in the task-space).  The traveling time between any two points A and B in the task-space is known.  The manipulator keeps track of the trajectory perfectly, i.e. the controller is considered as ideal. Hence, all the terms are computed with desired joint positions, velocities and accelerations ( , , ) ddd qqq .  In the dynamic models, we do not consider the effects of frictions (FV, FS).  Because ORTH manipulator has two equivalent postures which are close to each other (Fig. 5) and the PUMA

Optimization procedure Trajectory generation

Several methods for exciting trajectory generation exist (most of them have been developed for dynamic parameter identification) [START_REF] Pressé | Identification des Paramétres Dynamiques des Robots[END_REF]. It is proposed here to adapt one of them in order to generate the desired exciting trajectories that enable the manipulators to attain their maximal input torques.

A desired trajectory is composed of several segments which are determined by a subset of optimized points. Each via-point on the trajectory are constrained to belong to different areas in the workspace. These areas are chosen thanks to the static analysis that defined some zones where the static torques are high. It is necessary to choose different zones that allow covering all the workspace. Hence a suitable method to put the desired points is needed. Two methods can be used to link these points :  Point-to-Point: the trajectory will be constructed from several segments, i.e. P 1 -P 2 , P 2 -P 3 , P 3 -P 4 , P 4 -P 5 , etc.

The manipulator starts at the beginning of each segment and stops at its extremity. Different motion profiles can be applied on these paths (5 th order polynomial laws, bangbang profiles, etc.)  Via-Points: the trajectory is constructed using a continuous spline defined by a subset of several points, i.e. P 1 -P 2 -P 3 -P 4 -P 5 -...-P n . The manipulator starts at P 1 and only stops at point P n . In the remainder of the paper, cubic splines are used [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF].

Choosing an optimization tool

All the analysis steps were carried out by using the software MATLAB which provides several optimization toolboxes for solving single-objective problems or multiobjectives problems (i.e. fmincon, fgoalattatin...).

In this case, we choose fgoalattain function to perform the optimization tasks since several objectives will be considered.

Optimization algorithm

With the purpose of comparing, the optimization procedure is performed in order to define the trajectories in the task-space of the ORTH or PUMA manipulator. At first, it is important to select which posture of each manipulator will be compared. After achieving the optimized profile (that consists of several points and travelling time) for one manipulator, the corresponding inverse-solution of the other manipulator can be found.

By doing so, the constraints of the end-effector's position inside the ORTH workspace can be set in order not to cross the singularity boundary (especially for PUMA manipulator). Checking whether or not the trajectory is inside the working space is also needed for desired via-points which are close to the boundary.

The procedure to construct exciting trajectories is detailed as follows: Objective:

  k max max, for 1, 2,3 j k abs T j  
where the subscript k corresponds to the k th segment of the trajectory (Fig. 12).

Searching parameters:

= ( , ) 4 -1 X Nm  X Ps t Assume that the optimized trajectory has m points whose coordinates are stored in the vector Ps. Then the decision variables X is a column vector with the length of N X (the first (m-1) elements are travelling time and the following 3m elements are coordinates of m via-points in the task-space). Constraints:

              ikmax imax max Σ min max min max min max T T i = 1, 2,3 k = 1, 2, , 1 t
travelling time must be positive x, y, z ORTH's workspace (constraint limits in task-space)

k k k k k k k k k k k k m                          L VV 0 t P x x x y y y z z z
where:

N X is the total number of searching parameters, T imax is the maximal torque value of joint i th , T ikmax is the maximal torque value of joint i th at segment k th , V max is the (1 × m-1) vector of maximal velocities of the endeffector at all segments, L(k) is the total travelling length of the end-effector at segment k th :

      k = n +1 -n  L P P     Σ t = k
is the total traveling time at segment , , , is position of the end-effector in the task-space.

th k x y z

t P Determine the limit ranges

Choosing initial starting points is one important factor of the optimization procedure. Here, the static torques profiles (see section 3) are used to determine ( min max , kk

XX).

The idea is to select the limit ranges that could cover the entire workspace and also that could pass through the most exciting zones defined by static analysis (i.e. from low distributing torques areas to high distributing torques areas). One way to do so is to divide the workspace into small elements (i.e. small cubes).

The workspace can be divided into layers (separated by the segments 1 to 32 on Fig. 13) regrouped in 4 regions (I,II,III,IV). Each layer is a cross-section of the ORTH's workspace and is indicated with an index i (i = 1,2,…,32indices on the largest external circle of the torus). A crosssection is divided into 32 squares. Each square is indicated with an index j (j = 1,2,…,32). A center point P(i,j) of a searching zone is defined as the center of the corresponding square with index j on the layer with index i. Coordinates of this point in the task-space can be obtained by:   A trajectory with desired passing locations and directions can then be constructed by choosing suitable indices (i k ,j k ) for any via-point P k .

In optimization, the starting points are initialized as random values within offsets around the center points P k (x k ,y k ,z k ). Then the procedure will search for optimized points inside the limit cubes defined by (18).

Optimization results

Four different optimized trajectories are generated using 5 via-pointsi.e. 4 segments -(total 19 searching parameters) depending on the four regions of the workspace (I,II,III,IV). In table 5:  The matrix P contains the initial parameters used for the localization of the center via-points. Initial parameters of point P k are given by the k th column of P. The first row contains indices for layersi k , the second row contains indices for squaresj k .  The k th term of vector Δ represents the size Δ k of the searching interval for the final position of point P k (Fig. 13).  The maximum torques (T1, T2, T3) are chosen based on the static torque profiles with a payload of 5kg. Figure 13 shows an example of constructing the range's limit for any trajectory with desired locations and directions. For region I, we have: In this test, the optimization procedure for the ORTH manipulator is first performed to bring out its performance in the worst possible cases (reach constraint boundaries of maximum torques or maximum speed). Then, the same trajectory of the TCP is applied to the PUMA manipulator. Figures 141516show the results of the optimized trajectories with joint accelerations and joint torques of each manipulator by using cubic spline function method (via-points method). Test 2: optimization using 5 th degree polynomial (point-topoint): optimization for PUMA manipulator first

In this test, the optimization procedure for the PUMA manipulator is first performed to bring out its performance in the worst possible cases. Then, the same trajectory of the TCP is applied to the ORTH manipulator. In term of joint torques, the ORTH manipulator has better performance. The joint torques values of the PUMA in both cases are larger, especially much larger in case of joint 2 (Fig. 16,19). In term of joint accelerations, both manipulators have fair performances (Fig. 15,18). Because the torque profiles for each manipulator are complex surfaces (the maximal torque values are varied while changing either the speed Vs or payload), the contours plots of the differences between maximum torque values of the ORTH and the PUMA are presented. In the areas that have the negative values, the ORTH manipulator performs better than the PUMA manipulator.

limit cube of center via-point Pk

P k 2Δ k i k ,j k k k k k k k k k k k k k x x x y y y z z z                        
In practice, it is enough to consider the areas with a velocity of the end-effector less than 1m/s. The performances of joint torques 2 and 3 of the ORTH manipulator are better than for the PUMA manipulator. Although joint torque 1 of the ORTH manipulator is larger than for the PUMA, the differences are small with to the differences of the second joint torques.

CONCLUSIONS

The aim of this work was to analyze the dynamic performances of 3R orthogonal manipulators, and to compare them with anthropomorphic manipulators, which are very popular in industry.

Static and dynamic analyses based on the evaluation of the maximal input torques required for moving the manipulator have been achieved. Moreover, a systematic procedure for analyzing the dynamic performances of manipulators in term of input torques has been proposed. It has been shown that the dynamic performances of the serial 3R orthogonal manipulators are better than those of the PUMA manipulator.



  are the maximum and minimum singular values of J.

  d is the Chebyshev distance from the center point M 0 to a singular point A s :

Figure 4 .

 4 Figure 4. Optimized ORTH manipulator, (r3=0, r2 = d4 , d2 = 0)

Figure 5 .

 5 Figure 5. Postures of ORTH_VER and PUMA at an example point in the workspace

Figure 3 .Figure 6 .Figure 7 .

 367 Figure 3. RDW with k -1 min=0.25 (ORTH type C: d4 = 1.5 , r2 = 1) [4]

Table 3 .Γ

 3 Dynamic parameters of ORTH manipulator: payload = is the   1 n input electromagnetic force/torque vector of the drive chain, in SI units on the joint side, f Γ is the   1 n vector of the loss force/torque due to frictions and motor iron losses, eddy currents and hysteresis effect. Usually, it is approximated with a viscous friction coefficient and a Coulomb friction force/torque:

Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .

 89101112 Figure 8. Profiles of the static torques absolute values for ORTH and PUMA manipulators for the first IKS solution (payload = 5kg)

,

  P 0,j (x,0,z) is the center point of the j th square on the layer with index 0/32, and Rot(z, i  ) denotes the rotation matrix along the z-axis by an angle i  .

  This means that the center desired via-points are: P 1 (1,22) -P 2 (7,16) -P 3 (3,5) -P 4 (2,27) -P 5 (6,13) Test 1: optimization using Cubic Spline function (viapoints): optimization for ORTH manipulator first

Figure 13 .Figure 14 .

 1314 Figure 13. Constructing limit ranges for desired trajectory

Figure 15 .

 15 Figure 15. Joint Accelerations : Cubic Spline, task-space (Test 1optimization for ORTH first)Table 5. Initial parameters for optimization procedure

Figure 16 .Figure 17 .Figure 18 .Figure 19 .Figure 20 .

 1617181920 Figure 16. Actuator Torques : Cubic Spline, task-space (Test 1optimization for ORTH first)

Figure 21 .

 21 Figure 21. Difference in maximum of the absolute values of the input torques between ORTH and PUMA manipulators (verifying with 5 th degree polynomial method for via-points in test 2)

Table 1 . Geometric parameters of ORTH manipulator

 1

	Joint	σ	α	d	θ	r
	1	0	0	0	θ 1	0
	2	0	-90 0	d 2	θ 2	r 2
	3	0	90 0	d 3	θ 3	r 3

Figure 1. Geometric structure of 3R orthogonal manipulators Figure 2. The ten families of manipulators with at least one parameter equal to zero [2]

dyn in f