A Klimchik

Y Wu

C Dumas

S Caro

B Furet

A Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Identification of geometrical and elastostatic parameters of heavy industrial robots

Keywords: industrial robot, stiffness modeling, elastostatic calibration, gravity compensator, design of experiments

The paper focuses on the stiffness modeling of heavy industrial robots with gravity compensators. The main attention is paid to the identification of geometrical and elastostatic parameters and calibration accuracy. To reduce impact of the measurement errors, the set of manipulator configurations for calibration experiments is optimized with respect to the proposed performance measure related to the end-effector position accuracy. Experimental results are presented that illustrate the advantages of the developed technique.

I. INTRODUCTION

Further developments in aeronautic industry require highprecision and high-speed machining of huge aircraft components where industrial robots successfully replace conventional CNC-machines. For these applications, robots are quite attractive due to their large workspace (that can be easily extended) and high-speed motion capacity, as well as capability to process the parts with complex shape and geometry. However, because of high forces required for processing of modern aeronautic materials, the influence of robot stiffness becomes one of the key factors that essentially affect the position accuracy. To enhance robot stiffness, manufacturers tend to increase the manipulator link cross-sections that obviously leads to the augmentation of the robot mass. So, the gravity forces applied to the manipulator components become nonnegligible and also contribute to the position errors. To overcome this difficulty, the robots are usually equipped with different types of gravity compensators, which, however, considerably complicate the stiffness modeling of those heavy manipulators.

The problem of stiffness modeling for the heavy manipulators with gravity compensators has been in the focus of rather limited number of works. In contrast, for conventional serial manipulators without gravity compensators, the problem has been studied by a number of authors that considered both industrial and medical robots with essential compliance in the links and joints. [START_REF] Meggiolaro | Geometric and elastic error calibration of a high accuracy patient positioning system[END_REF][START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF][START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF]. Relevant works are mainly based on the virtual joint method (VJM), which was firstly proposed by Salisbury [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF] and lumped elastostatic properties of robot components in virtual springs. To our knowledge, the stiffness modeling for the manipulators with gravity compensators has not been studied in detail yet. Currently, the main activity in this area focuses on the gravity compensator design, which differs in kinematics [START_REF] Takesue | Design and Prototype of Variable Gravity Compensation Mechanism (VGCM)[END_REF] and/or may also employ some software tools embedded in the robot controller [START_REF] De Luca | A PD-type Regulator with Exact Gravity Cancellation for Robots with Flexible Joints[END_REF]. On the other hand, since the considered robots include closed loops induced by the compensators, some technique developed for the parallel manipulators can be adopted [START_REF] Company | A Method for Modeling Analytical Stiffness of a Lower Mobility Parallel Manipulator[END_REF][START_REF] Merlet | Parallel mechanisms and robots[END_REF][START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF]. This paper proposes a VJM-based stiffness model for a serial manipulator with a compensator attached to the second joint. It is assumed that the gravity compensation torque is generated by a spring incorporated in an additional link, which creates the closed loop to be included in the stiffness model. The main attention is paid to the identification of the model parameters and calibration experiment planning. The developed approach is confirmed by the experimental results that deal with compliance error compensation for robotic cell employed in manufacturing of large dimensional aircraft components.

To address these problems the remainder of the paper is organized as follows. Section 2 presents the stiffness modeling background. In Section 3, the geometrical and elastostatic models of the gravity compensator are presented. Sections 4 and 5 are devoted to the geometrical and elastostatic calibration. Section 6 summarizes the main contributions.

II. STIFFNESS MODELING BACKGROUND

Stiffness of a serial robot highly depends on its configuration and is defined by the Cartesian stiffness matrix that should be computed in the neighborhood of the loaded equilibrium configuration. Therefore, let us first present technique for computing the static equilibrium configuration and after focus on the computing of the stiffness matrix.

Using the VJM-based approach adopted in this paper, the manipulator can be presented as the sequence of rigid links separated by the actuators and virtual flexible joints incorporating all elastostatic properties of compliant elements [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF]. Since for the considered robots weights of the link is not negligible, it should be also incorporated in the model. In the frame of the VJM-based technique, it is convenient to replace the link weights i by the equivalent pair of the forces applied to the both link ends i and i . Further, after aggregation of all weight components applied to the same virtual joint, the external loading caused by the link weights can be described by the set of forces concentrated in the virtual

P  P  P i G
joint centers. For computational convenience, they can be collected in the matrix   1 ... n  GG G that will be referred to as "the auxiliary loading". This allows us to distinguish impact of the gravity from the influence of the conventional loading applied to the robot end-effector (that can be caused by the cutting forces in milling, for instance). Using those notations and the above defined assumptions the elastostatic model of the heavy serial robot can be presented as shown in Fig. 1. , where denote the vectors of the actuator and virtual joint coordinates respectively. It can be proved that applying the principle of virtual work and linearizing the geometrical model, the static equilibrium equations can be written in a matrix form as [11]

 , gqθ  θ , q θ      (G)T (F)T θθ    JG JF K θ
where θθ θ and θθ denote the Jacobians computed with respect to the relevant nodes where the forces are applied to, i.e.

(G) (1) ( ) T TT ... n    JJ J (F) ( ) n  JJ   () j  1n θ θ ,..., di KK K θ j
; and the matrix aggregates stiffnesses of all virtual springs.

,   Jg q θθ  θ ag 1 1 i   θ 
To find the desired static equilibrium configuration, the above system should be solved subject to the geometrical constraint , where the end-effector location is assumed to be given and the force is treated as an unknown. It is obvious that it is a dual problem compared to the traditional one where F is given and the vector is unknown, but this approach essentially reduces the computational effort. Relative iterative algorithm for solving this system can be presented as

 ,  tg q θ i  F t              1 (F) 1 (F)T θθ (F) (F) 1 (G)T θθ θ 1( G ) T ( F ) T θθ 1 θ 1 , • ii i ii i               FJ K J tg q θ J θ JK J G KJ G J F
This algorithm allows us to find the vectors F , defining deflections in the virtual springs under the loading F , which will be required for computing the stiffness matrix.

  () θ F ,   q θθ Jg , (F) (G) (G)T θθ θθ θθ θ    HHH J G θ , (G) 2 θθ 1 2 T j j j n     Hg G θ , (F θθ 2 )2 T   Hg F θ .
It should be noted that the above presented expressions have been derived for the case of manipulators without gravity compensators, where the matrix θ is constant. In the following Section, to take into account the compensator influence, these expressions will be modified by using the configuration dependent joint stiffness matrix describing properties of the virtual springs.

III. MECHANICS OF GRAVITY COMPENSATOR

The mechanical structure of the gravity compensator under study is presented in Fig. 2. The compensator incorporates a passive spring attached to the first and second links, which creates a closed loop that generates the torque applied to the second joint of the manipulator. This design allows us to limit the stiffness model modification by incorporating in it the compensator torque c M and adjusting the virtual joint stiffness matrix that here depends on the second joint variable only. The compensator geometrical model includes three node points P 0 , P  2 sq This mechanical design allows to balance the manipulator weight for any given configuration by adjusting the compensator spring preloading. It can be taken into account by introducing the zero-value of the compensator length 0 s corresponding to the unloaded spring. Under this assumption, the compensator force applied to the node P 1 can be expressed as follows

 0 •( ) sc s FK s      
where is the compensator spring compliance. c k Further, the angle  between the compensator links P 0 P 1 and P 1 P 2 (see Fig. 2) can be found from the expression

 2 sin s •( n i q as )       
which allows us to compute the compensator torque c M applied to the second joint

     0 •(1 )• • • ( /s i n cc ss L Ma K    2 ) q 2 q
Upon differentiation of the latter expression with respect to 2 , the equivalent stiffness of the second joint (comprising both the manipulator and compensator stiffnesses) can be expressed as:

q  22 0 •• c KKK L a       
where the coefficient

 2 2 0 22 2 sin cos co •( s )( ) ( q s aL qq s s            2 ) q  
highly depends on the value of joint variable 2 and the initial preloading in the compensator spring described by q 0 s . To illustrate this property, Fig. 3

0 0[m] s  0 0.2 [m] s  0 0.2 [m] s  0 0.4 [m] s  0 0.6 [m] s  0 0.8 [m] s  2 ,[deg] q  Figure 3
. Variation of the gravity compensator impact on the equivalent stiffness of the second joint Hence, using expression [START_REF] Merlet | Parallel mechanisms and robots[END_REF], it is possible to extend the classical stiffness model of the serial manipulator presented in Section II by modifying the virtual spring parameters in accordance with the compensator properties. While in the paper this approach has been used for the particular compensa-tor type (spring-based, acting on the second joint), the similar idea can be evidently applied to other compensator types.

Summarizing this Section, it is worth mentioning that the geometrical and elastostatic models of a heavy manipulator with a gravity compensator should include some additional parameters (  , , and c a L K , 0 s for the presented case) that are usually not included in datasheets. For this reason, the following Sections focus on the identification of the extended set of manipulator parameters.

IV. IDENTIFICATION OF COMPENSATOR GEOMETRY

In contrast to the serial manipulator that can be treated as a principal mechanism of the considered robots (whose geometry is usually defined in datasheets and can be perfectly tuned by means of calibration [START_REF] Daney | Interval method for calibration of parallel robots: Vision-based experiments[END_REF][13]), geometrical data concerning gravity compensators are usually not included in the technical documentation provided by the robot manufacturers. For this reason, this Section focuses on the identification of the geometrical parameters for the described above compensator mechanism (see Fig. 2).

A. Methodology

