S Caro

Projet COROUSSO

Chapter 1 Introduction

Industrial robots were originally dedicated to "pick and place" operations. They start to be used for machining operations such as trimming, deflashing, degating, sanding and sawing. Machining robots are currently developed to help the operator realize machining operations, which have been mainly performed with Computer Numerical Control (CNC) machines so far. Industrial robots can be used, first to reduce scrap rates and production costs, and secondly to increase the volume and flexibility of production lines. However, industrial robots are not as accurate as CNC machine tools, but the total cost of a machining robot is 30% less compared with an equivalent CNC machine. Therefore, it makes sense that the larger the parts to be machined, the more interesting robot machining. Some research works on robotic-based manufacturing have been conducted in the last two decades. The main drawback of industrial serial robots for machining operations turns to be their lack of stiffness and as a result, large manufacturing errors occur [START_REF] Matsuoka | High-Speed End Milling of an Articulated Robot and its Characteristics[END_REF][START_REF] Zhang | Machining with Flexible Manipulator: Toward Improving Robotic Machining Performance[END_REF][START_REF] Abele | Modeling and Identification of an Industrial Robot for Machining Applications[END_REF]. Indeed, a conventional machine tool is about fifty times stiffer than an industrial serial robot [START_REF] Zhang | Machining with Flexible Manipulator: Toward Improving Robotic Machining Performance[END_REF]. Matsuoka et al. [START_REF] Matsuoka | High-Speed End Milling of an Articulated Robot and its Characteristics[END_REF] studied the robotic-based milling. They claimed that the machining quality can be improved by decreasing the cutting forces, namely, by decreasing the diameter of the tool and by increasing the spindle speed. Zhang et al. [START_REF] Zhang | Machining with Flexible Manipulator: Toward Improving Robotic Machining Performance[END_REF] proposed a real-time compensation of the robot end-effector displacement, based on the evaluation of each joint deflection. To this end, they wrote an elastostatic model of the robot at hand and used a wrench sensor to measure the forces and moments applied on the robot end-effector at any time. Besides, Tisius et al. [START_REF] Tisius | An Empirical Approach to Performance Criteria for Manipulation[END_REF] described 50 operational criteria and 50 potential criteria for serial manipulator operation, much of which applies to precision machining. Machining quality and robot energy consumption may also depend on the workpiece placement. Thus, Ur-Rehman et al. focused on the multiobjective path placement optimization of a parallel kinematics machine based on energy consumption, shaking forces and maximum actuator torques [START_REF] Ur-Rehman | Multiobjective Path Placement Optimization of Parallel Kinematics Machines Based on Energy Consumption, Shaking Forces and Maximum Actuators Torques: Application to the Orthoglide[END_REF]. The vibration of the robot itself is another issue to deal with in robotic machining [START_REF] Abele | Modeling and Identification of an Industrial Robot for Machining Applications[END_REF][START_REF] Pan | Chatter Analysis of Robotic Machining Process[END_REF].

Projet COROUSSO Livrable L1.5

As a matter of fact, the first natural frequency of an industrial robot is about 10 to 100 times smaller than the one of conventional machine tools [START_REF] Pan | Chatter Analysis of Robotic Machining Process[END_REF]. Pan et al. also noticed that the vibrations occurring during robotic machining operations are not regenerative but structural [START_REF] Pan | Chatter Analysis of Robotic Machining Process[END_REF]. Therefore, they came up with a stability criterion from the elastodynamic modeling of the robot and the cutting forces. Then, they formulated some rules for robotic machining from several experimental tests. This report introduces a methodology that aims to determine the best placement of the workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the robot. A machining quality criterion is proposed and an optimization problem is formulated and solved. The KUKA KR270-2 robot is used as an illustrative example throughout the report. Section 2 presents the workpiece test and the cutting forces applied by the workpiece on the tool along the machined path. Section 3 describes the robotic cell. Section 4 deals with the formulation of a mono-objective optimization problem to find the optimum placement of the workpiece with regard to a proposed machining quality criterion. Section 5 highlights the optimum and worst workpiece placements within the robotic cell. The best and worst redundancy planning schemes associated with those placements are also determined. Section 6 is about the conclusions of the report and the future work.

Projet COROUSSO Livrable L1.5

Chapter 2

Robotic Cell and Parameterization

