Keywords: Identification, closed-loop output error, least-squares methods, robot dynamics

Off-line robot dynamic identification methods are mostly based on the use of the inverse dynamic model, which is linear with respect to the dynamic parameters. This model is sampled while the robot is tracking reference trajectories that excite the system dynamics. This allows using linear least-squares techniques to estimate the parameters. The efficiency of this method has been proved through the experimental identification of many prototypes and industrial robots. However, this method requires the joint force/torque and position measurements and the estimate of the joint velocity and acceleration, through the bandpass filtering of the joint position at high sampling rates. The proposed new method requires only the joint force/torque measurement. It is a closed-loop output error method where the usual joint position output is replaced by the joint force/torque. It is based on a closed-loop simulation of the robot using the direct dynamic model, the same structure of the control law, and the same reference trajectory for both the actual and the simulated robot. The optimal parameters minimize the 2-norm of the error between the actual force/torque and the simulated force/torque. This is a non-linear least-squares problem which is dramatically simplified using the inverse dynamic model to obtain an analytical expression of the simulated force/torque, linear in the parameters. A validation experiment on a 2 degree-of-freedom direct drive robot shows that the new method is efficient.

I. INTRODUCTION

HE usual identification method based on the inverse dynamic identification model (IDIM) and leastsquares (LS) technique has been successfully applied to identify inertial and friction parameters of several robotic prototypes and industrial robots [START_REF] Khosla | Parameter Identification of Robot Dynamics[END_REF], [START_REF] Atkeson | Estimation of Inertial Parameters of Manipulator Loads and Links[END_REF], [START_REF] Gautier | Identification of Robot Dynamics[END_REF], [START_REF] Kawasaki | Terminal-Link Parameter Estimation and Trajectory Contro[END_REF], [START_REF] Ha | An efficient estimation algorithm for the model parameters of robotic manipulators[END_REF], [START_REF] Raucent | Identification of Barycentric Parameters of Robotic Manipulators from External Measurements[END_REF], [START_REF] Prüfer | Identification of Robot Dynamics with Differential and Integral Models : a Comparison[END_REF], [START_REF] Gautier | Identification of the Dynamic Parameters of a Closed-Loop Robot[END_REF], [START_REF] Gautier | Dynamic Identification of Robots with Power Model[END_REF], [START_REF] Swevers | Optimal Robot Excitation and Identification[END_REF], [START_REF] Swevers | Dynamic model identification for industrial robots -Integrated experiment design and parameter estimation[END_REF], [START_REF] Gautier | Extended Kalman Filtering and Weighted Least-squares Dynamic Identification of Robot[END_REF], [START_REF] Khalil | Modeling identification and control of robots[END_REF], [START_REF] Vandanjon | Comparison of two identification techniques: Theory and application[END_REF], [START_REF] Hollerbach | Model Identification[END_REF], amongst others. Good results can be obtained provided a well-tuned derivative bandpass filtering of joint position to calculate the joint velocities and accelerations is used.

Another approach is to minimize a quadratic error between an actual output and a simulated output of the system, assuming both the actual and simulated systems have the same input. This is known as an output error (OE) identification method [START_REF] Richalet | The global approach in identification protocol optimization[END_REF], [START_REF] Walter | Identification of parametric models from experimental data[END_REF]. The optimal values of the parameters are calculated using non-linear programming algorithms to solve a non-linear least-squares problem. The output is given by a state-space model output equation, which is typically the joint position for mechanical systems. Difficulties arise from the choice of initial conditions, resulting in multiple, local solutions [START_REF] Tohme | A Methodology to Enhance the Convergence of Output Error Identification Algorithms[END_REF]. The OE method has been used to identify electrical parameters of a synchronous machine, and a A new closed-loop output error method for parameter identification of robot dynamics M. Gautier (1) , A. Janot (2) and P-O Vandanjon (3) T This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 2 comparison with the IDIM-LS method showed very similar results [START_REF] Khatounian | Parameters Estimation of the Actuator Used in Haptic Interfaces: Comparison of two Identification Methods[END_REF].

Both IDIM and OE methods require the joint position and the joint force/torque measurements.

The proposed new identification method needs only the joint force/torque measurements. It is based on a closed-loop simulation using the direct dynamic model while the optimal parameters minimize the 2norm of the error between the actual force/torque and the simulated force/torque, assuming the same control law and the same reference trajectory. This non-linear least-squares problem is dramatically simplified using the inverse dynamic model to formulate the simulated force/torque as an algebraic function linear in relation to the parameters. This paper describes the new identification method and experimental results obtained using a 2 DOF robot.

A condensed version of this work has been presented in [START_REF] Gautier | DIDIM: A new method for the dynamic identification of robots from only torque data[END_REF]. This paper contains detailed proofs to enlighten the theoretical understanding of the method and gives additional experimental results to show the practical efficiency of the method.

The paper is organized as follows: section II reviews the usual identification technique of the dynamic parameters of the robot. Section III presents the output error method. The new identification method is presented in section IV. The modeling of the SCARA prototype robot is presented in section V. This direct drive prototype is very well suitable for the study of the method because it emphasizes non linear coupling while it is divided by the squared high gear ratio for industrial robots. The experimental results are given in section VI. Finally, section VII is the conclusion.

II. IDIM: INVERSE DYNAMIC IDENTIFICATION MODEL TECHNIQUE

The inverse dynamic model (IDM) of a rigid robot composed of n moving links calculates the motor torque vector idm τ , 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], [START_REF] Featherstone | Dynamics[END_REF]. It is given by the following relation:

= ( ) + ( , ) idm τ M q q N q q && & (1)
Where q , q & and q && are respectively the ( )

x n 1 vectors of generalized joint positions, velocities and accelerations, ( ) M q is the ( )

x n n robot inertia matrix, and ( , )

N q q
& is the ( )

x n 1 vector of centrifugal, Coriolis, gravitational and friction forces/torques.

The choice of the modified Denavit and Hartenberg frames attached to each link allows a dynamic model that is linear in relation to a set of standard dynamic parameters, st χ [START_REF] Gautier | Identification of Robot Dynamics[END_REF], [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF]:

( ) idm st st τ IDM q,q,q χ = & && (2)
Where ( ) st IDM q,q,q & && is the ( )

x s n N jacobian matrix of idm τ , with respect to the ( )

x1 s N
vector st χ of the This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 3 standard parameters given by: T T T T ...

χ τ   =   (3)
where:

• j j j j j j XX , XY , XZ , YY , YZ , ZZ are the six components of the inertia matrix, j j J , of link j at the origin of frame j,

• j j j MX , MY , MZ are the components of the first moments, j j MS , of link j,

• j M is the mass of link j,

• j
Ia is a total inertia moment for rotor and gears of actuator j.

• j Fv , j Fc are the viscous and Coulomb friction parameters of joint j. with respect to the sign of the velocity and tj Of is due to the current amplifier offset which supplies the motor.

• 14 s N * n =
is the number of standard parameters. The base parameters are the minimum number of dynamic parameters from which the dynamic model can be calculated. They are obtained from the standard inertial parameters by eliminating those which have no effect on the dynamic model, and by regrouping some others by means of linear relations. They can be determined using simple closed-form rules [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF] or a numerical method based on the QR decomposition [START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF].

The minimal inverse dynamic model can be written as:

( ) idm τ IDM q,q,q χ = & && (4)
Where:

( ) IDM q,q,q & && is the ( )

x n b matrix of the minimal set of basis functions of the rigid body dynamics, [START_REF] Ha | An efficient estimation algorithm for the model parameters of robotic manipulators[END_REF] χ is the ( )

x b 1 vector of the b base parameters. ( 6
)
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 q,q,q χ e τ = + = + & && (7)
Equation ( 7) represents the Inverse Dynamic Identification Model (IDIM).

We consider the off-line identification of the base dynamic parameters χ , given measured or estimated off-line data for τ and ( ) q, q, q & && , collected while the robot is tracking some planned This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 4 trajectories.

Usually, the signals available from the robot controller are the joint position measurement and the

( )

x n 1 control signal vector v τ , calculated according to the control law.

Then ( )

q, q, q
& && in [START_REF] Prüfer | Identification of Robot Dynamics with Differential and Integral Models : a Comparison[END_REF] are estimated with ( ) ˆq, q, q & && , respectively, obtained by bandpass filtering the measure of q [START_REF] Gautier | Dynamic Identification of Robots with Power Model[END_REF]. The derivatives are calculated off-line without phase shift, using a central difference algorithm of the lowpass filtered position q . The filtered position q is calculated off-line with a non causal zero-phase digital filter by processing the input data, q , through a lowpass Butterworth filter in both the forward and reverse direction using the filtfilt procedure from Matlab.