The geometrical structure of the considered gravity compensator is presented in Fig. 4. Its principal geometrical parameters are denoted as ,

L x a , y , where x a •cos aa   , •sin y aa  
(see notation in Section III). As follows from the figure, the identification problem can be reduced to the determination of relative locations of points P 0 and P 1 with respect to P 2 .

It is assumed that the measurement data are provided by the laser tracker whose "world" coordinate system is located at the intersection of the first and second actuated manipulator joints. The axes Y, Z of this system are aligned with the axes of joints #1 and #2 respectively, while the axis X is directed to ensure right-handed orthogonal basis. To obtain required data, there are several markers attached to the compensator mechanism (see Fig. 4). The first one is located at point P 1 , which is easily accessible and perfectly visible (the center of the compensator axis P 1 is exactly ticked on the fixing element). In contrast, for the point P 0 , it is not possible to locate the marker precisely. For this reason, several markers P 0i are used that are shifted with respect to P 0 , but located on the rigid component of the compensator mechanism (these markers are rotating around P 0 while the joint coordinate 2 is actuated). It should be noted that for the adopted compensator geometrical model (which is in fact a planar one), the marker location relative to the plane XY is not significant, since the identification algorithm presented in the following sub-section will ignore Z-coordinate. q Using this setting, the identification problem is solved in two steps. The first step is devoted to the identification of the relative location of points P 1 and P 2 . Here, for different values of the manipulator joint coordinates ,

2 { i q 1, } im  ,
the laser tracker provides the set of the vectors   1 i p describing the points that are located in an arc of the circle. After matching these points with a circle, one can obtain the desired value of (circle radius) and the Cartesian coordinates The second step deals with the identification of the relative location of points P 0 and P 2 . Relevant information is extracted from two data sets  01 and  02 that are provided by the laser tracker while targeting at the markers P 01 and P 02 .

Here, the points are matched to two circle arcs with the same center (explicitly assuming that the compensator model is planar), which yields the Cartesian coordinates 0 of the point P 0 , also with respect to the laser tracker coordinate system. Finally, the desired values

 i p  i p p x a ,
are computed as a projection of the difference vector on the corresponding axis of the coordinate system. As follows from the presented methodology, a key numerical problem in the presented approach is the matching of the experimental points with a circle arc. It looks like a classical problem, however, there is a particularity here caused by availability of additional data  2i describing relative loca- tions of the points  . This feature allows us to reformu- late the identification problem and to achieve higher accuracy compared with the traditional approach.

 q  1 i p

B. Identification algorithm

The above presented methodology requires solution of two identification problems. The first one aims at approximating of a given set of points (with additional arc angle argument) with an arc circle, which provides the circle center and the circle radius. The second problem deals with an approximation of several sets of points by corresponding number of circle arcs with the same center. Let us consider them sequentially.

To match the given set of points  with additional set of angles   i p  i q with a circle arc, let us define the affine mapping

 ii   pR u t  
where ii i denotes the set of reference points located on the unit circle whose distribution on the arc is similar to i , [cos , sin , 0] T qq  u p  is the scaling factor that defines the desired circle radius, R is the orthogonal rotation matrix, t is the vector of the translation that defines the circle center. It worth mentioning that such a formulation has an advantage (in the sense of accuracy) comparing to a traditional circle approximation and it is a generalization of Procrustes problem known from the matrix analysis.

Using equation [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF], the identification can be reduced to the following optimization problem

       1 ,, min T ii i ii m F           Rt pR u t pR u t
which should be solved subject to the orthogonality constraint T  RR I. After differentiation with respect to t , the latter variable can be expressed as ;

ii i i m i i i i m mm      pp p uu u    
Further, differentiation with respect to  yields to

 11 m TT ii i i m i     pR u uu i     
So, finally, after relevant substitutions the objective function can be presented as

   11 1 12 min mm m ii i TT T ii ii i i F                 R pp uu pR u    
where the unknown matrix must satisfy the orthogonality constraint and can be solved using SVD-decomposition of the matrix

 1 m i i T i    up UΣ V    
where the matrices are orthogonal and is the diagonal matrix of the singular values. Further, using the same approach as for the Procrustes problem, it can be proved that the desired rotation matrix can be computed as ;

m 1 1 m j jj j j j i i l i ii ii l T ms m       pp p p p p p j j l T       
Further, after differentiation with respect to 0 p , one can get the underdetermined system of linear equations

   11 0 11 1 2 km km j jj ii i i ji ji T j s      pp p p      whose solution  0 c   pp n  
is expressed via the rotation axis vector and a point belonging to this axis (here

n c p
 is an arbitrary scalar factor).

To solve this ambiguity, an additional objective should be defined It should be mentioned that practical application of the latter expression is essentially simplified by the adopted assumption concerning orientation of the reference coordinate system (see previous sub-section), the direction of the identified rotation axis is close to Z-direction. Hence, the developed algorithms allows us to identify the compensator geometrical parameters , L x a , y that are directly related to the above mentioned rotation center points P 0 , P 2 and corresponding radii. Below they will be applied to the processing of the experimental data.

a C. Experimental results

To demonstrate efficiency of the developed technique, the experimental study has been carried out. The experimental setup employed the robot KR-270 and the Leica laser tracker, which allowed us to measure the Cartesian coordinates of the markers attached to the compensator elements (see Figs 2,4). Six different manipulator configurations where considered that differed in the value of the joint angle 2 and three markers has been used. The experimental data are presented in Table I. q These data has been processed using the identification algorithm presenting in the previous sub-section. The obtained values for the parameters of interest , L x a , y are given in Table II. It also includes the confidence intervals computed as a 3  , where the standard deviation  has been evaluated using the Gibbs sampling. In the next section, the obtained model will be used for some transformations required during elastostatic calibration

V. ELASTOSTATIC PARAMETERS IDENTIFICATION

In contrast to geometrical calibration, where the manipulator and compensator can be considered independently, in elastostatic calibration the corresponding equations cannot be separated and the model parameters should be identified simultaneously. This section gives general ideas of the developed methodology and relevant identification algorithms, and also presents experimental results validating the proposed technique.

A. Methodology

In the frame of the adopted VJM-based modeling approach the desired stiffness parameters describe elasticity of the virtual springs located in the actuated joints of the manipulator, and also elasticity and preloading of the compensator spring (see Section III for details). Let us denote them as ,1 , 6 kj  j  for the manipulator joint compliances and 0 , c k s for the compliance and preloading of the compensator.

To find the desired parameters, the manipulator sequentially passes through several measurement configurations where the external loading is applied to the special endeffector described in Fig. 5 (it allows to generate both forces and torques applied to the manipulator). Using the laser tracker, the Cartesian coordinates of the reference points are measured twice, before and after loading. To increase identification accuracy, it is preferable to use several reference points (markers) and to apply the loading of the maximum allowed magnitude. It is worth mentioning that in order to avoid numerical singularities, the direction of the external loading should not be the same for all experiments (in spite of the fact that the gravity-based loading is the most attractive from the practical point of view). Thus, the calibration experiments yield the dataset that includes values of the ma-

B. Identification algorithm

To take into account the compensator influence while retaining our previous approach developed for serial robots without compensators, it is proposed to include in the second joint an equivalent virtual spring with non-linear stiffness depending on the joint variable 2 (see eq. ( 8)). Using this idea, it is convenient to consider several independent parameters 2i q k  corresponding to each value of 2 q . This allows us to obtain linear form of the identification equations that can be easily solved using standard least-square technique.

Let us denote the set of desired parameters 12 1 2 2 3 6 as the vector k and linearize (1) with respect to this vector. This allows us to present the relevant force displacement relations in the form , ( , ...), ,...,

kkk k k    () p ii  pB k
where matrices are composed of the elements of the matrix