The robotic cell shown in Fig. 2.1 is composed of a KUKA KR270-2 robot, a FISCHER milling spindle mounted on the robot end-effector and a rotary table. The geometric parameters of the robot are presented in [START_REF] Dumas | Joint Stiffness Identification of Industrial Serial Robots[END_REF]. An identification procedure was developed in [START_REF] Dumas | Joint Stiffness Identification of Six revolute Industrial Serial Robots[END_REF] in order to determine all joint stiffness values of the KUKA KR270-2 robot. Accordingly, the translational and rotational displacements of the tool can be predicted knowing the wrench exerted on it. It is noteworthy that the flexibility of the links is considered through the identified joint stiffness values. The FISCHER milling spindle is of type MFW 1709/24. Its maximum rotational speed is equal to 24 000 rpm and its power is equal to 20 kW. The orientation of the workpiece can be fully represented with the variables (Q 2 , Q 3 , Q 4 ), a subset of the quaternions coordinates. Indeed, the quaternions represent the rotations of the workpiece with a rotation axis u = [u x u y u z ] T and an angle θ. The relation between the quaternions and the axis and angle representation can be found in [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

Q 1 = cos(θ/2) (2.1) Q 2 = u x sin(θ/2) (2.
2)

Q 3 = u y sin(θ/2) (2.3) Q 4 = u z sin(θ/2) (2.4)
where u 2 x + u 2 y + u 2 z = 1 and 0 ≤ θ ≤ π. Nevertheless, the workpiece can be only rotated about the vertical axis of the robot base frame F 0 due to the geometry of the rotary table. Therefore,

Q 2 = Q 3 = 0 (2.5)
Projet COROUSSO Livrable L1.5

Robot Spindle frame F 0 , the frame F 6 attached to the robot end-effector, the frame F 7 attached to the spindle and the workpiece frame F W . From the previous closed loop chain, we can write:

Rotary table

≤ Q 4 ≤ 1 (2.6) Z 0 Y 0 X 0 F 0 Z 6 Y 6 X 6 F 6 Z 7 Y 7 X 7 F 7 Z W Y W X W F W
0 T W W T P i = 0 T 6 6 T 7 7 T P i (2.7)
where 0 T W is the homogeneous transformation matrix from frame F 0 to frame F W expressed as:

0 T W =            2Q 2 1 -1 -2Q 1 Q 4 0 0 x OW Q 1 Q 4 2Q 2 1 -1 0 0 y OW 0 0 1 0 z OW 0 0 0 1           
(2.8)

0
x OW , 0 y OW and 0 z OW being the Cartesian coordinates of point O W expressed in frame F 0 .

W T P i is the homogeneous transformation matrix from frame F W to frame F P i attached to the ith point of the tool path as shown in Fig. 3.1.

0 T 6 is the homogeneous transformation matrix from frame F 0 to frame F 6 .

6 T 7 is the homogeneous transformation matrix from frame F 6 to frame F 7 and depends on the geometry of the spindle and how the latter is mounted on the robot end-effector. Here,

6 T 7 =            0 - √ 2/2 √ 2/2 0 1 0 0 0 0 √ 2/2 √ 2/2 0.684 0 0 0 1           
(2.9)

7 T P i is the homogeneous transformation matrix from frame F P i to frame F 7 . Note that the KUKA KR270-2 robot has six degrees of freedom, whereas the milling operation sets only five degrees of freedom as the rotation of the spindle about the tool axis is not fixed. Therefore, the robot is redundant with respect to the task and the kinematic redundancy is equal to one as explained in [START_REF] Conkur | Clarifying the definition of redundancy as used in robotics[END_REF][START_REF] Wenger | Curve-following for redundant manipulators with obstacles : feasibility analysis and solutions[END_REF]. Here, the kinematic redundancy is characterized by the angle β i , which corresponds to the rotation angle of the spindle about the tool axis

Projet COROUSSO Livrable L1.5

at the ith point of the tool path as shown in Fig. 4.6. As a consequence,

7 T P i =            cos(β i ) sin(β i ) 0 0 sin(β i ) -cos(β i ) 0 0 0 0 -1 0 0 0 0 1            (2.10)
Projet COROUSSO Livrable L1.5

Chapter 3

Machining Task and Cutting Force Model Frame F P i is attached to the ith point of the tool path, i = 1, . . . , n. X P i is along the feed direction. Z P i is along the tool axis and points toward the robot. The machining quality is affected by the robot deviation due to the cutting forces applied on the tool [START_REF] Dumas | Workpiece Placement optimization of Six-revolute Industrial Serial Robots for Machining Operations[END_REF]. The cutting conditions are given in Tab. 3.1 where f z , a p and a e denote the feed rate, the depth of cut and the width of cut, respectively. The cutting forces are evaluated thanks to the cutting force model described in [START_REF] Tlusty | Dynamics of cutting forces in end milling[END_REF]. As a result, Fig. 3.2 depicts the force components applied by the tool on the workpiece along axes X W , Y W and Z W of frame F W as a function of the tool path point number.

Cutting Force Model

The tool is composed of several teeth. We consider in this report a tool composed of four teeth. applied by each tooth on the workpiece.

F W F P i X W Y W Z W Y P i X P i Z P i 1 2 3 4 5
F x F y F z Point Number [N]
Feed Feed direction direction

ϕ in ϕ out
Piece to mill

F y F x F r F t f z sin(ϕ) f z ϕ Figure 3.

-Chip geometry and definition of ϕ in and ϕ out

The instantaneous forces exerted by the tooth on the workpiece are defined in [START_REF] Tlusty | Dynamics of cutting forces in end milling[END_REF] and are expressed as follows:

δF t = K t a p f z sin(ϕ) dϕ δF r = K r δF t (3.1)
The mean forces F x and F y along the x and y axes associated with one tooth are expressed as:

F x = ϕ in ϕout ( F t cos(ϕ) + F r sin(ϕ)) dϕ F y = ϕ in ϕout (-F t sin(ϕ) + F r cos(ϕ)) dϕ (3.2) Therefore, if ϕ ∈ [ϕ in ϕ out ], F x = K t a p f z (- cos(2ϕ out ) 4 + K r ϕ out 2 -K r sin(2ϕ out ) 4 ) + K t a p f z ( cos(2ϕ in ) 4 -K r ϕ in 2 + K r sin(2ϕ in ) 4
)

F y = K t a p f z ( sin(2ϕ out ) 4 - ϕ out 2 -K r cos(2ϕ out ) 4 ) + K t a p f z (- sin(2ϕ in ) 4 + ϕ in 2 + K r cos(2ϕ in ) 4 ) (3.3)
else,

F x = 0 F y = 0 (3.4)
f z is the feed rate and a p is the depth of cut. Then, the global cutting force exerted by the Projet COROUSSO Livrable L1.5

entire tool on the piece can be calculated. Coefficients K t and K r were determined from cutting experiments with a force measurement system.

Projet COROUSSO Livrable L1.5

Chapter 4

Experimentations

Identification of the Joint Stiffness Values

The translational and rotational displacements of the robot end-effector can be evaluated if the forces and moments applied on it are known. To this end, the stiffness model of the robot defined by Eq. (4.1) was determined based on the CCT Theory [START_REF] Chen | The 6x6 Stiffness Formulation and Transformation of Serial Manipulators via the CCT Theory[END_REF]. The links of the robot are supposed to be quite stiffer than the joints and not known as it happens usually for industrial robots [START_REF] Dumas | Développement de méthodes robotisées pour le parachèvement de pièces métalliques et composites[END_REF]. Then, the relationship between the wrench applied on the robot end-effector and its small displacement screw is defined as follows:

w = K X δX (4.1)
where K X is the Cartesian stiffness matrix of the robot. w is the 6-dimensional vector of forces and moments exerted on the robot end-effector and derived from Sec. 3.2. δX is the obtained 6-dimensional small displacement screw of the end-effector. The Cartesian stiffness matrix K X that depends on the robot configuration, the joint stiffness values and the wrench applied on the end-effector takes the form:

K X = J -T (K θ -K C ) J -1 (4.2)
Projet COROUSSO Livrable L1.5

with

K θ =                   k θ 1 0 0 0 0 0 0 k θ 2 0 0 0 0 0 0 k θ 3 0 0 0 0 0 0 k θ 4 0 0 0 0 0 0 k θ 5 0 0 0 0 0 0 k θ 6                   (4.3)
and

K C = ∂J T ∂θ 1 w ∂J T ∂θ 2 w ∂J T ∂θ 3 w ∂J T ∂θ 4 w ∂J T ∂θ 5 w ∂J T ∂θ 6 w (4.4)
J is the kinematic Jacobian matrix of the robot. K θ is its diagonal joint stiffness matrix and K C is the complementary stiffness matrix, which depends on the wrench applied on the end-effector. An identification procedure was developed in [START_REF] Dumas | Joint Stiffness Identification of Six revolute Industrial Serial Robots[END_REF] in order to determine all joint stiffness values of the robot. Accordingly, the translational and rotational displacements of the end-effector can be predicted knowing the wrench exerted on it. It is noteworthy that the flexibility of the links is considered through the identified joint stiffness values. The joint stiffness values of the KUKA KR270-2 robot are given in Table 4.1. They were obtained by using the experimental setup shown in Fig. 4.1 and described in [START_REF] Dumas | Joint Stiffness Identification of Industrial Serial Robots[END_REF][START_REF] Dumas | A Methodology for Joint Stiffness Identification of Serial Robots[END_REF].

Influence of the Workpiece Placement and the Kinematic

Redundancy on the Machining Quality

This section aims to show the influence of the workpiece placement and the kinematic redundancy on the machining quality of a test workpiece through some experimental results.

Evaluation of the tool displacement during a milling operation

Knowing the joint stiffness values of the robot from Sec. 4.1 and the forces exerted on the robot end-effector from Sec. 3.2, it is possible to predict the tool displacement using the stiffness Projet COROUSSO Livrable L1.5 matrix defined with Eq. (4.1). Four machining operations were realized in order to validate the theoretical elastostatic model of the KUKA KR270-2 robot: (i) a workpiece made up of a resin material called LAB was milled along a 300 mm linear path. It is noteworthy that the LAB material exerts a small force on the tool; (ii) a workpiece made up of an Aluminium alloy was milled along a 300 mm linear path with the cutting conditions P3 given in Table 4.2; (iii) a workpiece made up of an Aluminium alloy was milled along a 300 mm linear path with the cutting conditions P4 given in Table 4.2; (iv) a workpiece made up of an Aluminium alloy was milled along a 300 mm linear path with the cutting conditions P5 given in Table 4.2. N , a p and V f denote the rotational speed of the tool, the depth of cut and the feed speed, respectively. paths for the other machining operations is about 5% and is also partly due to the inherent straightness of the robot. Consequently, it is possible to evaluate the path that will be followed by the tool considering its displacement due to the cutting forces from the elastostatic model of the robot defined with Eq. (4.1).

Influence of the task placement on the machining quality

Influence of the kinematic redundancy on the machining quality

The KUKA KR270-2 robot has six degrees of freedom. The milling operations sets only five degrees of freedom because the rotation of the spindle about the tool axis is not fixed. Therefore, the robot is redundant with respect to the task, the robot redundancy being equal to one as explained in [START_REF] Conkur | Clarifying the definition of redundancy as used in robotics[END_REF][START_REF] Wenger | Curve-following for redundant manipulators with obstacles : feasibility analysis and solutions[END_REF]. Here, the robot redundancy is characterized by the angle β, which denotes the rotation angle of the spindle about the tool axis and is depicted in Fig. 4.6. Figure 4.7 illustrates the quality of two test workpieces made up of an aluminium alloy and milled in the same placement, but with two different redundancy planning schemes. It is apparent that the milling quality also depends on the redundancy planning scheme chosen for the machining operation.

Projet COROUSSO Livrable L1.5 Chapter 5

β i Spindle Z 7 X P i

Optimum Workpiece Placement

This section deals with a methodology that aims to determine the best placement of the workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the robot. The proposed methodology is highlighted through the definition of a machining quality criterion and the formulation of a mono-objective optimization problem.

Machining quality criterion

It makes sense that the tool displacement c i , normal to both the feed direction along axis X P i and the tool axis Z 7 , is mainly responsible for machining errors. 0 x P i and 0 z 7 are the unit vectors of axes X P i and Z 7 expressed in frame F 0 , respectively. c i can be expressed as follows:

c i = 0 δd T i 0 z 7 × 0 x P i (5.1)
where 0 δd i is the point-displacement of the tool at the ith point of the tool path expressed in frame F 0 and × denotes the cross product. The point-displacement 0 δd i is evaluated thanks to the Cartesian stiffness matrix K of the KUKA KR270-2 robot obtained in [START_REF] Dumas | Joint Stiffness Identification of Industrial Serial Robots[END_REF] and from the wrench w applied by the tool on the workpiece that are depicted in Fig. 3.2, namely,

0 δd i = K -1 w (5.2)
The machining quality criterion f M QC , which is the objective function of the optimization problem at hand, is defined as the mean of c i values, namely,

f M QC = 1 n n i=1 c i → min (5.
3)

The smaller f M QC , the better the machining quality.

Projet COROUSSO Livrable L1.5

Formulation of the optimization problem

From (5.3) to (5.16), the optimization problem to solve in order to find the best workpiece placement can be formulated as follows:

minimize

f M QC over x = 0 x OW 0 y OW 0 z OW Q 4 β T µ T T subject to max (|c 1 |, |c 2 |, • • • , |c n |) ≤ c max -185 deg ≤ θ 1i ≤ 185 deg 0 ≤ θ 2i ≤ 146 deg -245 deg ≤ θ 3i ≤ 29 deg -350 deg ≤ θ 4i ≤ 350 deg -125 deg ≤ θ 5i ≤ 125 deg -350 deg ≤ θ 6i ≤ 350 deg 1/κ F (J N i ) > 0.2 x min ≤ 0 x OW ≤ x max y min ≤ 0 y OW ≤ y max z min ≤ 0 z OW ≤ z max -1 ≤ Q 4 ≤ 1 i = 1, . . . , n (5.17)
Optimization problem (5.17) aims to find the optimum workpiece placement and the corresponding optimum redundancy planning scheme that minimize the machining quality criterion f M QC defined by (5.3) while respecting the set of constraints.

Result Analysis

A hybrid optimization algorithm was used to solve optimization problem (5.17). As a matter of fact, a genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search, Optimization & Machine Learning[END_REF] was used to find the initial guess of an interior-point algorithm for large-scale nonlinear programming [START_REF] Byrd | An Interior Point Algorithm for Large-Scale Nonlinear Programming[END_REF]. The Matlab ga function combined with the Matlab fmincon function were used.

The genetic algorithm converged after 51 generations as shown in The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement Figure 5.2 illustrates the tool displacement c i with respect to the tool path point number and β angle for the optimal workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement. The white areas can not be reached by the robot because of its joint limits. It is noteworthy that 0 x OW , 0 y OW , 0 z OW , Q 4 and µ are the only decision variables considered by the ga and fmincon functions in this optimization problem solving. As a matter of fact, an optimal β vector is searched at each iteration of the genetic algorithm and at each iteration of the interior-point algorithm. This vector is obtained in such a way that it minimizes the objective function f M QC and avoids discontinuities in the robot joint space along each segment.

Projet COROUSSO Livrable L1.5 It is noteworthy that the worst workpiece placement is obtained by maximizing the objective function f M QC while respecting the constraints of optimization problem (5.17 An video of the milling operation corresponding to Case 1 can be downloaded at 1 . An video of the milling operation corresponding to Case 2 can be downloaded at 2 . An video of the milling operation corresponding to Case 3 can be downloaded at 3 . An video of the milling operation corresponding to Case 4 can be downloaded at 4 5 . 3) for the previous four milling conditions. It is apparent that both the workpiece placement and the kinematic redundancy affect the machining quality of the workpiece. Table 5.1 shows that it is important to pay attention to the workpiece placement and kinematic redundancy. Indeed, f M QC value for Case 4 is more than 14 times higher than f M QC value for Case 1.

Figure 2 . 1 -

 21 Figure 2.1 -Robotic cell located at University of Nantes, France

Figure 2 . 2 -

 22 Figure 2.2 -Closed loop chain and frames F 0 , F 6 , F 7 and F W

Figure 2 . 3 -

 23 Figure 2.3 -Kinematic redundancy characterized with angle β i

Figure 2 .

 2 Figure 2.2 represents the closed loop chain composed of the KUKA KR270-2 robot, the FISCHER milling spindle, the workpiece and the rotary table. It also depicts the robot base

3. 1 Figure 3 . 1

 131 Figure 3.1 illustrates the workpiece to be machined, which is made up of aluminum alloy. F W of origin O W is the frame attached to the workpiece. The five segments AB, BC, CO W , O W D and DE, of length equal to 200 mm each, have to be milled. The tool path is offset by the tool radius from the five segments to be milled. The tool path is discretized into n points and is shown in dashed line in Fig.3.1. Frame F P i is attached to the ith point of the tool path, i = 1, . . . , n. X P i is along the feed direction. Z P i is along the tool axis and points toward the robot.

Figure 3 .

 3 3 presents the chip and the directions of the tangential force F t and the radial force F r

Figure 3

 3 Figure 3.1 -Workpiece

Figure 3 . 2 -

 32 Figure 3.2 -Force components applied by the tool on the workpiece expressed in frame F W as a function of the tool path point number

Figure 4 . 2 -

 42 Figure 4.2 -Calculated and measured paths associated with four milling conditions

Figure 4 . 3 -

 43 Figure 4.3 -Test workpiece

Figure 4 .

 4 [START_REF] Ur-Rehman | Multiobjective Path Placement Optimization of Parallel Kinematics Machines Based on Energy Consumption, Shaking Forces and Maximum Actuators Torques: Application to the Orthoglide[END_REF] shows the measures of the four resulting workpieces, namely, the milling quality of the four resulting workpieces. It is apparent that the resulting workpieces made up of a LAB material are closed to the desired workpiece represented with the black straight line in Fig.4.5, while there is a non negligible gap betweeen the resulting workpieces made up of an Aluminium alloy and the desired workpiece. Moreover, the difference between the two

Figure 4 . 5 -

 45 Figure 4.5 -Measurement of the test workpiece milled in two different placements

Figure 4 . 6 -

 46 Figure 4.6 -Kinematic redundancy characterized with angle β

Figure 4 . 7 -

 47 Figure 4.7 -Quality of the test workpiece milled in a given placement, but with two different redundancy planning schemes

5 T

 5 Fig. 5.1, each population containing 120 individuals. Then, a local optimum decision variable vector x opt of optimization Projet COROUSSO Livrable L1.5 problem (5.17) was obtained with the Matlab fmincon function and is expressed as: x opt = 0.0022 -1.6435 -0.5564 0.9258 β opt µ opt T (5.18) with µ opt = 4 4 4 5 meaning that the robot posture changes between segment 3 and segment 4 only.

Figure 5 . 2 -

 52 Figure 5.2 -Tool displacement c i [m] as a function of the tool path point number and β angle for the best workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement

Figure 5 .

 5 Figure 5.3 represents the optimum workpiece placement found by solving optimization problem (5.17).It is noteworthy that the worst workpiece placement is obtained by maximizing the objective function f M QC while respecting the constraints of optimization problem (5.17).

  Figure 5.3 represents the optimum workpiece placement found by solving optimization problem (5.17).It is noteworthy that the worst workpiece placement is obtained by maximizing the objective function f M QC while respecting the constraints of optimization problem (5.17).

Figure 5 . 4 -

 54 Figure 5.4 -Tool displacement c i [m] as a function of the tool path point number and β angle for the worst workpiece placement. The green curve characterizes the optimum redundancy planning scheme, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement

Figure 5 .Case 1 :Case 2 :Case 3 :Case 4 :

 51234 Figure 5.4 illustrates the tool displacement c i with respect to the tool path point number and β angle for the worst workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement. Let us compare the machining quality of the workpiece obtained with the four following milling conditions: Case 1: Optimum workpiece placement with the best redundancy planning scheme (see Fig. 5.2) Case 2: Optimum workpiece placement with the worst redundancy planning scheme (see Fig. 5.2) Case 3: Worst workpiece placement with the best redundancy planning scheme (see Fig. 5.4)

Table 3 .

 3

		1 -Cutting Conditions
	Spindle speed Feed speed	f z	a p	a e
	20 000 rpm	4 m/min 0.05 mm/tooth/rev 5 mm varies

Table 4 .

 4

	k θ 1	k θ 2	k θ 3	k θ 4	k θ 5	k θ 6
	0.237 3.32 2.79 0.486 0.521 0.38

1 -Apparent joint stiffness values in [MNm/rad] of the KUKA KR270-2 robot

Table 4 .

 4

		2 -Cutting conditions
	Reference N [tr/min] a p [mm] V f [mm/min]
	P3	20 000	5	3 600
	P4	10 000	5	3 600
	P5	6 670	5	3 600

Table 5 .

 5 1 -Machining Quality Criterion for Four Milling Conditions Case 1 Case 2 Case 3 Case 4

	f M QC [mm] 0.083	0.35	0.82	1.2

Table 5 .

 5 1 gives the values of the maching quality criterion f M QC defined in(5.

Author: S. CARO

Decision variables

The decision variable vector of the optimization problem at hand contains the Cartesian coordinates of point O W expressed in frame F 0 , i.e., 0 x OW , 0 y OW and 0 z OW and the orientation parameter Q 4 of the workpiece. The decision variable vector also contains the redundancy planning scheme. Let β i be the rotation angle of the spindle about the tool axis at the ith point of the discretized tool path as shown in Fig. 4.6. Accordingly, the redundancy planning scheme is defined by vector β expressed as:

The KUKA KR270-2 robot may have eight solutions to the inverse geometric model, namely, eight postures. The robot has to keep the same posture along a given segment in order to avoid signs of wear. However, the robot posture can change from one segment to the next one. µ j denotes the solution number to the inverse geometric model of the robot along the jth segment of the tool path, j = 1, . . . , 5. As a consequence, the posture scheme is defined by vector µ expressed as:

Constraints

Four types of constraints arise in the optimization problem at hand:

1. The joint angles of the KUKA KR270-2 robot should be bounded between its joint limits, θ ji being the jth joint angle of the robot, j = 1, . . . , 6, when the tool is located at the ith point of the tool path, i = 1, . . . , n. Note that the robot was parameterized by using the modified Denavit-Hartenberg convention [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF].

Projet COROUSSO Livrable L1.5

2. The tool displacement c i expressed with (5.1) should be smaller than a given value c max , which is defined with regard to the expected machining quality, i.e.,

where |.| denotes the absolute value.

3. The workpiece should be located within a work volume that depends on the rotary table size and location within the robotic cell. Here,

y min ≤ 0 y OW ≤ y max (5.14) 4. The robot should be far from singularities while the tool follows the path, namely,

where J N i is the normalized kinematic Jacobian matrix of the KUKA KR270-2 robot defined in [START_REF] Dumas | Joint Stiffness Identification of Industrial Serial Robots[END_REF] and evaluated at the ith point of the tool path. κ F (J N i ) is condition number of matrix J N i based on the Frobenius norm. Author: S. CARO Chapter 6

Conclusions and Future Work

This report introduced a methodology to determine the optimum placement of any workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the KUKA KR270-2 robot. First, a workpiece test was introduced and the forces applied by the workpiece to the tool were evaluated by using an existing cutting force model. Then, the robotic cell was presented. It contains a KUKA KR270-2 robot, a FISCHER milling spindle and a rotary table. The parameterization of the closed loop chain composed of the foregoing components and the workpiece was also described. A machining quality criterion was introduced and a mono-objective optimization problem was formulated in order to determine the optimum placement of a workpiece to be machined and the associated best redundancy planning scheme.

The worst workpiece placement and the associated worst redundancy planning were also obtained by maximizing the objective function of the optimization problem at hand. It is noteworthy that the machining quality criterion corresponding to the optimum workpiece placement and the associated best redundancy planning scheme is more than 14 times smaller than the machining quality criterion corresponding to the worst workpiece placement and the associated worst redundancy planning scheme. Therefore, it is important to pay attention to the workpiece placement and kinematic redundancy for machining operations with a KUKA KR270-2 robot. The theoretical results presented in this report will be checked experimentally soon. The extension of the proposed methodology to robot machining operations with a higher kinematic redundancy of the robot with regard to the task, such as a machining robot mounted on a rail or a mobile platform, are also part of the future work.

Projet COROUSSO Livrable L1.5

Projet COROUSSO Livrable L1.5


Projet COROUSSO Livrable n°1.2 Outils pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO

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 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 | Stiffness 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 onedimensional 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 onedimensional 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 ANR-10-SEGI-003-LI1.

COROUSSO

(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. 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 VJMmodeling 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 recomputing 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.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. ANR-10-SEGI-003-LI1.

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.1), even a four-beam approximation provides an accuracy of only 30-50% [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. 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 2.2. 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 Chapter, 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. Hence, the goal of this Chapter 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 Chapter. 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 δ

⎡ ⎤ ⎡ ⎤ = ⋅ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ F p K M φ (1.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 FEAbased techniques, which hypothetically produce rather accurate result. . .

m m m m - ⎡ ⎤ ⎡ ⎤ = ⋅ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ F F p p k M M φ φ (1.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 x z T T T T y y T T T T y z F M F M F M = = = = = = = = = = = = F M F M F M F M F M F M . (1.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 φ φ φ φ φ φ (1.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 (1.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. 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, } i i i n Δ = p p
(see Figure 4.1) 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 i n + Δ = ⋅ + = p p R φ p p , (1.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 (1.5) with respect to the six scalar variables ( , , , δ , , ) δ δ

x y z x y z p p p ϕ ϕ ϕ incorporated in the displacement vector p and the rotation matrix R . 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

SVD-based solution

For the comparison purposes, first let us briefly present the exact solution of the optimization problem (1.6). 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 (1.6) to the standard form. First, for the nonconstrained variable p, straightforward differentiation and equating to zero gives an expression

1 1 1 ( ) n n i i i i n = = ⎛ ⎞ = Δ -- ⎜ ⎟ ⎝ ⎠ ∑ ∑ p p R I p (1.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 , (1.8)
the original optimization problem is reduced to the orthogonal Procrustes formulation that is expressed via the singular value decomposition (SVD) of the matrix

1 • n T T i i i= = ∑ p g U Σ V , (1.11)
which requires some computational efforts. Hence, the above expressions (1.7), (1.10) allow to solve the optimization problem (1.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

δ δ δ δ 1 1 δ 1 δ z y z x y x ϕ ϕ ϕ ϕ ϕ ϕ - ⎡ ⎤ ⎢ ⎥ ≅ - ⎢ ⎥ - ⎢ ⎥ ⎣ ⎦ R (1.
[ ] ij r = R
in several ways (Table 4.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 3.5.2.

COROUSSO

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 (1.12). This allows rewriteing the equation of the 'rigid transformation' (1.5) in the form δ ; 1,

i i i n Δ = × + = p p φ p (1.13)
that can be further transformed into a linear system of the following form [ ]

; 1 δ , i i i n ⎡ ⎤ = Δ = ⎢ ⎥ ⎣ ⎦ p I P p φ (1.14)
where i P is a skew-symmetric matrix corresponding to the vector :

i p 0 0 0 zi y i i z i x i yi xi p p p p p p - ⎡ ⎤ ⎢ ⎥ = - ⎢ ⎥ - ⎢ ⎥ ⎣ ⎦ P (1.15)
Then, applying the standard least-square technique with the objective

2 1 mi δ n n i i , i f = = Δ - - → ∑ φ t p P φ p (1.16)
one can get the solution

1 1 1 1 1 1 δ n n i i i i n n n T T T i i i i i i i i n - = = = = = ⎡ ⎤ ⎡ ⎤ Δ ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ Δ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ∑ ∑ ∑ ∑ ∑ I P p p φ P P P P p (1.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 c i i n - = = ∑ p p leading to expression 1 1 1 1 1 ˆδ ˆn i i n n T T i i i i i i n - - = = = ⎡ ⎤ ⎡ ⎤ Δ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⋅ ⎛ ⎞ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎜ ⎟ Δ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦⎢ ⎥ ⎣ ⎦ ∑ ∑ ∑ I 0 p p φ 0 P P P p (1.18)
that requires inversion of the matrix of size 3×3. Here, ˆi P denotes a skew-symmetric matrix corresponding to the vector ˆi Case 1: Symmetrical field (Figure 5.1a). If the field is symmetrical with respect to its centre c p , the solution (1.18) can be presented in a compact analytical form as

1 1 n i i n - = = Δ ∑ p p ; 1 1 δ n T i i i - = = Δ ∑ φ D P p (1.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 (1.20)
is diagonal and easily inverted. 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 can 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.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 firstorder 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' (1.5) with the parameters ( , , ) T a a a = p and ( , ,

δ ) T b b b = φ
presented in Table 6.1, Table 6.2. All calculations have been performed using the double precision floating-point arithmetic.

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.

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 (1.5). Hence, it is reasonable to assume that the 'rigid transformation' (1.5) incorporates some additive random errors

( ) ; 1, δ i i i i i n + Δ = ⋅ + + = p p R φ p p ε (1.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 (1.18) can be rewritten as

1 o 1 ο 1 1 1 ˆδ δ ; n n n T T i i i i i i i i n - - = = = ⎛ ⎞ = + = + ⎜ ⎟ ⎝ ⎠ ∑ ∑ ∑ p p ε φ φ P P P ε (1.22)
where the superscript 'o' corresponds to the 'true' parameter value. This justifies usual properties of the adopted point-type estimator (1.18), which is obviously unbiased and consistent. Furthermore, the variancecovariance matrices for t, δϕ may be expressed as allowing to evaluate the estimation accuracy using common confidence interval techniques.

As it follows from (1.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 a a a × × , the standard deviation of the rotation angles may be approximately expressed as / 6 a n σ .

To evaluate the standard deviation σ describing the random errors ε , one may use the residual-based estimator obtained from the expression .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 2 1 E ( ) ( 3 6 ) δ n i i i i n σ = ⎛ ⎞ + Δ - ⋅ - = - ⎜ ⎟ ⎝ ⎠ ∑ p p R φ p p . ( 1
2 (3 6) n mσ -
, 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 3.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 6.3). 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 6.1) but their s.t.d. should be evaluated taking into account some additional issues. Thus, the s.t.d. of the translational error is about

COROUSSO

In contrast, for the rotational deflections, there exists strong dependence on the amplitude (Figure 6.2). In particular, for the angular deflection 0.1°, the s.t.d. of the identification error is about 6.1). 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% Auteurs : AK, AP COROUSSO

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 7.1) with the Young's Modulus . These data correspond to geometry and material properties of a typical robot link studied in this work.

For this element, an analytical expression for the compliance matrix can be presented as

k k k k k k k k k k ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ k (1.25)
where non-zero elements are:

2 35 53 / 2 y k k L EI = =- , 11 / k L EA = , 3 22 / 3 z k L EI = , 3 33 / 3 y k L EI = , 44 / k L GJ = , 55 / 3 y k L EI = , 66 / 3 z k L EI = , 2 26 62 / 2 z k k L EI = = .
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 (1.23) 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 F N = , 1 z F N = , 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 (1.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 7.1 and Figure 7.2) 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. Accuracy of these results can also be evaluated by the difference between the deflections computed using analytical expression (1.25) 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 7.3, 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 3.1), 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.

Figure 7.3 Influence of the stiffness matrix errors on the position 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 7.2) 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.

ANR-10-SEGI-003-LI1. 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 7.2) 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.1). 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. (1.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 7.4) 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 7.5). This technique essentially improved the identification accuracy, the maximum error for the compliance matrix elements reduced from 0.1% to 0.05%. 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 (1.23), the desired compliance matrix includes a number of zero elements (26 of 36), but the proposed identification procedure may produce some small nonzero 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 (1.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 (1.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.3, 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 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 8.1, 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 7.2) 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 7.1, such settings ensure accuracy of about 0.1%. Identification results are presented in Table 8.1, where the obtained matrices include a number of zero elements that were detected using the technique presented in Section 3.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 multibeam approximations, the FEA-modeling with linear meshing). They are presented in 0 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 8.1), 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. ANR-10-SEGI-003-LI1.

CONCLUSIONS

This Chapter 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.

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO Page 5/32 Auteurs : AK, AP COROUSSO 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.

Figure 1 . 1

 11 Figure 1.1 Integration of VJM modeling approach with FEA-based identification technique of the stiffness model parameters

COROUSSO

Figure 2 . 1

 21 Figure 2.1 CAD model of the link with complex shape

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

Figure 3 . 1 2 24

 312 Figure 3.1 Identification experiment for the element with a complex shape

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

Figure 4 . 1

 41 Figure 4.1 Field of points from FEA modeling

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO Page 12/32 Auteurs : AK, AP COROUSSO 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.

  the constraint T = R R I . The latter yields the solution[START_REF] Gower | Procrustes Problems[END_REF]

2 24 5 . 1

 251 that the obtained expression (1.18) is more computationally efficient compared to (1.7) and (1.10). Besides, for some typical cases corresponding to regular meshing patterns of the FEA, expression (1.18) can be further simplified as shown in the following subsection. ANR-10-SEGI-003-LI1.pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de Analytical expressions for typical case studiesIn 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 analytically and to obtain very simple expressions for p, δϕ. Let us consider several patterns that are useful for practice and are presented in Figure5.1.

Case 2 :Case 3 :Figure 5 . 1 2 24

 23512 Figure 5.1 Typical patterns of the deflection field

2 24

 2 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. ANR-10-SEGI-003-LI1.pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

ANR- 10

 10 pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

6 1 . 2 24

 12 36 10 - ⋅ mm and depends only on the FEA-induced component that is evaluated as 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. ANR-10-SEGI-003-LI1.pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO Page 19/32

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

Figure 6 . 1 Figure 6 . 2 2 24 2 24

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

Figure 7 . 1

 71 Figure 7.1 Examined cantilever beam

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

32 COROUSSOFigure 7 . 2

 3272 Figure 7.2 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)

ANRFigure 7 . 4

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

Figure 7 . 5 2 24

 752 Figure 7.5 Filtering of the deflection field outliers

Figure 8 . 1

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

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de

COROUSSO

ANR

  pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO Page 31/32

Real shape Material properties Joints particularities Couplings between translational\rotational deflections Deflections field Stiffness matrix

		Projet COROUSSO Livrable n°1.2	ANR-10-SEGI-003-LI1.2
		Outils pour l'identification des paramètres
		de raideur des robots à l'aide d'un logiciel	24/02/2012
	COROUSSO	de CAO	indice A
			Page 8/32
		CAD model
		FEA-based
		virtual experiments
		Stiffness matrix
		identification
		Figure 2.2 Algorithm for stiffness matrix identification procedure

Table 4

 4

	.1	Evaluation of the rotation angles from matrix R
	Method				ϕ	x								ϕ	y	ϕ	z
	SVD+				32 r								13 r	21 r
	SVD-				-	23 r								-	31 r	-	12 r
	SVD±	(	32 r	-	23 r	)/2		(	13 r r 31 -	)/2	(	21 r r 12 -	)/2
	SVD+asin		32 asin r				13 asin r	21 asin r
	SVD-asin		-	23 asin r				-	31 asin r	-	12 asin r
	SVD±asin	( asin (	32 r	-	23 r	) / 2	)	( asin (	13 r r 31 -	) / 2	)	( asin (	21 r r 12 -	) / 2	)

Auteurs : AK, AP

Table 6

 6

	.2	Identification errors for the translation [mm]
	Method	0.01 mm	Translation amplitude a 0.1 mm 1.0 mm	10 mm
	SVD	10 -16	10 -16	10 -16	2•10 -14
	LIN	10 -16	10 -16	10 -16	3•10 -15

Table 6

 6

	.3	Parameters of the FEA-modeling noise for different mesh type
	Type of the mesh	Meshing options Abbreviation	Size of the finite element	σ, mm
			3L	3 mm	4.10•10 -5

Table 7 .

 7

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

Table 7 .

 7

	Projet COROUSSO Livrable n°1.2	2
	Outils pour l'identification des paramètres
	de raideur des robots à l'aide d'un logiciel	24/02/2012
	de CAO	indice A
		Page 24/32

Table 7 .

 7

	3	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

Table 8 .

 8 Outils pour l'identification des paramètres de raideur des robots à l'aide d'un logiciel de CAO 2 Comparison of the compliance matrix elements obtained from different methods

	Projet COROUSSO Livrable n°1.2	2
	24/02/2012
	indice A
	Page 29/32