The control signal, v τ , is connected to the input current reference of the current closed-loop of the amplifiers which supplies the motors. Assuming that the current closed-loop has a bandwidth greater than 500Hz, then its transfer function is equal to its static gain, c K , in the frequency range (less than 10Hz) of the rigid robot dynamics. Then, the actual force/torque, τ , is calculated with the relation:

g v τ τ τ = (8)
where:

g τ , is the ( ) xn n
diagonal matrix of the drive gains, with:

r c g K K K τ τ = (9)
where: Those parameters have a priori values, given by manufacturers, which can be checked with special tests [START_REF] Restrepo | Calibration of drive chain of robot joints[END_REF].

The inverse dynamic identification model (IDIM) ( 7) is calculated at a frequency measurement m f , using samples of ( )

ˆq, q, q
& && to calculate ( ) ˆÎDM q,q,q & && and samples of v τ to calculate τ with (8), at

different times k t , m k 1,...,n =
, while the robot is tracking a reference trajectory ( ) r r r q ,q ,q & && , during the time length obs T , of the trajectory.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 5

The equations of each joint are regrouped together on all the trajectory to get an over-determined linear system such that:

( ) ( )

fm fm fm ˆŶ τ W q,q,q χ ρ = +
& && [START_REF] Swevers | Optimal Robot Excitation and Identification[END_REF] With:

( )

m 1 fm j 1 j fm fm n fm j n Y ( t ) Y τ ... , Y ... Y ( t ) τ τ         = =             (11)
( )

( ( ) ( ) ( )) ( ( ) ( ) ( ) m m m 1 j fm 1 1 1 j fm fm n j fm n n n
ˆŴ IDM q t ,q t ,q t ˆŴ q,q,q ... , W ... ˆŴ ÎDM q t ,q t ,q t

        = =             & && & && & && (12)
where:

( ( ) ( ) ( )) j k k k ˆÎDM q t ,q t ,q t & &&
is the jth row of the ( )

x n b matrix of the basis functions, fm fm ˆˆˆŴ IDM q,q,q W q,q,q = & && & && , will be used to recall that fm W , is calculated with a sampling of ( ) ˆÎDM q,q,q & && .

( ( ) ( ) ( )) k k k ˆÎDM q t ,q t ,q t & && , (5)
In order to eliminate high frequency force/torque ripple in τ , and to window the identification

m d f n .
After the data acquisition procedure and the parallel decimation of [START_REF] Swevers | Optimal Robot Excitation and Identification[END_REF], we obtain the overdetermined linear system:

( ) ( ) ˆŶ τ W q,q,q χ ρ = + & && (13)
where:

• ( )

Y τ is the x1 ( r ) vector of measurements, built from the actual force/torque τ ,

• ( ) ˆŴ q,q,q & && is the x ( r b ) observation matrix, built from the estimated values ( )

ˆq,q,q
& && of ( )

q, q, q & && .
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 6

• ρ is the x1 ( r ) vector of errors.

• m d r n* n / n = is the number of rows in [START_REF] Khalil | Modeling identification and control of robots[END_REF]. In Y and W , the equations of each joint are grouped together such that: ,

1 1 n n Y W Y ... W ... Y W         = =             (14)
where j Y and j W represent the m d n / n equations of joint j .

The ordinary LS (OLS) solution χ minimizes the squared 2-norm 2 ρ of the vector of errors.

Using the base parameters and tracking "exciting" reference trajectories [START_REF] Gautier | Exciting trajectories for the identification of the inertial parameters of robots[END_REF], we get a full rank and well conditioned matrix W . The LS solution χ is given by: ( ) ( )

1 T T χ W W W Y W Y - + = = (15)
It is computed using the QR factorization of W . Standard deviations i χ σ , are estimated using classical results from statistics under the assumptions that W is a deterministic matrix, according to the data filtering procedure described above, and ρ , is a zero-mean additive independent Gaussian noise, with a covariance matrix ρρ C , such that:

T 2 ρρ ρ r ( ) σ C E ρρ I = = ( 16
)
where E is the expectation operator and r I , the ( )

x r r identity matrix.

An unbiased estimation of the standard deviation ρ σ is:

2 2 ρ σ (r b) ˆŶ Wχ = - - (17)
The covariance matrix of the estimation error is given by:

T 2 T 1 χχ ρ [( )( ) ] σ ( ) ˆˆĈ E χ χ χ χ W W - = - - = (18)
i 2 χ χχ σ C ( ) ˆˆˆi ,i = is the i th diagonal coefficient of χχ ˆĈ . The relative standard deviation ri χ %σ ˆ is given by: ri i χ χ i %σ 100σ χ ˆˆ= , for i χ ˆ≠ 0 ( 19
)
The OLS can be improved by taking into account different standard deviations on joint j equations errors [START_REF] Gautier | Dynamic Identification of Robots with Power Model[END_REF]. Each equation of joint j in (13), [START_REF] Vandanjon | Comparison of two identification techniques: Theory and application[END_REF], is weighted with the inverse of the standard deviation of the error calculated from OLS solution of the equations of joint j , given by: ( ) ( ) ( )

j j j j j ˆŶ τ W IDM q,q,q χ ρ = + & && (20)
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 7

This weighting operation normalises the errors in [START_REF] Khalil | Modeling identification and control of robots[END_REF] and gives the weighted LS (WLS) estimation of the parameters.

This identification method is illustrated in Fig. 1.

Compared with the OE method described in the following section III, the use of IDIM, which is an analytical function of ( ) q,q,q & && , is particularly interesting because it does not require the integration of the direct dynamic model [START_REF] Featherstone | Dynamics[END_REF]. Moreover, χ is a one step linear LS solution which does not need initial conditions. However, the calculation of the velocities and accelerations are required using well-tuned bandpass filtering of the joint position [START_REF] Gautier | Dynamic Identification of Robots with Power Model[END_REF].

Robot ( )

In v erse D yn am ic

Id en tificatio n M o d el ID IM ˆÎD M q ,q ,q & && ˆ, , & && q q q ( ) ( ) fm ˆŴ IDM q ,q ,q & && Linear LS 2 ˆm in - Y W χ χ χ = ( ) ( ) ˆ( ), , , Y W IDM q q q τ & && ( ) q t ( ) t τ Control law χ ( ) r r r q ,q ,q & &&

III. THE OUTPUT ERROR METHOD (OE)

The OE identification methods minimize a quadratic error between an actual output y , and a simulated output s y , of the system, assuming both the actual and the simulated systems have the same input. This approach can be implemented in an open-loop form, [START_REF] Walter | Identification of parametric models from experimental data[END_REF], [START_REF]Identification of Continuous-time Models from Sampled Data[END_REF], or in a closed-loop form, [START_REF] Carrillo | Output error identification algorithms for continuous-time systems operating in closed-loop[END_REF], [START_REF] Landau | Closed-Loop Output Error Identification Algorithms for Nonlinear Plants[END_REF]. Considering a closed-loop controlled robot, the input, in the open loop scheme shown in Fig. 2, is the actual force/torque τ , and the input, in the closed-loop scheme shown in Fig. 3, is the reference trajectory ( ) r r r q ,q ,q & && . Because the open loop simulation of unstable robotic systems is very sensitive to the initial state conditions and to the errors in numerical algorithms which solve the differential equations, it is more suitable to choose the closed-loop form.

In both cases, the output is given by a state-space model output equation. Considering a robot and This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 8

taking the measured joint position as the output, the actual output is, y q = , and the simulated output is, s ddm y q = , as shown in Fig. 2 and Fig. 3, where ddm q ( t ) , is the solution of the differential equation given by the Direct Dynamic Model (DDM).

The DDM can be obtained by writing the IDM equation (1), as following:

( , ) = -( , , )
ddm ddm ddm ddm ddm M q q τ N q q χ χ && & (21)
where:

( , )

ddm M q χ and ( , , )

ddm ddm N q q χ &
depend on an estimation of the base parameters χ , ddm τ , is the force/torque input of the DDM.

The function ( ) ddm q t ,χ , is the result of the integration of the linear implicit differential equation [START_REF] Featherstone | Dynamics[END_REF] which can be written as a non-linear state-space model:

s s s s G( x )x f ( x ,u ) = & (22)
where:

ddm s ddm q x q   =     & , is the ( ) 2 x1
* n state-space vector,

s ddm u τ = , is the ( x1) n control input, x x G ( , ) -( , , ) n n n ddm s s s n n ddm ddm ddm ddm I 0 q ( x ) , f ( x ,u ) 0 M q τ N q q χ χ     = =         & & (23)
where, x n n 0 , is a ( x ) n n , matrix of zeros.

The linear output equation is given by:

s s s s s y C x D u = + (24)
Taking the measure of joint position as the output, s ddm y q = , we get:

[ ] x2 s n n *n C I 0 = , is the, ( x2* ) n n , output matrix, ( 25
) x s n n D 0 = , is the, ( x ) n n , direct feedthrough matrix. ( 26
)
Hence, for robotic systems, an OE identification method is based on the integration of the Direct Dynamic Model.

The optimal solution χ , minimizes the quadratic criterion ( ) J χ , given by:

( ) ( ) ( ) T 2 s s s J Y Y Y Y Y Y χ = - = - - (27)
where:

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 9

Y and s Y , are vectors obtained by filtering the vectors of samples fm Y and Sfm Y , respectively, where the equations of each joint are grouped together, with:

m 1 fm j 1 j fm fm n fm j n Y q ( t ) Y ... , Y ... Y q ( t )         = =         =     , 1 Sfm ddm 1 j Sfm Sfm n Sfm ddm k Y q ( t ) Y ... , Y ... Y q ( t )         = =             (28)
The minimization of ( ) 27), is a non-linear least-squares problem. The estimation of the parameters can be computed using algorithms such as the gradient method, the Newton methods or the Levenberg Marquardt method. These methods are based on a first or second order Taylor's expansion of ( ) ,( ) ,( )