() p i B  11 ,... 1, ) ,( TT i i i i ni ni i im    F AF JJ JJ   
that is usually used in stiffness analysis of serial manipulators [START_REF] Klimchik | Optimal Selection of Measurement Configurations for Stiffness Model Calibration of Anthropomorphic Manipulators[END_REF]. Here, ni denotes the manipulator Jacobian column, i F is the applied external force, and superscript '(p)' stands for the Cartesian coordinates (position without orientation). It is clear that transformation from i to J A () p i B is rather trivial and is based on the extraction from i the first three lines and inserting in it several zero columns.

A

Using these notations, the elastostatic parameters identification can be reduced to the following optimization problem

 0 ,, () () 1 () ()
   1 () () () 11 • T mm pp p ii i i ii           kB B B p T   
where the parameters 13 6 describe the compliance of the virtual joints #1,#3,...#6, while the rest of them 21 22 present an auxiliary dataset allowing to separate the compliance of the joint #2 and the compensator parameters c , ,..., kk k , ... kk 0 , k  . Using eq. ( 8), the desired expressions can be written as

     2 2 1 0 0 11 • qq i T cc i mm TT ii ii KK K K s         CC C  
where q is the number of different angles in the experimental data,

m 2 q   22 1• / • / • cos sin cos ii i LL a s aL as i           C   here 2 ii q     .
Thus, the proposed modification of the previously developed calibration technique allows us to find the manipulator and compensator parameters simultaneously. An open question, however, is how to find the set of measurement configurations that ensure the lowest impact of the measurement noise.

C. Design of calibration experiments

The main idea of the calibration experiment design is to select a set of robot configurations  i (and corresponding external loadings  q   i ) that ensure the best identification accuracy. The key issue here is the ranging of different plans in accordance with the prescribed performance measure. This problem has been already studied in the classical regression analysis [START_REF] Atkinson | Optimum Experi-ment Designs[END_REF], however the results are not suitable for the elastostatic calibration and require additional efforts.

F

In this work, it is proposed to use the industry oriented performance measure that evaluates the calibration plan quality. Its physical meaning is the robot positioning accuracy (under the loading), which is achieved after compliance error compensation based on the identified elastostatic parameters. It should be noted that usual approach (used in the classical design of experiments) evaluates the quality of experiments via covariance matrix of the identified parameters, which is does not have sense for our application.

Assuming that each experiment includes the additive measurement error i , the covariance matrix for the desired parameters can be expressed as

ε k       1 () () 1 1 () () () () 11 cov( ) E T TT m pp ii i mm pT p p p ii i i i i ii          kB B B εεBB B   m i n         Bk p Bk p

 

Following also usual assumption concerning the measurement errors (independent identically distributed, with zero expectation and standard deviation 2

 for each coordinate), the above equation can be simplified to

    1 2( ) ( ) 1 cov( ) T m pp ii i      kB B 
Hence, the impact of the measurement errors on the accuracy of the identified parameters k is defined by the matrix (in regression analysis it is known as the information matrix).

It is evident that in practice the most essential is not the accuracy of the parameters identification, but the accuracy of the robot positioning achieved using these parameters. Taking into account that this accuracy highly depends on the robot configuration (and varies throughout the workspace), it is proposed to evaluate the calibration accuracy in a certain given "test-pose" provided by the user. For the considered application, the test pose is related to the typical machining configuration 0 and corresponding external loading 0 related to the technological process. Let us denote the mean square value of the mentioned positioning error as 0 q F 2  and the matrix () p i A (see eq. ( 31)) corresponding to this test pose as .

() 0

p A It should be noted that that the proposed approach operates with a specific structure of the parameters included in the vector , where the second joint is presented by several components 21 22 while the other joints are described by a single parameter 13 6 . This motivates further rearrangement of the vector k and replacing it by several vectors  can be presented as

   1 () () 2 0 2) 0 1 1 () trace T q T m jp jp m p ii i j      0 ( p        AA A A  
Based on this performance measure, the calibration experiment design can be reduced to the following optimization problem

     1 () () 00 1 {,} () () 1 trace min T q T ii m jp jp i pp j i m i        qF AA A A subject to  max ,1 . i . F im  F  
whose solution gives a set of the desired manipulator configurations and corresponding external loadings. It is evident that its analytical solution can hardly be obtained and a numerical approach is the only reasonable one. More detailed description of the developed technique providing generation of the optimal calibration plans is presented in our previous publication [START_REF] Klimchik | Optimal Selection of Measurement Configurations for Stiffness Model Calibration of Anthropomorphic Manipulators[END_REF], where the problem of gravity compensator modeling had not been addressed yet.

D. Experimental results

Similar to the previous section, the developed technique has been applied to the robot KR-270. The parameters of interest were the compliances j k of the actuated joints and the elastostatic parameters 0 , c k s of the gravity compensator. To generate elastostatic deflections, the gravity forces 250-280 kg have been applied to the robot end-effector (see Fig 6). The Cartesian coordinates of three markers located on the tool (see Fig, 5) have been measured before and after the loading. To find optimal measurement configurations, the design of experiments has been carried out for five different angles 2 that are distributed between the joint limits. For each 2 three optimal measurement configurations have been found taking into account physical constraints that are related to the joint limits and the possibility to apply the gravity force (work-cell obstacles and safety reasons). The results of the calibration experiment design are presented in Table III. Here, values of 1 have been chosen to ensure good visibility of the markers for the laser tracker. To ensure identification accuracy for each configuration, the experiments were repeated three times. In total, 405 equations were considered for the identification, from which 7 physical parameters have been obtained (because of page limit, the experimental data cannot be presented in the paper ).

Gravity compensator

q q q TABLE III. OPTIMAL MEASUREMENT CONFIGURATIONS

Joint angles, [deg]

q 1 q 2 q 3 q 4 q 5 q 6 The obtained experimental data have been processed using the identification algorithm presented in sub-section V.B. Corresponding values of the gravity compensator and the manipulator elastostatic parameters are presented in Table IV. It also includes the confidence intervals computed as 3

 ,

where the standard deviation  has been evaluated using Gibbs sampling. Using the obtained results it is possible to identify an equivalent non-linear spring 22 (see Fig, 7) that is used in the stiffness modeling of the manipulator with the gravity compensator (see Section II).

() kq

The identified joint compliances can be used to predict robot deformations under the external loading. The identification errors of the joint compliances vary from 1.3% to 4.9%, where the lowest errors have been achieved for the joints #2 and #3. The highest errors are in the joints #4-#6. So, higher precision has been achieved for the joints that are closer located to the manipulator base. This fact is due to the different joint errors contribution to the total positioning error, which have been minimized while design of calibration experiments. Comparatively low identification accuracy for the compensator spring is caused by limited number of different angles 2 q . Fig. 7 shows that due to the gravity compensator the equivalent compliance of the second joint has been decreased comparing to the compliance of the second joint of the corresponding serial manipulator. The paper presents a new approach for the identification of the elastostatic parameters of heavy industrial robots with the gravity compensator. It proposes a methodology and data processing algorithms for the identification of both geometrical and elastostatic parameters of gravity compensator and manipulator. To increase the identification accuracy, the design of experiments has been used aimed at proper selection of the measurement configurations. In contrast to other works, it is based on a new industry oriented performance measure that is related to the robot accuracy under the loading.

The advantages of the developed techniques are illustrated by experimental study of the industrial robot Kuka KR-270, for which the joint compliances and parameters of the gravity compensator have been identified.

FFigure 1 .

 1 Figure 1. Virtual joint model of a serial industrial robot For this serial chain, the end-effector location and locations of the node centers t

θ

  To find the corresponding stiffness matrix C K , the forcedeflection relation (1) should be linearized in the neighborhood of the equilibrium configuration . After relevant transformations [11], one can get the following expression (, F θ,) t which includes both the first and second order derivatives (Jacobians and Hessians) of the functions   , j gq θ describing the manipulator geometry:

Figure 2 .

 2 Figure 2. Gravity compensator and its model

  ,  and the distances x a and y , whose geometrical meaning is clear from Fig.2. Using these notations, the variable a s describing the compensator spring deflection can be computed from the expression 

presents a set of curves 2 (

 2 ) q  obtained for different values of 0 s (the remaining parameters  , , correspond to robot KUKA KR-270 studied in the experimental part of the paper, see Section IV).

4 RFigure 4 .

 44 Figure 4. Geometrical parameters of the gravity compensator and location of the measurement points labeled with markers

  Since the matrix R is included in the second term only, the problem can be further simplified to

  to find the scaling factor  defining the arc radius and the vector defining the arc center.t The second problem aims at approximating of several point sets     1 ,..., k pp ii by corresponding number of concen- tric circle arcs with the same center 0 p . It should be noted that here the data set   i is not useful, since the required angles q   i  are not measured directly and cannot be computed without having exact compensator geometry. In this case, the objective function can be written in a straightforward way be proved that this optimization problem does not lead to a unique solution (in fact, it gives the rotation axis passing through the desired center). For this problem, differentiation of the objective function F with respect to2

FFigure 5 .

 5 Figure 5. End-effector used for elastostatic calibration experiments

F

  which leads to the following solution

Figure 6 .

 6 Figure 6. Experimental setup for the identification of the elastostatic parameters

  A. Klimchik, D. Bondarenko, A. Pashkevich, S. Briot, B. Furet, of tool deflection in robotic-based Milling," the 9th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2012), July 28-31, 2012, Rome, Italy, pp. 113-122.

Figure 7 .

 7 Compliance of equivalent non-linear spring in the second joint VI. CONCLUSIONS

TABLE I

 I

		.	EXPERIMENTAL DATA FOR GEOMETRICAL CALIBRATION
	q 2		P 1		P 01		P 02
	[deg]	x, [mm]	y, [mm]	x, [mm]	y, [mm]	x, [mm]	y, [mm]
	-0.01		-31.84	183.86 -872.10 -125.38 -813.50 -255.59
	-30 -118.44	143.42 -872.30 -126.07 -813.33 -256.18
	-60 -173.30	65.12 -872.50 -109.90 -825.09 -244.64
	-90 -181.76	-30.14 -868.43	-78.20 -844.66 -219.04
	-120 -141.45 -116.82 -858.90	-47.60 -859.43 -190.44
	-145		-78.10 -165.47 -852.53	-33.68 -864.66 -176.01
	TABLE II.	GEOMETRICAL PARAMETERS OF GRAVITY COMPENSATOR
				L, [mm]	a x , [mm]	a y , [mm]
	value		184.72	685.93	120.30
	CI		±0.06	±0.70		±0.69

of the point P 2 (circle center) with respect to the laser tracker coordinate system.

ACKNOWLEDGMENT

The work presented in this paper was partially funded by the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO). The authors also thank Fabien Truchet, Guillaume Gallot, Joachim Marais and Sébastien Garnier for their great help with the experiments.

*The work presented in this paper was partially funded by the Agence Nationale de la Recherche (ANR), France (Project ANR-2010-SEGI-003-02-COROUSSO).


Alexandr Klimchik
email: alexandr.klimchik@mines-nantes.fr

Anatol Pashkevich

Damien Chablat

CAD-based approach for identification of elasto-static parameters of robotic manipulators

Keywords: Robotic manipulator, Stiffness modeling, Elasto-static parameters, CAD-based Identification, Finite element analysis, Stiffness matrix

The paper presents an approach for the identification of elasto-static parameters of a robotic manipulator using the virtual experiments in a CAD environment. It is based on the numerical processing of the data extracted from the finite element analysis results, which are obtained for isolated manipulator links. This approach allows to obtain the desired stiffness matrices taking into account the complex shape of the links, couplings between rotational/translational deflections and particularities of the joints connecting adjacent links. These matrices are integral parts of the manipulator lumped stiffness model that are widely used in robotics due to its high computational efficiency. To improve the identification accuracy, recommendations for optimal settings of the virtual experiments are given, as well as relevant statistical processing techniques are proposed. Efficiency of the developed approach is confirmed by a simulation study that shows that the accuracy in evaluating the stiffness matrix elements is about 0.1%.

Introduction

Current trends in mechanical design of robotic manipulators are targeted at essential reduction of moving masses, in order to achieve high dynamic performances with relatively small actuators and low energy consumption. This motivates using advanced kinematical architectures (Orthoglide, Delta, Gantry-Tau, etc.) and light-weight materials, as well as minimization of the cross-sections of all manipulator elements [START_REF] Aviles | Optimization of linkages for generalized rigid-body guidance synthesis based on finite element models[END_REF][START_REF] Siciliano | Springer Handbook of Robotics[END_REF][START_REF] Merlet | Parallel Robots[END_REF][START_REF] Tyapin | Kinematic and elastostatic design optimization of the 3-DOF Gantry-Tau parallel kinamatic manipulator[END_REF]. The primary constraint for such minimization is the mechanical stiffness of the manipulator, which is directly related with the robot accuracy defined by the design specifications. For this reason identification (or evaluation) of the manipulator elasto-static parameters becomes one of the key issues in development and optimization of modern robotic systems Similar to general structural mechanics, the robot stiffness characterizes the manipulator resistance to the deformations caused by an external force or torque applied to the end-effector [START_REF] Duffy | Statics and Kinematics with Applications to Robotics[END_REF][START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF]. Numerically, this property is usually defined through the stiffness matrix K , which is incorporated in a linear relation between the translational/rotational displacement and static forces/torques causing this transition (assuming that all of them are small enough). The inverse of K is usually called the compliance matrix and is denoted as k . As it follows from related works, for conservative systems K is a 6 6  semi-definite non-negative symmetrical matrix but its structure may be non-diagonal to represent the coupling between the translation and rotation [START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF].

In stiffness modeling of robotic manipulator, because of some specificity, there are some particularities in terminology. In particular, the matrix K is usually referred to as the "Cartesian Stiffness Matrix" C K and it is distinguished from the "Joint-Space Stiffness Matrix" θ K that describes the relationship between the static forces/torques and corresponding deflections in the joints [START_REF] Ciblak | Synthesis of Cartesian Stiffness for Robotic Applications[END_REF]. Both of these stiffness matrices can be mapped to each other using the Conservative Congruency Transformation [START_REF] Kao | Properties of grasp stiffness matrix and conservative control strategy[END_REF], which is trivial if the external (or internal) loading is negligible. In this case, the transformation is entirely defined by the corresponding Jacobian matrix. However, if the loading is essential, it is described by a more complicated equation that includes both the Jacobian as well as the Jacobian derivatives and the loading vector [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF][START_REF] Sh | Geometrical Approach to The Conservative Congruence Transformation (CCT) for Robotic Stiffness Control[END_REF]. Other specific cases, where the above transformation is non-trivial (non-linear or even singular), are related to manipulators with passive joints and over-constrained parallel architectures [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF].

In the most general sense, existing approaches to the manipulator stiffness modeling may be roughly divided into three main groups: (i) the Finite Elements Analyses (FEA), (ii) the Matrix Structural Analyses (MSA), and (iii) the Virtual Joint Modeling method (VJM). Their advantages and disadvantages are briefly presented below.

An evident advantage of the FEA-modeling is its high accuracy that is limited by the discretization step only. For robotic application it is very attractive, since the links/joints are modeled with its true dimension and shape [START_REF] Corradini | Evaluation of a 4 degree of freedom parallel manipulator stiffness[END_REF][START_REF] Avilés | A procedure based on finite elements for the solution of nonlinear problems in the kinematic analysis of mechanisms[END_REF][START_REF] Fernández-Bustos | A new finite element to represent prismatic joint constraints in mechanisms[END_REF]. However, while increasing of the number of finite elements, the problem of limited computer memory and the difficulty of the high-dimension matrix inversion become more and more critical. Besides the high computational efforts, this matrix inversion generates numerous accumulative round-off errors, which reduce accuracy. In robotics, this causes rather high computational expenses for the repeated re-meshing and re-computing, so in this area the FEA method is usually applied at the final design stage only [START_REF] Long | Optimal structural design of a planar parallel platform for machining[END_REF][START_REF] Rizk | A comparative stiffness analysis of a reconfigurable parallel machine with three or four degrees of mobility[END_REF]. Nevertheless, in this work this method is applied for the links stiffness matrix identification that can be used further in the frame of the VJM and MSA techniques. Such combination allows us to use the advantages of the FEA while avoiding intensive computations for different manipulator configurations.

Matrix Structural Analysis method incorporates the main ideas of the FEA but operates with rather large compliant elements such as beams, arcs, cables, etc. [START_REF] Martin | Introduction to matrix methods of structural analysis[END_REF]. This obviously leads to the reduction of the computational expenses and, in some cases, allows us to obtain an analytical stiffness matrix for the specific task. For parallel robots, this method has been developed in works [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF][START_REF] Ecorchard | Elasto-geometrical modeling and calibration of redundantly actuated PKMs Original Research Article[END_REF], where a general technique for stiffness modeling of the manipulator with rigid/flexible links and passive joints was proposed. It has been illustrated by stiffness analysis of parallel manipulator of Delta architecture where the links were approximated by regular beams. The latter causes some doubts in the model accuracy compared to the combination of the FEA and VJM techniques that are being developed here