J χ , (
J χ . Actual Robot Direct Dynamic Model (DDM) ddm ddm ddm ddm ddm M( q ) q N( q , q ) τ = - && & Non Linear LS 1 min 2 2 s ˆY Y χ χ = - ( ) t τ Control law r r r q q q           & && χ Sampling and filtering ( )
s ddm Y q Y q ( ) s ddm y q t = obs T ( ) y q t = ( ) t τ Fig. 2. Open-loop OE identification scheme. Actual Robot Direct Dynamic Model (DDM) ddm ddm ddm ddm ddm M( q ) q N( q , q ) τ = - && & Non Linear LS 1 min 2 2 s ˆY Y χ χ = - ( ) t τ Control law r r r q q q           & && χ Sampling and filtering ( )
s ddm Y q Y q ( ) s ddm y q t = obs T ddm τ ( ) y q t = Control law
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 10 Fig. 3. Closed-loop OE identification scheme.

In [START_REF] Gautier | DIDIM: A new method for the dynamic identification of robots from only torque data[END_REF], we used the Gauss-Newton method to calculate the optimal solution. It is a Newton method where approximations of the gradient and the hessian of ( ) J χ are calculated with the jacobian matrix of s y with respect to χ . The Gauss-Newton regression is a simpler way to calculate the optimal solution [START_REF] Davidson | Estimation and Inference in Econometrics[END_REF]. It is based on a Taylor series expansion of s y , at a current estimate k χ , of the parameters at

iteration k : ( ) ( ) +1 +1 k k k k k S S S χ y ( χ ) ˆŷ ( χ ) y ( χ ) χ χ o χ  ∂  = + - +   ∂   (29)
where:

( )

s y / k S χ χ y ( χ ) δ χ  ∂  =   ∂   (30)
s y /χ δ is the ( )

x n b , jacobian matrix of s y , with respect to χ , evaluated at k χ .

Each coefficient of s y /χ δ , defines a sensitivity function.

These sensitivity functions characterize the variation of the output function S y , with respect to a variation of the parameter χ . The sensitivity functions are the solutions of a differential system calculated from [START_REF] Featherstone | Dynamics[END_REF]. However, this technique is more time-consuming compared to the IDIM method.

Indeed, the DDM and the sensitivity functions must be integrated many times at each step of the iterative non-linear optimization method. Moreover, it is necessary to have good initial conditions in order to avoid multiple and local solutions.

Let us define:

+1 y y k S ( χ ) e = + (31)
From ( 29), it becomes:

( ) ( ) +1 / s k k k s y χ ˆŷ y ( χ ) δ χ χ o e - = - + + (32)
An over-determined linear system is obtained by filtering and sampling (32) over the time window obs T :

+1 k δ Y W χ ∆ ∆ ρ = + (33) with:
( )

+1 +1 k k k χ χ χ ∆ = -
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 11 Y ∆ , δ W , and ρ are, respectively, the sampling and filtering of ( )

k s ŷ y ( χ ) - , / s y χ
δ , and of ( )

o e + . +1 k χ ∆ is the LS solution of (33). This process is iterated with a new estimate, +1 +1 k k k ˆˆχ χ χ ∆ = + , until: k 1 k k 1 tol ρ ρ ρ + - ≤ , and, +1 k k i i 2 k i 1,...,b i ˆχ χ max tol χ = - ≤ (34)
where, 1 tol and 2 tol , are values ideally chosen to be small numbers to get fast convergence with good accuracy.

IV. DIDIM: DIRECT AND INVERSE DYNAMIC IDENTIFICATION MODEL TECHNIQUE

A. Theoretical approach

In the OE method as shown in Fig. 3, the actual output is the measured joint position, y q = .

We propose to change the output, y , from the actual joint position q , to the actual joint force/torque τ , and the simulated output s y , from the simulated joint position, ddm q , to the simulated joint force/torque, ddm τ . Then, we take y τ = , and s ddm y τ = , according to Fig. 4.

Actual Robot

Direct Dynamic Model ddm ddm ddm ddm ddm M( q , ) q N( q , q , ) This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 12

χ τ χ = - && & Non Linear LS min 2 s ˆY Y χ χ = - ( ) t τ Control law r r r q q q           & &&
The optimal solution, χ , minimizes the quadratic criterion, ( ) ( )

J χ , (
m 1 fm j 1 j fm fm n fm j n Y ( t ) Y τ ... , Y ... Y ( t ) τ τ         = =             , j j m 1 ddm 1 Sfm j Sfm Sfm n Sfm ddm n ( t ) Y Y ... , Y ... Y ( t ) τ τ         = =             (35)
This non-linear LS problem is solved by the Gauss-Newton regression as explained in section III.

The input force/torque of the DDM, ddm τ , can be calculated with the analytical expression of the inverse dynamic model ( 4), such as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) s ddm idm ddm ddm ddm y χ τ χ τ χ IDM q χ ,q χ ,q χ χ = = = & && (36)
The Taylor series expansion [START_REF] Davidson | Estimation and Inference in Econometrics[END_REF], with s ddm y τ = , at a current estimate, k χ , of the parameters χ , at iteration k , is calculated with the jacobian matrix of ( ) ddm τ χ , given by: (

/ k k k k k k ddm idm χ ddm ddm ddm ˆχ χ ˆˆˆδ IDM q ( χ ),q ( χ ),q ( χ ) χ χ χ χ τ τ   ∂ ∂   ∂ = = =     ∂ ∂ ∂     & && ) ( ) ( ) ( ) s y
Then, it becomes:

( ) ( ) ( ) ( ) ( ) k k k k k k k ddm ddm ddm ddm ddm ddm k k k k ddm ddm ddm ˆˆˆˆˆˆÎDM q ( χ ),q ( χ ),q ( χ ) χ IDM q ( χ ),q ( χ ),q ( χ ) ... χ ˆˆˆÎDM q ( χ ),q ( χ ),q ( χ ) χ χ ∂ = + ∂ ∂ ∂ & && & && & && (38)
The calculation of the second term on the right side of (38) needs to calculate the expression:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
k k k k k k ddm ddm ddm ddm ddm ddm ddm ddm k k k ddm ddm ddm ddm ddm k k k ddm ddm ddm ddm ddm q ˆˆˆˆˆÎDM q ( χ ),q ( χ ),q ( χ ) IDM q ( χ ),q ( χ ),q ( χ ) ... χ q χ q ˆˆÎDM q ( χ ),q ( χ ),q ( χ ) ... q χ q ˆˆÎDM q ( χ ),q ( χ ),q ( χ ) q χ ∂ ∂ ∂ = + ∂ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ∂ ∂ (39)
Let us recall that the joint force/torque y τ = , is obtained while the robot is tracking a reference trajectory, ( ) r r r q ,q ,q & && , with a closed-loop control law. The closed-loop simulation uses the direct dynamic model, the same control law and the same reference trajectory ( ) r r r q ,q ,q & && , as the actual one, to calculate S y .

In the following section IV.B, we show how to tune the control law of the closed-loop simulation in This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 13 order to keep the same bandwidth and stability margin as the actual closed-loop for any k χ , obtained at iteration k. This assumes for the simulated tracking error to keep close to the actual one for any k χ , that is to say:

( ) ( ) k k k ddm ddm ddm ˆˆq ( χ ),q ( χ ),q ( χ ) q,q,q & && & && , for any k χ (40)
This means that ( )

ddm ddm ddm q ( χ ),q ( χ ),q ( χ ) & &&
, have little dependence on χ , such that:

ddm ddm ddm q q q 0 χ χ χ ∂ ∂ ∂ ∂ ∂ ∂ & &&
Then (39) is simplified as:

( ) ( ) 0 k k k ddm ddm ddm ˆˆÎDM q ( χ ),q ( χ ),q ( χ ) χ ∂ ∂

& &&

Taking into account this simplification, we have in (38):

( ) ( ) ( ) k k k k k k ddm ddm ddm k ddm ddm ddm ˆˆˆˆˆˆÎDM q ( χ ),q ( χ ),q ( χ ) χ IDM q ( χ ),q ( χ ),q ( χ ) χ ∂ ∂ & && & &&
As a result, the jacobian matrix (37) can be approximated by: (

/ k k k k k k k χ ddm ddm ddm ddm ddm ddm ˆˆˆˆˆˆδ IDM q ( χ ),q ( χ ),q ( χ ) χ IDM q ( χ ),q ( χ ),q ( χ ) χ ∂ = ∂ & && & && ) ( ) ( ) s y
Each sensitivity function in the jacobian matrix is approximated by an algebraic equation. This is much more simpler than for usual OE method where the sensitivity functions are the solutions of complicated differential equations. This is the reason why it is much more simpler to minimize the error between the measured force/torque and the simulated force/torque than to minimize the error between the actual position and the simulated position.

Taking the approximation (41) of the jacobian matrix into the Taylor series expansion (32), it becomes:

( )( ) ( ) +1 k k k k k k s ddm ddm ddm ˆˆˆˆŷ y ( χ ) IDM q ( χ ),q ( χ ),q ( χ ) χ χ o e τ = = + - + + & && (42)
From (36), it becomes:

( )

k k k k k k s idm ddm ddm ddm ˆˆˆˆˆŷ ( χ ) ( χ ) IDM q ( χ ),q ( χ ),q ( χ ) χ τ = = & && (43)
Taking (43) in (42), it becomes:

( ) ( ) +1 k k k k ddm ddm ddm ˆˆŷ IDM q ( χ ),q ( χ ),q ( χ ) χ o e τ = = + + & && (44)
This is the Inverse Dynamic Identification Model, IDIM, [START_REF] Prüfer | Identification of Robot Dynamics with Differential and Integral Models : a Comparison[END_REF], where ( )

q, q, q
& && are estimated with ( )

ddm ddm ddm q ,q ,q & &&
, simulated with k DDM ( χ ) [START_REF] Featherstone | Dynamics[END_REF]. At each iteration k , the IDIM method is applied as This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 14

described in section II.

The sampling of (44) at a sampling rate m f , gives an over-determined linear system such as:

( ) ( ) , k fm fm ddm ddm ddm fm Ŷ τ W q ,q ,q χ ρ δ χ = + & && (45)
With:

( )

m 1 fm j 1 j fm fm n fm j n Y ( t ) Y τ ... , Y ... Y ( t ) τ τ         = =             (46) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) 1 1 1 , , , m m m j k 1 ddm ddm ddm fm k j fm ddm ddm ddm fm n j k fm ddm n ddm n ddm n ÎDM q t ,q t ,q t χ W Ŵ q ,q ,q χ ... , W ... W ÎDM q t ,q t ,q t χ δ δ δ δ           = =             & && & && & && (47)
The parallel decimation of (45) gives:

( ) ( ) , k ddm ddm ddm Ŷ τ W q ,q ,q χ χ ρ δ = + & && (48)
The LS solution of (48

) gives k 1 χ + , at iteration + k 1.
This process is iterated until:

k 1 k k tol 1 ρ ρ ρ + - ≤ , and, +1 k k i i 2 k i 1,...,b i ˆχ χ max tol χ = - ≤ ,
where, 1 tol and 2 tol are values ideally chosen to be small numbers to get fast convergence with good accuracy.

This new identification method is based on a closed-loop simulation using the direct dynamic model (DDM) while the optimal parameters minimize the 2-norm of the error between the actual force/torque τ , and the simulated force/torque ddm τ , over an observation window time obs T . This new technique overcomes the problems of non-linear optimization in OE method, section III, using the IDIM to calculate the simulated force/torque vector, s ddm idm y τ τ = = . Because this method uses both models DDM and IDIM, it is named the DIDIM method: Direct and Inverse Dynamic Identification Models technique.

The DIDIM method with the Gauss-Newton regression is illustrated Fig. 5.

This approach is particularly interesting thanks to the following reasons: • It needs only the actuator force/torque measurement or estimation,

• It avoids tuning the bandpass filter in the IDIM method by using the integration of the DDM in a closed-loop simulation where the tuning of the bandwidth automatically defines the same frequency This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 15 range for the dynamics of the actual and of the model to be identified. • It combines the inverse and the direct dynamic model and validates, in the same identification procedure, both models for computed torque control and for simulation.

• It dramatically simplifies the computation of the matrix of the sensitivity functions which is given by an algebraic equation (the inverse dynamic identification model) whereas it is given by the resolution of a complicated system of differential equations in the usual OE method.

The drawback is that the structure and the tuning of the actual closed-loop control law must be known to be implemented in the closed-loop simulation of the robot. Most often, this is not a real problem, because working on identification for simulation or control of the robot, needs a minimal knowledge on the robot controller.

r r r q q q           & && Actual Robot Direct Dynamic Model k ddm ddm k ddm ddm ddm M( q , ) q N( q , q , ) χ τ χ = - && & Linear LS +1 min 2 k ˆY W δ χ χ χ = - ( ) t τ Control law χ ( ) ( ) ( ), W , , , k ddm ddm ddm Y IDM q q q δ τ χ & && obs T ddm τ ( ) q t
Control law ( )

Inv erse D ynam ic Id entificatio n M o d el k d dm dd m dd m ÎD M q ,q ,q , χ & && ( ) k ddm ddm ddm ÎDM q ,q ,q ,χ & && ( ) t y τ = ddm q
Actual closed loop robot Simulated closed loop robot sampling ( ) lowpass filtering downsampling fm Fig. 5. DIDIM with the Gauss-Newton regression, identification scheme.

B. Initialization of the algorithm

A problem is how to choose the initial values 0 χ .

We can use CAD values, or identified values with the IDIM method, but we show that there is no need at all of a priori values.

We propose an algorithm not sensitive to the initial conditions, which assumes that the condition ( ) ( )

ddm k ddm k ddm k ˆˆq ( χ ),q ( χ ),q ( χ ) q,q,q & && & &&
, is satisfied at any iteration k , and especially for k =0. This is possible by taking the same control law structure for the actual robot and for the simulated one with the same performances given by the bandwidth, the stability margin or the closed-loop poles.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 16

Because the simulated robot parameters k χ , change at each iteration k , the gains of the simulated control law must be updated according to k χ .

For example, let us consider a PD control law for each joint j . The inverse dynamic model IDM [START_REF] Khosla | Parameter Identification of Robot Dynamics[END_REF] for the joint j , can be written as a decoupled double integrator perturbed by a coupling force/torque, such that:

= ( ) + ( , )= ( ) ( ) + ( , )= ( ) j n n j idm j ,i i j j , j j j ,i i j j , j j j i 1 i j
τ τ M q q N q q M q q M q q N q q M q q p

= ≠ = + - ∑ ∑ && & && && & && (49)
where j p is considered as a perturbation given by: ( ) ( , )

n j j ,i i j i j p M q q N q q ≠ = - - ∑ && & (50) ( ) j ,i
M q which depends on q , is approximated by a constant inertia moment j J , given by: ( )

( ) j j j j a j , j j a q J ZZ I max M q ZZ I = + + - - (51)
j J , is the maximum value, with respect to q , of the inertia moment around joint j z axis. This gives the smallest damping value and the smallest stability margin of the closed-loop second order transfer function (55), while q varies.

It can be calculated from a priori CAD values of inertial parameters and must be taken at least as j j a ZZ I + .

The joint j dynamic model is approximated by a double integrator, where j p , is a perturbation, as following:

( ) ( ) ( ) j j j j j j , j j 1 1 q τ p τ p M q J = + + && (52)
Let us consider the joint j PD control of the actual robot which is illustrated Fig. 6:

+ - + - j a g τ j r q j a v k 1 a j J 1 s 1 s + + j a p k j p j v τ j q & j q j τ j q &&
Fig. 6. Joint PD control of the actual robot.

The control input calculated by the robot controller is given by: This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 17 ( )

j j j j j a a a p v r j v j v k k q q k q τ = - - & (53) j v τ
is the current reference of the current amplifiers which supplies the motor.

The joint j , force/torque is given by:

j j a j g v τ τ τ = (54)
where: j a g τ is the actual drive gain, calculated with the actual parameters in [START_REF] Gautier | Dynamic Identification of Robots with Power Model[END_REF].

a j
J is the actual value of j J .

In order to tune the tracking performances of the reference position j r q , the transfer function j r j q q is calculated with =0

j p : =0 j j j j j j j j a 2 a 2 j j r p a a a a a 2 a v p p nj nj q 1 1 H J s 2 q 1 s s 1 s 1 g k k k τ ζ ω ω    = = =      + + + + (55)
where:

a nj ω is the actual natural which characterizes the closed-loop bandwidth,