The core of Virtual Joint Modeling method is an extension of the conventional rigid-body model of the robotic manipulator, where the links are treated as rigid bodies but the joints are assumed to be compliant (in order to accumulate all types of existing flexibilities in the joints only). Geometrically, such approximation is equivalent to adding to the joints some auxiliary virtual joints (with embedded virtual springs). It is obvious that such lumped presentation of the manipulator stiffness (that in reality is distributed) essentially simplifies the model. So, at present it is the most popular stiffness analysis method in robotics. This method was first introduced by Salisbury [START_REF] Salisbury | Active Stiffness Control of a Manipulator in Cartesian Coordinates[END_REF] and Gosselin [START_REF] Gosselin | mapping for parallel manipulators[END_REF], who assumed that the main flexibility sources were concentrated in the actuator joints. The derived expression defining relation between the joint and Cartesian stiffness matrices (Conservative Congruency Transformation) became a basis for the manipulator stiffness analysis in many research works. Later, these results have been further developed in order to take into account some specific geometrical constrains [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF] and external loading [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]. Nevertheless, external loading is assumed small enough to detect any non-linear effects discovered in this work. Due to its computational simplicity, the VJM method has also been successfully applied to the analytical stiffness analysis of a translational parallel manipulator [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF]. A key issue of this method is how to define the virtual spring parameters. At the beginning, it was assumed that each actuated joint is presented by a single one-dimensional virtual spring [START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF]. Further, to take into account the links flexibilities, the number of virtual joints was increased and in each actual actuated or passive joint several translational and rotational virtual springs were included [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF]. The latest developments in this area operate with 6-dimensinal virtual springs identified using the FEA-based method [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. This leads to essential increasing of the VJM method accuracy that becomes comparable with the accuracy of the FEA-based techniques, but with much lower computational expenses. Hence, both for the MSA and VJM techniques the problem of accurate stiffness matrices of the links with a complex shape become a critical.

In the first works, it was explicitly assumed that the main sources of elasticity are concentrated in actuated joints. Correspondingly, the links were assumed to be rigid and the VJM model included one-dimensional springs only. In other recent works, compliance of the links has been taken into account by introducing additional virtual joints describing their longitudinal elasticity [START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF] or stiffness properties in several directions [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF]. Recent development in this area use 6-dimensional virtual joints to describe elasticity of each link [START_REF] Li | Stiffness analysis for a 3-PUU parallel kinematic machine[END_REF].

In the most of precious works, the stiffness parameters of the virtual joints describing the link elasticity (and incorporated in the matrix θ K ) were evaluated rather roughly, using a simplified representation of the link shape by regular beams. Besides, it was assumed that all linear and angular deflections (compression/tension, bending, torsion) are decoupled and are presented by independent one-dimensional springs that produce a diagonal stiffness matrix of size 66 for each link. Latter, this elasticity model was enhanced by using complete 66 non-diagonal stiffness matrix of the cantilever beam [START_REF] Li | Stiffness analysis for a 3-PUU parallel kinematic machine[END_REF]. This allowed taking into account all types of the translational/rotational compliance and relevant coupling between different deflections. Other enhancements include the link approximation by several beams, but it gives rather modest improvement in accuracy.

Figure 1

Integration of VJM modeling approach with FEA-based identification technique of the stiffness model parameters Further advance in this direction (applicable to the links of complicated shape) led to the FEA-based identification technique that involves virtual loading experiments in CAD environment and stiffness matrix estimation using dedicated numerical routines [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. The latter essentially increased accuracy of the VJM-modeling while preserving its high computational efficiency. It is worth mentioning that usual high computational expenses of the FEA is not a critical issue here, because it is applied only once for each link (in contrast to the straightforward the FEA-modeling for the entire manipulator, which requires complete re-computing for each manipulator posture). As a result, this approach allowed the authors to integrate accuracy of the FEA-modeling into the VJM-modeling technique that provides high computational efficiency. General methodology of this hybrid approach is presented in Figure 1.

However, in spite of good results obtained for some case studies, there are still a number of open questions in the FEA-based identification technique. They include optimal setting of the virtual experiment (i.e. definition of the mesh parameters, the joint contact surfaces to apply forces, etc.) and enhancement of numerical algorithms used for computation of the stiffness matrix elements (increasing robustness with respect to FEA-modeling errors, distinguishing zero elements from small ones, etc.). These issues have never been given proper attention in previous works (only some preliminary results have been presented in our previous work [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]) and will be in the focus of this work.

Problem statement

In robotic literature, there are several ways to obtain the stiffness matrix of a link. The simplest of them assumes that a real link (with rather complex shape) can be approximated by a simple beam, for which the stiffness matrix can be easily expressed analytically. Another, more accurate, technique deals with multi-beam approximation where the link is presented as a serial chain of rigid bodies separated by virtual springs. For this presentation, the stiffness matrix of the whole link is computed using a common procedure known from stiffness analysis of serial manipulators. However, in spite of computational convenience and better accuracy, the second approach can be hardly applied to many industrial manipulators (where links may be nonhomogeneous, their shape is quite complex and cross-section is non-constant, etc). For instance, for the Orthoglide foot (Figure 2), even a four-beam approximation provides an accuracy of only 30-50% [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF].

Figure 2 CAD model of the link with complex shape

To achieve desired accuracy, here it is proposed to apply the FEA-based identification methodology that was firstly introduced in [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] but needs some further development. The corresponding algorithm is schematically presented in Figure 3. At the first step, a CAD model of the link is created, which properly describes the link shape, cross-sections, distribution and physical properties of the material (density, Young's modulus, Poisson's ratio, etc.). Then, at the second step, one of the link ends is fixed in accordance with contact surfaces of an adjacent element (for example by cylindrical surface for the revolute joint). At the second end of the link, certain distributed (or localized) force/torque is applied in accordance with contact surfaces of an adjacent element. For these settings, the FEA-modeling is performed which yields the deflection field for a huge number of points generated by the meshing procedure. From this set, a subset (so called 'deflection field') is extracted corresponding to the neighborhood of the reference point. This field contains desired information on the link compliance with respect to the applied force/torque. Such virtual experiments are repeated several times, for different directions of forces and torques (issue of their magnitude is discussed in the following sections). And finally, at the third step, the proposed identification procedure that gives the desired stiffness matrix of the link is applied.

It is worth mentioning that the above described algorithm has some essential advantages, which are not achievable while using any analytical technique. In particular, here it is possible to take into account (straightforwardly and explicitly) the joint elasticity that usually has a number of particularities, such as significant compliance of link/joint areas located closed to the contact points or surfaces [START_REF] Redon | An efficient, error-bounded approximation algorithm for simulating quasi-statics of complex linkages[END_REF]. For instant, for a case study presented in this paper, some of the stiffness matrix elements are reduced by the factor of 10-12 if the joint particularities are modeled properly. It is evident that existing analytical expressions do not take into account these issues.

CAD model

FEA-based virtual experiments

Stiffness matrix identification

Real shape Material properties

Joints particularities

Couplings between translational\rotational deflections

Deflections field

Stiffness matrix

Figure 3 Algorithm for stiffness matrix identification procedure Hence, the goal of this paper is the accuracy improvement of the stiffness matrix identification for manipulator elements. As it follows from the related analysis, particular problems that should be considered here are the following:

(i)
Further development of the FEA-based methodology for the link stiffness matrix identification, which allows us to take into account the link shape, coupling between rotational/translational deflections and joint particularities;

(ii) Development of numerical technique, which computes the stiffness matrix from the set of points (displacement field) extracted from the FEA-based virtual experiments;

(iii) Evaluation of this technique accuracy with respect to virtual experiment settings (meshing type\size, deflection range, fixing and force application method) and size\shape of virtual sensor;

(iv) Minimization of the identification errors by statistical processing of the experimental data (outlines elimination, determination of the confidence intervals and detecting "zero" elements of the stiffness matrix);

(v) Validation of the developed technique on application examples related to typical parallel manipulators and their comparison with conventional regular-shape approximation models.

To address these problems, the remainder of this chapter is organized as follows. Section 3 introduces the FEA-based methodology for identification of the link stiffness matrix. Section 4 proposes a numerical technique for evaluating the stiffness matrix elements from the field of points. Section 5 focuses on the accuracy estimation. Section 6 deals with minimization of the identification errors. In Section 7, the developed technique is applied to the links of Orthoglide manipulator. And finally, Section 8 summarizes main results and contributions of this paper.

FEA-based approach for identification of link stiffness matrix

Let us start from a detailed description of the FEA-based methodology for the link stiffness matrix identification. It is based on a number of the virtual experiments that are conducted in the CAD-environment (CATIA with ANSYS, for instance). Each of these experiments gives some information on the resistance of an elastic body or mechanism to deformations caused by an external force or torque. The desired stiffness model is obtained by integrating data obtained from several different experiments, which differ in the direction of applied forces/torques. For relatively small deformations, the stiffness properties are defined through the so-called stiffness matrix K, which defines the linear relation δ

             Fp K M φ (1)
between the three-dimensional translational/rotational displacements T ( , , )

x y z p p p  p ; T δ δ δ ( , , ) δ x y z     φ
and the static forces/torques T ( , , )

x y z F F F  F , T ( , , ) x y z M M M  M
causing this transition. As known from mechanics, K is a 66 symmetrical semi-definite non-negative matrix, which may include non-diagonal elements to represent the coupling between the translations and rotations [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. The inverse of K is usually called the compliance matrix and is denoted as k . For robotic manipulators, the matrix K can be computed semi-analytically provided that the stiffness matrices of all separate components (links, actuators, etc.) are known with desired precision [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. However, explicit expressions for the link stiffness matrices can be obtained in simple cases only (truss, beam, etc.). For more sophisticated shapes that are commonly used in robotics, the stiffness matrix is usually estimated via the shape approximation, using relatively small set of primitives [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF]. However, as it follows from the studies [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF], the accuracy of this approach is rather low. Hence, in a general case, it is prudent to apply to each link the FEA-based techniques, which hypothetically produce rather accurate result.

Figure 4

Identification experiment for the element with a complex shape Using the FEA, the stiffness matrix K (or its inverse k ) is evaluated from several numerical experiments, each of which produces the vectors of linear and angular deflections (p, ) corresponding to the applied force and torque (F, M) (see Figure 4). Then, the desired matrix is computed from the linear system ..

mm mm               F F p p k MM φφ ( 2
)
where m is the number of experiments ( 6 m  ) and the matrix inverse is replaced by the pseudo-inverse in the case of 6 m  . It is obvious that the case of 6 m  with special arrangement of the forces and torques is numerically attractive

                        1 1 4 4 2 2 5 5 3 3 6 6
, 0, 0 0, 0, 0 ; 0, 0, 0 , 0, 0 0, , 0 0, 0, 0 ; 0, 0, 0 0, , 0 0, 0, 0, 0, 0 ; 0, 0, 0 0, 0,

T T T T xz TT TT yy T T T T yz FM FM FM             F M F M F M F M F M F M . ( 3
)
corresponding to the diagonal structure of the matrix to be inverted. In this case, each FEA-experiment produces exactly one column of the compliance matrix

1 2 3 4 5 6 1 2 3 4 5 6 δ δ δ δ δ δ / / / / / / / / / / / / x y z x y z x y z z y z F F F M M M F F F M M M     p p p p p p k φ φ φ φ φ φ (4)
and the values (p i ,  i ) may be clearly physically interpreted. On the other hand, by increasing the number of experiments ( 6 m  ) it is possible to reduce the estimation error. Besides, it is worth mentioning that the classical methodology of the optimal design of an experiment cannot be applied here directly, because it is not possible to include in equation ( 4) the measurement errors induced by virtual experiments as additive components. Some aspects of this problem are studied in section 5.4.

Hence, to obtain the desired stiffness (compliance) matrix, it is required to estimate first the deflections (p, ) corresponding to each virtual experiment. This issue is in the focus of the following section.

Numerical technique for evaluating the stiffness matrix elements

Usually, in FEA-based experiments, the values (p, ) are computed from the spatial location of a single finite element enclosing the reference point. In contrast to this approach, here it is proposed to evaluate (p, ) from the set of points (displacement field) describing transitions of a rather large number of nodes located in the neighborhood of reference point. It is reasonable to assume that such modification will yield positive impact on the accuracy, since the FEA-modeling errors usually differ from node to node, exposing almost quasi-stochastic nature.

Related optimization problem

To formulate this problem strictly, let us denote the displacement field by a set of vector couples { , | 1, } ii in  pp (see Figure 5) where the first component i p defines the node's initial location (before applying the force/torque), i p refers to the node's displacement due to the applied force/torque, and n is the number of considered nodes. Then, assuming that all the nodes are close enough to the reference point, this set can be approximated by a "rigid transformation"

( ) , 1, δ i i i in       p p R φ p p , ( 5
)
that includes as the parameters the linear displacement vector p and the orthogonal 33 matrix R that depends on the rotational displacement . Then, the problem of the deflection estimation can be presented as the best fit of the considered vector field by equation ( 5) with respect to the six scalar variables ( , , , δ , , ) δδ In practice, the FEA-modeling output provides the deflection vector fields for all nodes referring to all components of the mechanism. So, it is required to select a relevant subset corresponding to the neighborhood of the reference point p 0 . Besides, the node locations p i must be expressed relative to this point, i.e. the origin of the coordinate system must be shifted to p 0 . The latter is specified by the physical meaning of the deflections in the stiffness analysis.

To estimate the desired deflections (p, ), let us apply the least square technique that leads to minimization of the sum of squared residuals 2 1 ( min δ )

n i i i , i f         Rp p p R φ p p (6)
with respect to the vector t and the orthogonal matrix R representing the rotational deflections . The specificities of this problem (that does not allow direct application of the standard methods) are the orthogonal constraint T  R R I and non-trivial relation between elements of the matrix R and the vector . The following sub-sections present two methods for computing p, , as well as their comparison study.

SVD-based solution

For the comparison purposes, first let us briefly present the exact solution of the optimization problem [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF]. It relies on some results from the matrix algebra referred to as "Procrustes problem" [START_REF] Gower | Procrustes Problems[END_REF]. The corresponding estimation procedure is decomposed in two steps, which sequentially produce the rotation matrix R and the translation vector p. Then, the desired vector of rotation angles  is extracted from R.

Let us transform the original optimization problem [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF] to the standard form. First, for the non-constrained variable p, straightforward differentiation and equating to zero gives an expression 11 1 ()

nn ii ii n          p p R I p (7)
Then, after relevant substitution and denoting

1 1 ˆn i i i i n    p p p ; 1 1 ˆ() n i i i i i i n         g p p p p , (8)
the original optimization problem is reduced to the orthogonal Procrustes formulation

2 1 ˆˆmin n ii i f      R g R p . ( 9
)
with the constraint T  R R I . The latter yields the solution [31] T  R VU [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF] that is expressed via the singular value decomposition (SVD) of the matrix

1 • n TT ii i   p g U ΣV , (11)
which requires some computational efforts. Hence, the above expressions [START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF], [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF] allow to solve the optimization problem (6) in terms of variables R and p.

Further, to evaluate the vector of angular displacements , the orthogonal matrix R must be decomposed into a product of elementary rotations around the Cartesian axes x, y, z. It is obvious that, in a general case, this decomposition is not unique and depends on the rotation order. However, for a small  (that is implicitly assumed for FEA-experiments) this matrix may be uniquely presented in differential form δδ δδ Using this expression, the desired parameters δ δ δ ,,

x y z    may be extracted from [] ij r  R
in several ways (Table 1), which are formally equivalent but do not necessarily possess similar robustness with respect to round-off errors. A relevant comparison study will be presented in Section 5.2.

LIN-based solution

To reduce the computational efforts and to avoid the SVD, let us introduce linearization of the rotational matrix R at the early stage, using explicit parameterization given by expression [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF]. This allows rewriteing the equation of the "rigid transformation" (5) in the form

δ ; 1, ii in      pp φp (13)
that can be further transformed into a linear system of the following form

  ;1 δ , ii in       p I P p φ ( 14
)
where i P is a skew-symmetric matrix corresponding to the vector :

i p 0 0 0 zi yi i zi xi yi xi pp pp pp         P ( 15
)
Then, applying the standard least-square technique with the objective

2 1 mi δ n n ii , i f      
 φt pP φp [START_REF] Long | Optimal structural design of a planar parallel platform for machining[END_REF] one can get the solution

1 11 1 1 1 δ nn ii ii n n n T T T i i i i i i i i n                                       I P p p φ P P P P p (17)
that employs the 66 matrix inversion. This solution can be simplified by shifting the origin of the coordinate system to the point

1 1 n ci i n     pp leading to expression 1 1 1 1 1 ˆδ ˆn i i n n T T ii ii i i n                           I0 p p φ 0 P P Pp (18)
that requires inversion of the matrix of size 33. Here, ˆi P denotes a skew-symmetric matrix corresponding to the vector ˆi

i c  p p p .
It is evident that the obtained expression ( 18) is more computationally efficient compared to ( 8) and [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF]. Besides, for some typical cases corresponding to regular meshing patterns of the FEA, expression ( 18) can be further simplified as shown in the following subsection.

Analytical expressions for typical case studies

In practice, the points i p are distributed in the space in a regular way, in accordance with the meshing options chosen for the FEA-modeling. This allows us to inverse the matrix 1 ˆn T ii i  PP analytically and to obtain very simple expressions for p, . Let us consider several patterns that are useful for practice and are presented in Figure 6.

Figure 6 Typical patterns of the deflection field

Case 1: Symmetrical field (Figure 6a). If the field is symmetrical with respect to its centre c p , the solution (18) can be presented in a compact analytical form as

1 1 n i i n     pp ; 1 1 δ n T ii i     φ D P p ( 19
)
where the matrix

2 2 2 2 2 2 1 1 1 ˆˆˆˆˆ( ) ( ) ( ) n n n yi zi xi zi xi yi i i i diag p p p p p p              D ( 20
)
is diagonal and easily inverted.

Case 2: Cubic field (Figure 6b). If the field is symmetrical and, in addition, it is produced by uniform meshing of the cubic subspace aaa , the matrix D is expressed as

d  DI where 2 ( 1) 6( 1) d a n n n    .
Case 3: Planar square field (Figure 6c). For the field produced by uniform meshing of the square aa  located perpendicular to the x-axis, the expression for the matrix

  / 2 / 2 diag d d d  D
Summarizing Section 4, it should be noted that the proposed numerical technique is computationally attractive and allow simultaneous estimation both translational and rotational deflections p,  from the FEA-produced field. Below it is evaluated from the point of view of the precision and robustness.

Accuracy of the stiffness matrix identification

To be applied in practice, the accuracy of the developed numerical technique should be evaluated with respect to the virtual experiment settings (meshing type\size, deflection range, fixing and force application method) and size\shape of virtual sensor. Let us consider several case studies corresponding to typical industrial applications.

Error sources and their impact

It is obvious that the main source of estimation errors is related to the FEA-modeling that highly depends on the size and type of the finite elements, meshing options, incorporated numerical algorithms, computer word length and the round-off principle. Hypothetically, the accuracy be essentially improved by reducing the mesh size and increasing the number of digits in presentation of all variables. But there are some evident technical constraints that do not allow ignoring the FEA limitations.

Another type of errors arises from numerical differentiation incorporated in the considered technique. Strictly speaking, the linear relation ( 1) is valid for rather small deflections that may be undetectable against the FEA-modeling defects. On the other side, large deflections may be out of the elasticity range. Hence, it is prudent to find a compromise for the applied forces/torques taking into account both factors.

Influence of linearization and round-offs

Since both of the considered algorithms (SVD-based and the proposed one, see Section 4) involve numerous matrix multiplications, they may accumulate the round-off errors. Besides, they employ the first-order approximation of the matrix R that may create another source of inaccuracy. Hence, it is prudent to obtain numerical assessments corresponding to a typical case study. For these assessments, there were examined data sets corresponding to the cubic deflection field of size 101010 mm 3 with the mesh step of 1 mm (1331 points). The deflections have been generated via the "rigid transformation" [START_REF] Duffy | Statics and Kinematics with Applications to Robotics[END_REF] As it follows from the analysis, the influence of the linearization and round-offs is negligible for the translation (the induced errors are less than 10 -14 mm). In contrast, for the rotation, practically acceptable results may be achieved for rather small angular deflections that are less than 1.0° (the errors are up to 0.01°). The latter imposes an essential constraint on the amplitude of the forces/torques in the FEA-modeling that must ensure reasonable deflections. Another conclusion concerns comparison of the SVD-based and LIN-based methods. It justifies advantages of the proposed LIN-based technique that provides the best robustness and lower computational complexity.

Influence of FEA-modeling errors

By its general principle, the FEA-modeling is an approximate method that produces some errors caused by the discretization. Beside, even for the perfect modeling, the deflections in the neighborhood of the reference point do not exactly obey the equation [START_REF] Duffy | Statics and Kinematics with Applications to Robotics[END_REF]. Hence, it is reasonable to assume that the "rigid transformation" (5) incorporates some additive random errors

( ) ; 1, δ i i i i in        p p R φ p p ε (21)
that are supposed to be independent and identically distributed Gaussian random variables with zero-mean and standard deviation .

In the frame of this assumption, the expression for the deflections (18) can be rewritten as

1 o1 ο 1 1 1 ˆδδ ; n n n TT i i i i i i i i n                pp ε φ φ P P P ε (22)
where the superscript "o" corresponds to the "true" parameter value. This justifies usual properties of the adopted point-type estimator [START_REF] Martin | Introduction to matrix methods of structural analysis[END_REF], which is obviously unbiased and consistent. Furthermore, the variance-covariance matrices for t,  may be expressed as

    1 2 2 1 ˆĉov ; cov δ n T ii i n          pI φ P P ( 23
)
allowing to evaluate the estimation accuracy using common confidence interval techniques.

As it follows from ( 21), for the translational deflection p the identification accuracy is defined by the standard deviation n  and depends on the number of the points only. In contrast, for the rotational deflection, the spatial location of the points is a very important issue. In particular, for the cubic field of the size aaa , the standard deviation of the rotation angles may be approximately expressed as /6 an  .

To evaluate the standard deviation  describing the random errors ε , one may use the residual-based estimator obtained from the expression

2 2 1 E ( ) (3 6) δ n i i i i n              p p R φ p p . ( 24
)
where (.) E denotes the expectation of the random variable. The latter may be easily derived taking into account that, for each experiment, the deflection field consists of n three-dimensional vectors that are approximated by the model containing 6 scalar parameters. Moreover, to increase accuracy, it is prudent to aggregate the squared residuals for all FEA-experiments and to make relevant estimation using the coefficient

2 (3 6) nm  
, where m is the experiments number.

In addition, to increase accuracy and robustness, it is reasonable to eliminate outliers in the experimental data. They may appear in the FEA-field due to some anomalous causes, such as insufficient meshing of some elements, violation of the boundary conditions in some areas of the mechanical joints, etc. The simplest and reliable method that is adopted in this research is based on the "data filtering" with respect to the residuals (i.e. eliminating certain percentage of the points with the highest residual values). Another practical question is related to detecting zero elements in the compliance matrix or, in other words, evaluating the statistical significance of the computed values compared to zero. These issues will be studied in Section 6.

Optimal settings for virtual experiments

Finally, to evaluate the combined influence of various error sources, let us focus on the deflection identification from the cubic field of size 101010 mm 3 (1331 points, mesh step 1 mm). In particular, let us contaminate all deflections using the Gaussian noise with the s.t.d. 5×10 -5 mm that is a typical value discovered from the examined FEA data sets (Table 4). Similar to the previous case, all calculations were performed using the double precision floating-point arithmetic (16 decimal digits).

The obtained results confirmed the main theoretical derivations of the previous subsections. The identification errors obey the normal distribution (Figure 7) but their s.t.d. should be evaluated taking into account some additional issues. Thus, the s.



mm. The influence of the linearization and round-offs is negligible here (this component is less than 10 -14 mm). Also, this type of the error does not depend on the translation amplitude. In contrast, for the rotational deflections, there exists strong dependence on the amplitude (Figure 8). In particular, for the angular deflection 0.1°, the s.t.d. of the identification error is about 2). Moreover, the simulation results allow to define preferable values of the angular deflection that may be extracted from the FEA-data with the highest accuracy. They show that the deflection angles should be in the range 0.01 …0.2° to ensure the identification accuracy of about 0.2% Figure 8 Identification errors for different amplitudes of the rotational deflections Thus, the case studies considered in Section 5 confirm that the proposed LIN-based algorithm ensures the same accuracy as the SVD-based one, while possessing lower computational complexity. Besides, these results give some practical recommendations for setting the FEA-based experiments and evaluating the identification accuracy.

Minimization of the identification errors

To demonstrate efficiency of the developed technique and to evaluate its applicability to real-life situations, let us consider an example for which the desired compliance matrix can be obtained both numerically and analytically. A comparison of these two solutions provides convenient benchmarks for different FEA-modeling options and also gives some practical recommendations for achieving the required accuracy.

Benchmark example

As an example, let us consider a cantilever beam of size 10001010 mm 3 (see Figure 9) with the Young's Modulus . These data correspond to geometry and material properties of a typical robot link studied in this work.

where non-zero elements are:

2 35 53 /2 y k k L E I    , 11 / k L E A  , 3 22 /3 z k L E I  , 3 33 /3 y k L E I  , 44 / k L G J  , 55 /3 y k L E I  , 66 /3 z k L E I  , 2 26 62 /2 z k k L E I  .
Here L is the length of the beam, A is its cross-section area, y I , z I are the second moments, J is the cross-section property.

During modeling, the loads have been applied at one end of the beam with the other end fully clamped. The force/torque amplitudes were determined using expression [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF] and the optimal accuracy settings for the deflections 0.1…1.0 mm and 0.01…0.20°, which yielded the following values: 1000

x N F  , 1 y FN  , 1 z FN  , 1 • x M N m  , 1 • y M N m , 1 • z M N m 
. These loads were applied sequentially, providing 6 elementary FEA-experiments, each of which produced a single column of the compliance matrix k, in accordance with (4).

Optimal selection of FEA-modeling options

Since the identification errors, which are studied here, essentially depend on the discretization of the FEA-based model and definition of the deflection field, let us focus on the influence of the meshing parameters of the FEA-model and on dimensions of the virtual sensor that specifies this field.

Meshing options. The adopted software provides two basic options for the automatic mesh generation: linear and parabolic ones. It is known that, generally, the linear meshing is faster computationally but less accurate. On the other hand, the parabolic meshing requires more computational resources, but leads to more accurate results.

For the considered case study, both meshing options have been examined and compared with respect to the accuracy of the obtained compliance matrix. The mesh size was gradually reduced from 5 to 1 mm, until achieving the lower limit imposed by the computer memory size. The obtained results (Table 5 and Figure 10) clearly demonstrate advantages of the parabolic mesh, which allow achieving appropriate accuracy of 0.1% for the mesh step of 2 mm using standard computing capacities. In contrast, the best result for the linear mesh is 12% and corresponds to the step of 1 mm. 10% Accuracy of these results can also be evaluated by the difference between the deflections computed using analytical expression [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF] and the identified compliance matrix (assuming that the link is under similar loadings as in virtual experiments). This allows us to simplify physical meaning of the identification errors and also to obtain a weighted performance measure that neglects some identification failures, whose impact on the reference point deflection is not essential. Relevant results are presented in Figure 11, where there are presented the worst values selected from six virtual experiments. They confirm that the mesh 2P (parabolic with step 2 mm) ensures accuracy of 0.1% and is preferable for practice. On the other hand, as it follows from separate study, further reduction of the mesh step does not lead to essential reduction of the identification errors. Obviously, this conclusion is valid for this case study only, but it can be applied to other cases using proper scaling of dimensions. Defining the virtual sensor. The developed technique operates with the deflection field bounded by the virtual sensor (see Figure 4), which is located in the neighborhood of the reference point (RP). As stated above, this neighborhood should be large enough to neutralize the influence of the FEA-induced errors, but its unreasonable increase may lead to violation of some essential assumptions and, consequently, to a reduction in the accuracy. To get a realistic inference concerning this issue, a number of experiments have been carried out, for different definitions of the RP-neighborhood (i.e., the virtual sensor size). The obtained results (Table 6) show that the highest accuracy 0.1% is achieved for the one-layer configuration of the deflection field, which is composed of the nodes located on the rare edge of the examined beam. This pattern is very close to the square-type field 1010 mm 2 studied above (see Section 4.4). In contrast, increasing the sensor size up to the cubic-type of 101010 mm 3 leads to the identification error of about 0.08%. Hence, in practice, it is reasonable to estimate the deflection values from the field of points corresponding to the square-type virtual sensor.

It is worth mentioning that, in practice, the reference point of the link may be located outside the the link material (see Figure 2). Consequently, the RP-neighborhood does not include any finite elements that may be used for creating the deflection field required by the identification procedure. In this case, the link CAD-model should be complemented with an additional solid body centered in the reference point and restrained by the surfaces of the corresponding joint (this body should be rigid enough to insure correctness of FEA-based simulations). After such modifications, the virtual sensor can be defined in the usual way. As it follows from our experience, proper definition of the virtual sensor (and additional solid body) as well as applied loading (distributed on the joint surfaces) play a crucial role in conducting of virtual experiments.

Statistical processing of FEA-based data

Though the finite element analysis is based on strictly deterministic assumptions, it includes tedious computations that may generate some errors, which may be treated statistically. This idea is applied below in order to improve the identification accuracy and to detect the stiffness matrix elements that may be set to zero (in practice, the stiffness matrices with strictly zero elements are rather common, see eq. ( 25)).

Eliminating outliers. As noticed above, the FEA-modeling data may include some anomalous samples that do not obey the assumed statistical properties. This phenomenon has been detected in 2 of 6 experiments, (see Figure 12) where the histograms demonstrated obvious presence of the outliers changing the regular distribution shape (local maximums around the tails). For this reason, a straightforward filtering technique was applied that eliminated 10% the nodes corresponding to the highest residual values (see Figure 13). This technique essentially improved the identification accuracy, the maximum error for the compliance matrix elements reduced from 0.1% to 0.05%.

Figure 12 Residuals for stiffness model identification with parabolic mesh of 2 mm

It is worth mentioning that here, because of the 3-dimensional nature of the problem, each node has been evaluated by three residuals and it was eliminated if any of these residuals was treated as an outlier. In addition, the detailed analysis showed that the outliers were concentrated at the beam edges, which confirms previous assumptions concerning the FEA-induced errors.

Eliminating non-significant elements. According to [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF], the desired compliance matrix includes a number of zero elements (26 of 36), but the proposed identification procedure may produce some small non-zero values. To evaluate their statistical significance, for each element k ij the confidence interval was computed. Relevant computations were performed using expressions for the variances of the deflections (21) and the s.t.d. value of the FEA-modeling noise, which was estimated as (by averaging for all 6 experiments). Then, the computed confidence intervals were scaled in accordance with (4), to be adopted to corresponding elements of the matrix k (this procedure involves simple division by the magnitude of the applied force or torque).

Using this approach, the compliance matrix was revised by assigning to zero all non-significant elements. The employed "decision support" algorithm treated an element as a non-significant one if its confidence interval included zero. As shown in Table 7, this technique allowed us to detect all 26 zero elements mentioned above. It should be also noted that all non-zero elements were evaluated as "significant" ones, with essential "safety factor" of 10 2 and higher. Remarks and comments. The presented illustrative example that deals with a classical element (cantilever beam) confirmed validity of the developed method but also demonstrated some limitations of the FEA-modeling with respect to the stiffness analysis. In particular, some (not very essential but non negligible) disagreement between numerical values of the applied forces/torques and their values extracted from the modeling protocol were detected. Besides, there are a number of non-trivial issues in defining modeling options that are normally set by default. All these factors contribute to the accuracy, but a practically acceptable level 0.1% can be achieved rather easily, using standard computing facilities. Some additional enhancement can be achieved symmetrizing the obtained matrix : ( ) / 2 T  k k k that is motivated by the physical reasons.

Application example: stiffness matrices for Orthoglide links

Let us apply the proposed methodology to the identification of the link stiffness matrices for a parallel manipulator of Orthoglide family [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF]. The principal components of the manipulator are presented in Figure 14, where the elements (b, e, f) are treated as flexible ones and the element (c) is assumed to be rigid. In previous works, relevant stiffness matrices have been obtained either via single-beam approximation [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF] or by using FEA-based identification procedure with linear meshing option [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. Besides, joints particularities have not been taken into account and their influence on the elements of the stiffness matrix was not studied. Hence, it is quite possible that some stiffness matrix elements have been identified with essential errors. In accordance with Section 6.1, the FEA-based virtual experiments were performed using the parabolic mesh of size 2 mm (2P in Figure 10) and the test forces/torques of 1 N and 1 N•m respectively, which, for considered link size and material properties, correspond to the area of linear force-deflection relation. As it follows from Table 5, such settings ensure accuracy of about 0.1%. Identification results are presented in Table 8, where the obtained matrices include a number of zero elements that were detected using the technique presented in Section 6.2. It is also worth mentioning that these matrices are symmetrical, which confirms validity of the developed method. For comparison purposes, similar matrices have been computed using other methods (single-and multi-beam approximations, the FEA-modeling with linear meshing). They are presented in Table 9 and show essential dissimilarity in the evaluation of some matrix elements by different methods. For instance, for the parallelogram bar (i.e., link (f) in Figure 14), the torsional compliance defined by the element k 44 differs by a factor of 13. The main reason for this is that the developed technique takes into account the joint particularities that define the force distribution rule for the applied loadings (previous results assumed that the loading was localized in the reference point). However, the final conclusion concerning accuracy of the obtained stiffness matrices can be obtained after integration of these matrices in the stiffness model of the entire manipulator and comparing it with a straightforward FEA-modeling of the manipulator assembly.

Hence, presented examples confirm advantages of the developed stiffness matrix identification technique and demonstrate its ability to take into account complex link shapes as well as joint particularities related to the force/torque distribution. In addition, it produces symmetrical matrices with some zero elements in accordance with physical properties of the considered manipulator components. In the following chapters, these matrices will be used for the VJM-based stiffness modeling of parallel manipulators.

Conclusions

The paper proposes a CAD-based approach for identification of the elasto-static parameters of the robotic manipulators. The main contributions are in the areas of virtual experiment planning and algorithmic data processing, which allows to obtain the stiffness matrix with required accuracy. In contrast to previous works, the developed technique operates with the deflection field produced by virtual experiments in a CAD environment. The proposed approach provides high identification accuracy (about 0.1% for the stiffness matrix element) and is able to take into account the real shape of the link, coupling between rotational/translational deflections and joint particularities. To compute the stiffness matrix, the numerical technique has been developed, and some recommendations for optimal settings of the virtual experiments are given. In order to minimize the identification errors, the statistical data processing technique was applied. The advantages of the developed approach have been confirmed by case studies dealing with the links of parallel manipulator of the Orthoglide family, for which the identification errors have been reduced to 0.1%. Further development in this area will focus on automation of some operations that currently are performed in the interactive mode and full integration of the developed technique in existing CAD environment. In particular, it is reasonable to develop a dedicated routines that simplify data exchange between the CAD and the computational procedures as well as definition of the simulation parameters while performing the virtual experiments.

     incorporated in the displacement vector p and the rotation matrix R .

Figure 5

 5 Figure 5 Field of points from FEA modeling

Figure 7

 7 Figure 7 Histograms for the identification errors ( a = 1.0 mm, b=0.1°, = 5×10 -5 mm)

  linearization component is about 8.8•10 -5 deg (see Table

Figure 9

 9 Figure 9 Examined cantilever beamFor this element, an analytical expression for the compliance matrix can be presented as

Figure 10 Figure 11

 1011 Figure 10 Errors (%) in estimations of non-zero elements of the compliance matrix (3L,2L,1L are linear mesh with steps 3,2,1 mm; 5P,3P,2P are parabolic mesh with steps 5,3,2 mm)

Figure 13

 13 Figure 13 Filtering of the deflection field outliers

Figure 14

 14 Figure 14 CAD model of Orthoglide and its principal components: Orthoglide (a), foot (b), end-effector (c), parallelogram (d), parallelogram axis (e) and bar (f)

Table 1

 1 Evaluation of the rotation angles from matrix R

	R	1 δ  1   δ 1 zy zx yx        

Table 2 ,

 2 Table 3. All calculations have been performed using the double precision floating-point arithmetic.

	with the

Table 2

 2 Identification errors for the rotation  [deg]

	Method

Table 4

 4 Parameters of the FEA-modeling noise for different mesh type

	Meshing options

Table 5

 5 Maximum errors in estimation of compliance matrix elements

		Meshing options		Maximum errors in elements of
	Type of the mesh	Abbreviation	Size of the finite element	identified matrix k ij
		3L	3 mm	27%
	Linear mesh	2L	2 mm	20%
		1L	1 mm	12%
		5P	5 mm	3.3%
	Parabolic mesh	3P	3 mm	0.19%
		2P	2 mm	0.

Table 6

 6 Evolution of the identification errors for different virtual sensors

	h-sensor	k 11	k 22	k 33	k 44	k 55	k 66	k 62	k 53
	1-layer	0.10%	0.05%	0.05%	0.03%	0.10%	0.10%	0.04%	0.04%
	0.1%	0.11%	0.13%	0.13%	0.01%	0.08%	0.09%	0.09%	0.09%
	0.2%	0.09%	0.2%	0.2%	0.05%	0.04%	0.04%	0.14%	0.14%
	0.3%	0.04%	0.28%	0.28%	0.10%	0.02%	0.02%	0.19%	0.19%
	0.5%	0.06%	0.42%	0.42%	0.20%	0.13%	0.13%	0.29%	0.29%
	1%	0.33%	0.8%	0.80%	0.46%	0.41%	0.41%	0.54%	0.54%
	2%	0.86%	1.55%	1.55%	0.98%	0.91%	0.92%	1.03%	1.03%

Table 7

 7 Eliminating non-significant elements

	Parameter	Estimated value	CI	Exact value
	k 11	5.00•10 -5	±2.3•10 -8	5.00•10 -5
	k 22	2.00•10 0	±2.4•10 -5	2.00•10 0
	k 33	2.00•10 0	±2.4•10 -5	2.00•10 0
	k 44	9.00•10 -6	±6.0•10 -9	9.00•10 -6
	k 55	6.00•10 -6	±8.2•10 -9	6.00•10 -6
	k 66	6.00•10 -6	±8.2•10 -9	6.00•10 -6
	k 26	-3.00•10 -3	±2.4•10 -8	-3.00•10 -3
	k 62	-3.00•10 -3	±8.5•10 -6	-3.00•10 -3
	k 35	3.00•10 -3	±2.4•10 -8	3.00•10 -3
	k 53	3.00•10 -3	±8.5•10 -6	3.00•10 -3
	k 31	1.8•10 -8	±2.3•10 -8	0
	k 32	7.8•10 -6	±2.4•10 -5	0
	k 23	7.8•10 -6	±2.4•10 -5	0
	k 45	-6.2•10 -10	±8.4•10 -8	0
	k 54	-5.7•10-10	±5.6•10 -8	0
	k 63	-6.7•10 -9	±2.4•10 -8	0

Table 9

 9 Comparison of the compliance matrix elements obtained from different methods

				Compliance matrix elements
	Method	k 11	k 22	k 33	k 44	k 55	k 66
		[mm/N]	[mm/N]	[mm/N]	[rad/N•mm]	[rad/N•mm]	[rad/N•mm]
			Link (b): foot
	Single-beam approximation	3.45•10 -4	3.45•10 -4	18.1•10 -4	2.10•10 -7	2.10•10 -7	0.91•10 -7
	Four-beam approximation	2.77•10 -4	4.34•10 -4	17.9•10 -4	2.11•10 -7	1.95•10 -7	0.91•10 -7
	FEA-based evaluation (linear mesh)	2.45•10 -4	3.24•10 -4	15.9•10 -4	2.07•10 -7	2.06•10 -7	1.71•10 -7
	FEA-based evaluation (parabolic mesh)	2.77•10 -4	4.15•10 -4	19.4•10 -4	2.29•10 -7	2.30•10 -7	0.84•10 -7
			Link (e): parallelogram axis
	Single beam approximation	1.34•10 -6	2.65•10 -5	2.65•10 -5	4.29•10 -8	3.18•10 -8	3.18•10 -8
	FEA-based evaluation (linear mesh)	1.99•10 -6	1.29•10 -5	1.50•10 -5	6.81•10 -6	8.23•10 -6	2.67•10 -6
	FEA-based evaluation (parabolic mesh)	6.23•10 -6	2.83•10 -5	2.59•10 -5	2.77•10 -7	4.84•10 -7	1.20•10 -7
			Link (f): parallelogram bar
	Single-beam approximation	3.75•10 -5	4.38•10 -2	1.09•10 -1	3.96•10 -6	3.40•10 -6	1.37•10 -6
	FEA-based evaluation (linear mesh)	4.50•10 -5	3.64•10 -2	8.01•10 -2	3.76•10 -6	2.65•10 -6	1.09•10 -6
	FEA-based evaluation (parabolic mesh)	4.55•10 -5	5.08•10 -2	2.33•10 -1	2.88•10 -5	7.19•10 -6	1.50•10 -6

Acknowledgments

The work presented in this paper was partially funded by the Region "Pays de la Loire", France, by the project ANR COROUSSO, France and FEDER ROBOTEX, France.