a j
ζ is the actual damping coefficient which characterizes the closed-loop stability margin, with:

j i i a a a a nj p v a j g k k J τ ω = , 1 j i i a a v a j a a p j g k 2 k J τ ζ = (56)
Then it becomes: The variables ( )

, j j j j j ddm ddm ddm ddm ddm v , q , q , q τ τ & &&
, in Fig. 7, are computed by numerical integration of k DDM ( χ ) , [START_REF] Featherstone | Dynamics[END_REF].

The control law of the simulated robot has the same structure as the actual one, Fig. 6, where we take:

ddm k ddm k ddm k ˆˆq ( χ ),q ( χ ),q ( χ ) q,q,q & && & && , at each iteration k.
We propose to take a regular inertia matrix 0 ( , ) ddm M q χ , in order to have a good initialization for the numerical integration of the DDM ( 21) . This is named the "regular initialization".

It can be obtained with:

0 0 χ = , except for, 0 1, j Ia j 1,n = = (61)
The inertia of the rotor and gear of actuator j is generally taken into account in the IDM model ( 1) as:

τ j r j j
Ia q = && Then, the initial inertia matrix becomes the identity matrix, which is the best regular matrix:

0 ( , ) = ddm n M q I χ (62)
Another simple regular initialization is to take:

0 0 χ = , except for, 0 1, j ZZ j 1,n = = ( 63
)
The initial inertia matrix, 0 ( , ) ddm M q χ , is no more the identity matrix, but remains regular.

Another point is to choose the state initial condition of the state vector, ( )

(0) (0) ddm ddm q ,q &
, in order to integrate the DDM [START_REF] Featherstone | Dynamics[END_REF]. Because DIDIM doesn't need the joint position measurement, the actual values ( )

(0) (0) q ,q &
, are supposed to be unknown and we choose, ( ) ( )

(0) (0) (0) (0) ddm ddm r r q ,q q ,q = & & , which is close to ( ) (0) (0) q ,q & . Because

V. CASE STUDY: MODELING OF THE SCARA ROBOT

The identification method is carried out on a 2 degree-of-freedom planar direct drive prototype robot This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 20

without gravity effect, shown in Fig. 8. This direct drive prototype is very suitable for the study of DIDIM because it emphasizes non linear coupling torques while this non-linear effect is divided by at least 2500 for industrial robots with gear ratio greater than 50. Moreover, the dynamic model of this robot depends on eight parameters only, which facilitates the study of the identification efficiency with respect to several conditions. At last, this robot and its real parameters, called the nominal parameters, are well known. Thus, we can check the physical meaning of the identified parameters.

The description of the geometry of the robot uses the modified Denavit and Hartenberg (DHM) notations [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF] which are illustrated in Fig. 9. The robot is direct driven by 2 DC permanent magnet motors supplied by PWM amplifiers.

[ ] T 1R 1 1 2 R 2 2 2 2 ZZ Fv Fc ZZ LMX LMY Fv Fc χ = (64) 2 1R 1 1 2 ZZ ZZ Ia M L = + + 2 R 2 2 ZZ ZZ Ia = + L =0.5m
, is the length of the first link.

In the case of the SCARA robot, the parameters,

:,k IDM k = , of ( )
IDM q,q,q & && , in IDIM [START_REF] Prüfer | Identification of Robot Dynamics with Differential and Integral Models : a Comparison[END_REF], are the following:

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 21

0 0 0 1 R 1 1 2 R 2 1 1 :,1 ZZ :,2 Fv 1 2 1 :,3 Fc :,4 ZZ 1 2 ) 1 2 2 2 1 2 2 :,5 LMX 2 1 2 1 q q IDM IDM , IDM IDM ,
q q sign( q ) IDM IDM , IDM IDM , q q ( 2q q ) cos q -q ( 2q q sin q IDM IDM q cos q q sin

+     = = = =         +     = = = =     +     + = = + && & && && & && && && && & & & && & 0 0 2 2 2 2 ) 1 2 2 2 1 2 2 :,6 LMY 2 1 2 1 2 :,7 Fv :,8 Fc 2 2
, q ( 2q q ) sin q q ( 2q q cos q IDM IDM , q cos q q sin q IDM IDM , IDM IDM q sign( q )

+       - + -   = =   -       = = = =         (65)
The closed-loop control is a PD control law (53) , according to Fig. 6, with:

2 1 1R 2 R 2 J ZZ ZZ LMX = + + , and
2 2 R J ZZ = .
The actual gains are calculated with (58), taking a desired damping,

VI. EXPERIMENTAL IDENTIFICATION RESULTS

The new identification process is performed in different cases in order to compare the previous IDIM technique to the new DIDIM technique and to investigate the robustness of DIDIM with respect to the initialization, to the acquisition sampling rate, to the data filtering and to the closed-loop tuning.

All the results are given in SI units, on the joint side.

A. Comparison of IDIM and DIDIM with good initial values,

0 IDIM ˆˆ χ χ = .
At first, the algorithm is initialized with, IDIM χ , the vector of parameters identified with the IDIM LS estimator.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 22

The IDIM LS off-line estimation is carried out with a filtered position q , calculated with a 20Hz cutoff frequency forward and reverse Butterworth filter, and with the velocities q & , and the accelerations, q && , calculated with a central difference algorithm of q . The parallel decimation of fm Y and fm W , in [START_REF] Swevers | Optimal Robot Excitation and Identification[END_REF], is carried out with a sample rate divided by a factor, d n =20, and a lowpass filter cut-off frequency equal to, 0.8* /(2* )=4Hz

m d f n .
The results are given in Table 1. It needs only 2 steps to obtain the optimal solution which is very close to the IDIM solution. Hence, the DIDIM method does not improve the IDIM solution calculated with good bandpass filtered data.

A validation is plotted on Fig.

ˆŶ W q ,q ,q δ χ χ = & &&
, as defined in (45), ( 46), (47), are very close.

Both methods give a small relative norm error, Ŷ W / Y χ -<3%, which shows a good accuracy for the model and for the identified values.

It can be seen that the parameters,

1 Fv , and 2
Fv , have no significant estimations because of their large relative standard deviation (>30%). They have no significant contribution in the joint torques and they can be cancelled to keep a set of essential parameters of a simplified dynamic model, without loss of accuracy [START_REF] Khalil | Identification of the payload inertial parameters of industrial manipulators[END_REF].

However, we prefer to keep all the parameters in the following, for a better comparison of IDIM and DIDIM identification methods.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 23 The initial values of the dynamic parameters are given by (61), with:

[ ] T 1 0 0 1 0 0 0 0 0 ˆ χ =
The identified values given in Table 2, are very close to those given in Table 1. This result validates the regular initialization procedure , described in section IV.B.

Moreover the algorithm converges in only 3 steps and is not time consuming.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 24 The relative norm errors on joint position, velocity and acceleration are plotted in Fig. 11 ,with the following legend:

• norm error relative to the actual filtered joint position, j ddm j j ˆq q / q -, velocity j ddm j j ˆq q / q -, and acceleration, j ddm j j ˆq q / q -&& && && , where ( )

ˆq,q,q
& && , are calculated as given in section VI.A.

* norm error relative to the reference joint position, j j j ddm r r q q / q -, velocity,

j j j ddm r r q q / q - & & &
, and acceleration, j j j ddm r r q q / q -&& && && .

The assumption (40), made in section IV.B, (

) ( )

ddm k ddm k ddm k ˆq ( χ ),q ( χ ),q ( χ ) q,q,q & &&
& && , at each iteration k, is confirmed on Fig. 11. , with a constant relative norm error close to 0.5% for the position, 5%, for the velocity and 10%, for the acceleration.

These results validate the updating procedure (59), of the simulated PD control law gains.

It can be seen also on Fig. 11. , that the simulated trajectory, ( )

ddm k ddm k ddm k ˆˆq ( χ ),q ( χ ),q ( χ ) & &&
, is 3 to 5 times closer to the actual one, ( ) ˆq,q,q & && , than to the reference one, ( ) r r r q ,q ,q & && , with a relative norm error close to 1.5% for the position, 15%, for the velocity and 30%, for the acceleration. Moreover, this error depends on the closed-loop bandwidth. It means that computing the observation matrix in [START_REF] Khalil | Modeling identification and control of robots[END_REF] with the reference trajectory, ( ) r r r q ,q ,q & && , leads to a bad identification of the dynamic parameters.

Then, the right assumption is, (

) ( ) We have seen that ( ) ( )

ddm k ddm k ddm k ˆq ( χ ),q ( χ ),q ( χ ) q,q,q & && & && , ( 40
ddm k ddm k ddm k ˆq ( χ ),q ( χ ),q ( χ ) q,q,q & &&
& && , at each iteration k, with a constant small error. On the contrary, the relative torque norm error, given in Table 3, and plotted on Fig. 12, dramatically decreases in only 3 steps. This shows the fast algorithm convergence. The fast convergence of each parameter is shown in Table 4 , and is plotted on Fig. 13. The actual torque and the simulated data are resampled to obtain a low frequency measurement m f = 0.5Hz. This is a downsample procedure without lowpass antialiasing filtering which investigates a real problem on industrial robots where the available sample rate measurement given by the controller may be much lower than the control sample rate. All the actual and simulated data are sampled at m f = 0.5Hz.

The IDIM LS estimation is carried out with the measured joint position q , and with ( )

Ŷ W / Y χ - =0.5 3 Ŷ W / Y χ - =0.04
The identified values with IDIM are not good while the identified values with DIDIM are still good.

This shows the robustness of DIDIM with respect to the sampling rate measurement.

IDIM fails because there is an amplitude distortion in the estimation of ( )

ˆq, q
& && , with a central difference of q , sampled at a too low frequency m f . This point is illustrated in Table 6, which gives the relative norm errors on velocity (80%) and acceleration (80%).

( )

(200Hz) (200Hz) ˆq , q & &&
, is calculated from q , sampled at 200Hz and lowpass filtered at a 0.5Hz cutoff frequency, and derived with a central difference algorithm.

( )

(0.5Hz) (0.5Hz) ˆq , q & &&
, is calculated from q , sampled at 0.5Hz and derived with a central difference algorithm.

DIDIM succeeds because, ( )

ddm ddm ddm q ,q ,q & &&
, is computed with accuracy by the integration of the DDM with a well-tuned variable step solver, and it can be sampled without error at any frequency m f .

f = (200Hz) (0.5Hz) (200Hz) 1 1 1 ˆˆq q / q - 0.39 (200Hz) (0.5Hz) (200Hz) 1 1 1 ˆˆq q / q - && && && 0.73 (200Hz) (0.5Hz) (200Hz) 2 2 2 ˆˆq q / q 0.80 (200Hz) (0.5Hz) (200Hz) 2 2 2 ˆˆq q / q - && && && 0.81
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 28

D. Comparison of IDIM and DIDIM, without data filtering.

All the actual and simulated data are sampled at m f = 200Hz.

The IDIM LS estimation is carried out with the measured joint position q , and with ( )

ˆq,q
& && , calculated by a central difference algorithm of q , without lowpass Butterworth filtering. There is no parallel decimation. DIDIM starts with the regular initialization. Results are given in Table 7.

Ŷ W / Y χ - =0.8 2 Ŷ W / Y χ - =0.08
The identified values with IDIM are not good while the identified values with DIDIM are still good.

IDIM fails because of the too large noise in the observation matrix, ( ) fm ˆŴ q,q,q & && , coming from the derivation of q , without lowpass filtering. Then the LS estimation is biased.

DIDIM succeeds because the observation matrix, ( )

, k fm ddm ddm ddm Ŵ q ,q ,q δ χ & &&
, is calculated without noise with the simulated values ( )

ddm ddm ddm q ,q ,q & && .
This validation shows that DIDIM cancels the bias of IDIM estimation, coming from a noisy estimation of ( ) ˆq,q,q & && , which gives a too noisy observation matrix ( ) fm ˆŴ q,q,q & && . The DIDIM identification is performed taking half the values of the full ones given in section V, /2=1/2 (rd/s) The relative norm errors on joint position, velocity and acceleration are plotted on Fig. 14, with the same legend as those given for Fig. 11, section VI.B.

It can be seen that, (

) ( ) ddm k ddm k ddm k ˆq ( χ ),q ( χ ),q ( χ ) q,q,q & && & && , at each iteration k, with a constant norm error larger but close to the value obtained with the full bandwidth, Fig. 11, close to, 0.5% for the position, 5%, for the velocity and 10%, for the acceleration. This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 30

The relative torque norm error, given in Table 9, and plotted in Fig. 15, decreases in 6 steps, only twice more than with the full bandwidth, Table 3, Fig. 12. This shows that DIDIM is not very sensitive to error in the simulated closed-loop bandwidth, provided the control law structure is known.

However, DIDIM fails beyond 1/3 of the full bandwidth, with 3 Output Error approach, but considering the output is no more the joint position but the joint force/torque.

The optimal parameters are the solution of a non-linear least-squares problem which is solved with a Gauss-Newton method. Each step of the iterative procedure of the Gauss-Newton regression is dramatically simplified to a linear regression which is solved with the Inverse Dynamic Identification Model technique (IDIM). Then, DIDIM mixes the closed-loop OE technique and the IDIM technique.

DIDIM needs a closed-loop simulation of the robot using the direct dynamic model (DDM) and assuming the same structure of the control law and the same reference trajectory for both the actual and the simulated robot. Then, it needs to initialize the parameters and the state vector of the DDM.

The difficulties for the choice of the initial conditions for non-linear LS problem are overcome with a "regular initialization" of the parameters and an updating of the control law gains at each step of the iterative procedure. The initial state is given by the initial values of the reference trajectory.

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 31

An experimental validation is carried out on a 2 dof robot. The following points were checked: • DIDIM gives the same results as IDIM, provided well-tuned data filtering for IDIM, adapted to the system dynamics, • DIDIM is robust to the initialization of both parameters and state,

• DIDIM is robust to the closed-loop performances tuning errors between the simulated and the actual closed-loop robot, provided the same control law structure.

Compared to IDIM, DIDIM technique is particularly attractive thanks to the following reasons:

• It needs only the actuator force/torque measurement or estimation,

• It avoids tuning the bandpass filter in the IDIM method by using the integration of the DDM in a closed-loop simulation where the tuning of the bandwidth automatically defines the same frequency range for the dynamics of the actual system and of the model to be identified, • It cancels bias in IDIM due to errors in bandpass filtering data, or no filtering at all, or too low frequency measurement, • It combines the inverse and the direct dynamic model and validates, in the same identification procedure, both models for computed torque control and for simulation. Up to now, the DDM was validated a posteriori in simulation.

  parameter where Fsj Of is the dissymmetry of the Coulomb friction

rK,

  is the ( ) xn n gear ratios diagonal matrix of the joint drive chains ( the ( ) xn n static gains diagonal matrix of the current amplifiers, K τ , is the ( ) xn n diagonal matrix of the electromagnetic motor torque constants.

Fig. 1 .

 1 Fig. 1. IDIM LS identification scheme.

Fig. 4 .

 4 Fig. 4. DIDIM identification scheme.This means that the output equation (24) of the state-space model[START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF] reduces to a direct feedthrough equation such as, s

ζFig. 7 .

 7 Fig. 7. Joint PD control of the simulated robot.

k

  value of j J , (51), calculated with the estimation k χ , at iteration k. are the gains of the simulated control law. They are calculated in order to keep the same performances for the simulated closed-loop and for the actual closed-loop, that is to say to keep the same desired values, d nj ω and d j ζ , for the closed-loop poles. Then, it becomes: not depend at all on the parameters values, but the derivative gain in the simulator , each iteration k . It is important to note that only the gain in the simulated closed-loop, j s v k , is modified during the iterative procedure. The actual gain of the robot control law,

  the closed-loop transient response due to different initial conditions differs between the actual and the simulated signals during a transient period of approximately, 5 d n / ω , the corresponding joint force/torque samples are eliminated from the identification data in (45).

Fig. 8 .

 8 Fig. 8. The scara robot prototype. Fig. 9. DHM frames of the scara robot.The dynamic model depends on 8 minimal dynamic parameters, considering 4 friction parameters:

  d j ζ =1, for joint 1 and joint 2.The desired natural frequency, d nj ω , is chosen according to the driving capacity without saturation of the joint drive. For this robot we obtain a full bandwidth with, of the control and of the measurement are equal to, m f =200Hz.Torque data are obtained from (54), and from the current reference data v τ .The simulation of the robot is carried out with the same reference trajectory and with the same control law structure as the actual robot.The gains in the simulator are calculated with (59) and with the same values,

χ

  Fig. 10. DIDIM, validation,

Fig. 11 .

 11 Fig.11. • norm error relative to the filtered actual position, velocity, acceleration. * norm error relative to the reference position, velocity, acceleration.

TABLE 3 :Fig. 12 .-

 312 Fig. 12. DIDIM, convergence of the joint torque error, j

Fig. 13 .

 13 Fig. 13. DIDIM, parameters convergence

Eω

  . DIDIM robustness with respect to error in the simulated closed-loop tuning, d n This section investigates the effect of an error between the actual value, a n ω , and the simulated value d n ω , of the natural frequency which represents the closed-loop bandwidth.

Fig. 14 .

 14 Fig. 14. • norm error relative to the filtered actual position, velocity, acceleration. * norm error relative to the reference position, velocity, acceleration.

Fig. 15 .

 15 Fig. 15. DIDIM, convergence of the joint torque error, j

  27), where, Y , and s Y , are vectors obtained by filtering the vectors of samples, fm Y and Sfm Y , respectively, where the equations of each joint are grouped, with:

TABLE 1 :

 1 COMPARISON OF IDIM AND DIDIM METHODS IDIM DIDIM Parameter

				χ	IDIM		χ 2 σ ˆ	χ %σ ˆ	r	0 ˆˆ χ χ =	IDIM	ˆ χ	2	χ 2 σ ˆ	χ %σ ˆ	r
	ZZ 1R	3.44	0.034	0.50	3.44	3.45	0.036	0.52
	1 Fv		0.03	0.031	52. 0	0.03	0.04	0.032	40.0
	1 Fc		0.82		0.1	6.0		0.82	0.82	0.05	3.0
	ZZ	2		0.062	0.0006	0.51	0.062	0.061	0.0006	0.49
	LMX	2	0.121	0.0014	0.56	0.121	0.124	0.0013	0.52
	2 LMY	0.007	0.0007	5.0		0.007	0.007	0.0005	3.5
	2 Fv		0.013	0.006	23.0	0.013	0.014	0.0084	30.0
	2 Fc		0.137	0.006	2.30	0.137	0.133	0.0080	3.0
				-	W	IDIM Ŷ χ	/ Y	=0.024					-	W	2 Ŷ / Y χ	=0.021

TABLE 2 :

 2 DIDIM WITH THE REGULAR INITIALIZATION Parameter

				ˆ χ	0	ˆ χ	3	χ 2 σ ˆ	χ %σ ˆ	r
	ZZ 1R	1	3.45	0.014	0.2
	1 Fv		0	0.02	0.012	15
	1 Fc		0	0.85	0.016	1.0
	ZZ	2		1	0.061	0.0001	0.1
	LMX	2	0	0.124	0.0002	0.1
	2 LMY	0	0.007	0.0003	2.0
	2 Fv		0	0.01	0.003	10
	2 Fc		0	0.132	0.0008	0.3

TABLE 4 :

 4 PARAMETERS CONVERGENCE Parameters

					ˆ χ	0	ˆ χ	1	ˆ χ	2	ˆ χ	3
		ZZ 1R	1	3.46 3.45 3.45
		1 Fv		0	0.04 0.02 0.02
		1 Fc		0	0.86 0.85 0.85
		ZZ	2		1	0.06 0.061 0.061
		LMX	2	0 0.122 0.124 0.124
		2 LMY	0	0.05 0.07 0.07
		Fv	2		0 0.005 0.01 0.01
		2 Fc		0 0.130 0.132 0.132
	0 1	1		2		3

TABLE 5 :

 5 Copyright may be transferred without notice, after which this version may no longer be accessible p 27 decimation. DIDIM starts with the regular initialization. Results are given in Table 5 . IDIM AND DIDIM, LOW SAMPLING RATE, m

	ˆq,q

& && , calculated by a central difference algorithm of q , without lowpass Butterworth filtering. There is no parallel This work has been submitted to the IEEE for possible publication.

TABLE 6 :

 6 IDIM, JOINT DATA ERRORS AT 0.5Hz

	m

TABLE 7

 7

					: IDIM AND DIDIM ESTIMATION WITHOUT DATA FILTERING
						IDIM							DIDIM
	Parameter	χ	IDIM	χ 2 σ ˆ	χ %σ ˆ	r	χ	0	ˆ χ	2	χ 2 σ ˆ	χ %σ ˆ	r
	ZZ 1R	1.50	0.05	1.60	1.0	3.45	0.007	0.1
	1 Fv		0.095	0.15	80. 0	0	0.05	0.023	21
	1 Fc		0.55	0.26	23.3	0	0.81	0.004	0.24
	ZZ	2		0.14	0.018	6.7		1.0	0.061	0.0004	0.3
	LMX	2	0.63	0.035	2.7		0	0.124	0.0015	0.3
	2 LMY	0.1	0.023	11.8	0	0.008	0.0009	5.6
	2 Fv		0.001	0.143	700.0	0	0.023	0.0022	48
	2 Fc		0.19	0.244	68.40	0	0.13	0.0038	1.5
						IDIM

TABLE 8 :

 8 This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible p 29shown in Table2, that is to say a frequency measurement, m f =200Hz, and a parallel decimation with a factor, d n =20, and a lowpass filter cut-off frequency equal to 4Hz .The parameters, given inTable 8, converge in only 6 steps to values which are very close to those obtained in Table 2, with a full closed-loop bandwidth. DIDIM, WITH SIMULATED HALF FULL BANDWIDTH , /2

	d	ω	1 n	=	d	ω	1 f n	and	d	ω	2 n	=	d	ω	2 f n / 2	=	10/2 (rd/s)	, and the same procedure used to obtain results

& && & & && &

& & & & && & &

&& &

Future work concerns the validation of DIDIM on a 6 dof industrial robot.

VIII. REFERENCES


Christian Sallé

Frédérique Ferey

Philippe Grosseau

Bernard Guilhot

Philibert Iacconi

Mourad Benabdesselam

Gilbert Fantozzi

Detection of foreign phases in doped a-alumina powders by thermoluminescence

Keywords: Al 2 O 3, CaAl 12 O 19, calcium, composites, defects, thermoluminescence, ZrO 2

Introduction

The objective of this work is to show that the thermoluminescence (TL) allows the detection of foreign phases in the alumina matrix. The thermoluminescence (TL) is a technique used to study the fine structure of crystals by revealing the point defects. The purpose of this paper is to give two examples of foreign phases detected in alumina.

The TL of monoclinic zirconia has already been studied [1][2][3]. Six peaks were found at -170, -140, -85, -30, 10 and 110°C [3], (for a measurement at 30°C min -1 ) and one of them (-30°C) was attributed to the presence of tetragonal zirconia, which is present as foreign phase in the sample of monoclinic zirconia [3][4].

The Al 2 O 3 ⋅CaO system was studied a lot according to its importance in metallurgy, ceramics or cements. The low limit of solubility of calcium in alumina (according to the system binary diagram) leads to the presence of the hibonite phase in the alumina matrix [5]. The optical properties of alumina are then changed.

Experimental procedure

The mixture preparation α-Al 2 O 3 -ZrO 2 mixed powders

The γ-alumina used in this study is a highly pure (>99.99% Al 2 O 3 ) commercial powder (Baikowski A125) obtained by thermal decomposition of ammonium alum (NH 4 Al(SO 4 ) 2 ⋅12H 2 O). Then α-alumina is obtained by a thermal treatment of γ-alumina at 1400°C under controlled atmosphere (oxygen). The zirconia is a commercial powder (supplied by Rhône-Alpes Chimie, France) and its purity is higher than 99%. XRD analysis shows that this zirconia powder exists in the form of monoclinic phase (baddeleyite).

The mixtures of oxides, α-Al 2 O 3 and ZrO 2 , were wet ground in an organic solvent (acetone), with continuous stirring. Various mixtures are used: 1, 5, 10…mass% of ZrO 2 .5gofeach mixture were prepared.

Then half part of these mixtures undergoes a thermal treatment under controlled atmosphere (oxygen), in order to reveal and study, by the techniques of X-rays diffraction and thermoluminescence, the effect of doping by zirconia on the alumina crystalline structure.

This treatment is conducted by the procedure described in Fig. 1.

α-Al 2 O 3 :Ca samples

The α-alumina doped by calcium samples are prepared using a commercial ultra-pure γ-alumina powder (Baïkowski CR140) obtained by thermal decomposition of ammonium alum (NH 4 Al(SO 4 ) 2 ⋅12H 2 O). The doping is performed through the impregnation of the γ-alumina powder by an alcoholic solution of calcium chloride [6]. All the doping concentrations (in ppm) are calculated on the basis of a mole ratio of Ca/Al. A further thermal treatment allows the decomposition of the solvent and the salt, and the diffusion of the calcium inside the host material (drying at 100 and 600°C for 24 h). The doped γ-alumina is then transformed into doped α-alumina by a thermal treatment at 1450°C for 2 h under a pure gas flow (O 2 or Ar+2%H 2 ). The doping by Ca 2+ was performed for concentrations in calcium ranging from 0 to 10000 ppm.

The experimental techniques

The X-ray diffraction (XRD) experiments have been performed with a Siemens D5000 diffractometer using CuK α radiation (λ α K =0.15406 nm) in a θ/2θ configuration. A database on CD-ROM is used for searching and retrieving X-ray powder diffraction patterns.

The TL emission is detected by a photomultiplier Hamamatsu R562 with a spectral response between 160 and 850 nm. For semiconducting or insulating solid crystals, the valence and the conduction bands are separated by an energy gap (Fig. 2). When the crystal is submitted to ionising radiation (UV, X-ray, γ…), an electron-hole pair is formed in the valence band. The electron is free to travel through the conduction band and to fall into an electron trap. While, the hole travels through the valence band and falls into a hole trap.

After irradiation, the trapped charges may escape if the temperature of the crystal is raised and may recombine together on a centre, which de-excites by photon emission. If the heating rate of the sample is constant and if the emitted luminous intensity is measured using a photomultiplier tube (PMT) linked to a pico-ammeter, some peaks are observed, their temperature position being characteristic of the trap nature and their intensity depending on concentrations of both the traps and the recombination centres.

For the irradiation, we expose the samples one minute under an UV lamp (λ=253.7 nm), or five minutes under a W-target (45 kV, 2 mA), at -190°C, and then we keep the samples two minutes in the dark.

The samples are then heated at 60°C min -1 in the case of UV exposure, and at 30°C min -1 in the case of X-ray irradiation.

Experimental results

The case of tetragonal zirconia

Pure oxides (monoclinic zirconia and α-alumina) were studied by low temperature TL after UV irradiation, before and after the thermal treatment at 1400°C under oxygen. The studied α-alumina presents no peak after UV irradiation. The monoclinic zirconia, which is used in this study, presents 6 TL peaks after UV irradiation situated at -170, -145, -90, -35, 0 and 95°C, but the peak at -35°C disappears when the powder undergoes the thermal treatment.

In the case of alumina-zirconia mixtures (before thermal treatment), we observe the same peaks as for pure zirconia, except the one at -35°C which doesn't appear. On the other hand (after thermal treatment), the peak at -35°C appears on the TL glow curve (Fig. 3). Above 1170°C, the stable phase of zirconia is the tetragonal one [7]. The presence of tetragonal zirconia in the powder is detected by XRD, but only for mixtures realised with a zirconia amount higher or equal to 10 mass% and treated at 1400°C (Fig. 4): XRD being not sensitive enough to detect tetragonal zirconia for lower amount of zirconia in the mixtures. In Al 2 O 3 -ZrO 2 composites, the tetragonal phase is stable at room temperature only if the zirconia grain size lies between two critical values (from about 0.1 up to 0.5 µm) [8].

So, in order to verify that the peak at -35°C is due to the presence of the tetragonal phase, we have studied a commercial mixture (containing 10 vol% of m-ZrO 2 ), allowing the stabilisation of zirconia in tetragonal form after sintering. This mixture is prepared in the following manner: an alumina powder (purity>99.9% Al 2 O 3 ) is homogeneously mixed with a monoclinic zirconia (purity>99.9% ZrO 2 ) in distilled water, with a conventional mixing-milling technique. Subsequently, they are spray-dried. The final size of the grains in the mixture was measured between 0.33 and 0.40 µm.

This mixture was characterised by TL and XRD, before and after the thermal treatment (at 1400°C under oxygen). Before the thermal treatment, we observe by TL (Fig. 5) an intense peak at -145°C, a peak towards -50°C, and a peak at 0°C which is probably the characteristic peak of monoclinic zirconia at 0°C. After thermal treatment, the peaks at -50 and 0°C disappear, the one at -145°C decreases and a peak at -35°C appears.

By X-ray diffraction, the presence of the monoclinic phase is observed in the mixture. But after the thermal treatment, almost all the monoclinic zirconia disappears and we observe a strong prevalence of the tetragonal zirconia (Fig. 6).

After complementary experiments carried out by low temperature TL under X-rays irradiation, it is reasonable to attribute the appearance of the TL peak at -35°C to the presence of tetragonal zirconia stabilised by the alumina.

The case of hibonite (Ca)

The study of high temperature thermoluminescence (TL) on non-doped alumina reveals the presence of three peaks called C' (around 70°C), D' (200°C) and E' (360°C) (Fig. 7). These are well known TL peaks in alumina [START_REF] Iacconi | Phénomènes luminescents des matériaux isolants[END_REF]. When we analyse, by thermoluminescence, the effect of both doping alumina by Ca 2+ and an oxidizing atmosphere of samples calcination, for a dopant amount of 1500 ppm, we still observe the presence of D' and E' peaks, but the C' peak seems to be moved forward from 70 to 110°C.

But in the case of a reducing atmosphere of doped alumina calcination, for the same amount of dopant, the C' peak does not move.

In that case, TL's intensity of the C' peak is represented according to the amount of calcium in the Fig. 8. For weak amount of dopant, the intensity strongly increases until a 1500 ppm concentration of Ca 2+ ; beyond, TL's signal rapidly decreases. This peak, which seems rather unstable by its intensity as much as by its position in temperature, presents the phenomenon of quenching by concentration. The more the amount of a defect is important in the lattice, the more the distance between two defects is weak. In case these defects intervene as centre of recombination, they are susceptible to re-absorb emitted photons during TL's analysis by a point defect of the Beyond 2500 ppm amount of Ca 2+ , TL's intensity of the C' peak increases linearly with the concentration of dopant (Fig. 8). We supposed so that the second part of this curve reports really the evolution of the thermoluminescence of a peak of another phase than the alumina. This new peak takes place at 110°C.

An analysis by X-rays diffraction, for both atmospheres of calcination, showed the presence of a second phase, the calcium hexaluminate CaAl 12 O 19 , so called up hibonite or CA6 phase.

We can summarize the amount of calcium from which the CA6 phase is discovered by X-rays diffraction, according to the atmosphere:

-Calcination under oxygen: 2500 ppm -Calcination under Ar+H 2 (2%): 5000 ppm In order to characterize its thermoluminescence, the hibonite was synthesised from a mixture of α-alumina [5] and calcium carbonate CaCO 3 in the isopropanol.

The high temperature thermoluminescence of the hibonite presents at least six peaks, as one can see on the DSC curve of the Fig. 9: a double peak between 90 and 130°C, three less intense peaks at 230, 345 and 415°C and a very weak intense peak around 280°C.

We have then compared the TL of the hibonite and the TL of alumina doped by 5000 ppm of calcium with a calcination atmosphere of O 2 (Fig. 9). For the doped sample, one finds the peak around 100°C, as well as a less intense one at 130°C, which characterise the presence of the CA6 phase. In the case of the oxidizing atmosphere, there is no uncertainty as regards the allocation of peaks, because the C' peak of alumina is not visible.

The linear increase of C' peak intensity with the amount of calcium under reducing atmosphere (Fig. 7) can be attributed to the double peak of the hibonite. There is however an ambiguity in the case of the reducing atmosphere: why is the double peak of the hibonite discovered at the same temperature as the C' peak, around 70°C ?

The answer probably lies in the instability of the thermoluminescence of the CA6 phase. Indeed, several analyses of thermoluminescence were realized successively on the same sample of alumina doped by 5000 ppm of calcium with a calcination atmosphere of Ar+2% H 2 (Fig. 10). We observe then very clearly that the C' peak of alumina moves forward to the highest temperatures under the successive irradiations influence, until it corresponds to the double peak of the hibonite.

Finally, it is necessary to note that the thermoluminescence of the hibonite is more intense than that of the α-alumina. This observation allows to explain that we can detect the CA6 phase even when it is in very small quantity in alumina.

Conclusions

The presence of zirconia tetragonal phase stabilized by alumina can be detected in alumina-zirconia composites by the thermoluminescence technique. Moreover, the detection limit of this technique is appreciably better than with the X-rays diffraction technique. Indeed, for an amount of zirconia lower than 10 mass% in the composite, XRD does not allow to detect the presence of tetragonal zirconia, while the peak at -35°C is still appearing for only 1 mass% of zirconia by TL method.

In the same way, the thermoluminescence is a valid technique to detect the presence of the hibonite phase in an alumina lattice.

Fig. 1

 1 Fig. 1 Thermal cycle for the treatment of the α-Al 2 O 3 -ZrO 2 mixtures

Fig. 2

 2 Fig. 2 Schematic diagram of thermoluminescence mechanism

Fig. 3 Fig. 4 X

 34 Fig.3Thermoluminescence after UV irradiation of mixtures before and after the thermal treatment (for 1wt%, 5wt% and 10wt% of ZrO 2 )

Fig. 5

 5 Fig. 5 Thermoluminescence after UV irradiation of a commercial mixture (allowing the stabilisation of zirconia in tetragonal form after sintering), containing 10 vol% of ZrO 2 , before and after the thermal treatment

Fig. 6 XFig. 7

 67 Fig. 6 X-rays diffraction of the commercial mixture (10 vol% of ZrO 2 ) before and after the thermal treatment

Fig. 8

 8 Fig. 8 Evolution in function of the amount of dopant, of the TL intensity of the peak situated between 70°C and 110°C for á-Al 2 O 3 doped by Ca 2+ , prepared by a thermal treatment under Ar + 2% H 2

Fig. 9 2 Fig. 10

 9210 Fig. 9 TL curves of the hibonite phase and of á-Al 2 O 3 doped by 5000 ppm of Ca 2+ , prepared by a thermal treatment under O 2