Alexandr Klimchik

Alexandre Ambiehl

Sébastien Garnier

Benoît Furet

Anatol Pashkevich

Benoit Furet

Efficiency evaluation of robots in machining applications using industrial performance measure

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

High-speed machining is quite a new application of industrial robots since previously they were mainly used in automotive manufacturing, for part handling and welding [START_REF] Kochan | BMW chooses innovative robotic machining system for axle carriers[END_REF][START_REF] Terrier | Qualification of parallel kinematics machines in high-speed milling on free form surfaces[END_REF] . As follows form related study [START_REF] Chen | Robot machining: recent development and future research issues[END_REF], the machining segment represents less than 5% of the total market of the industrial robots, but the share is continuously increasing. So, replacement of conventional CNC machines by more competitive industrial robots becomes more and more attractive. The main restraint here is rather limited knowledge of robotics by potential customers and lack of competence of the of robotic cell end-users. Besides, replacing CNC machines by robots leads to additional management expenses. On the other side, the research labs have already confirmed that CNC machines replacement by robots gives essential benefits (reduced product cost, increasing of the work-cell flexibility), which must be clarified for practicing engineers. For this reason, this paper proposes an industry oriented technique for evaluation of the robot capacities in machining, which can be used as the base for the related comparison study.

Compared to conventional CNC machines, robots are able to process complex bulky 3D shapes and provide large and easily extendable workspace that can be modified by adding extra axes. Besides, the same workspace can be shared by several robots. However, the robot trajectory generation is much more complex task compared to the Cartesian machines since mapping from the actuator space to the operational space is highly non-linear. Nevertheless, the results obtained for the tool path optimization of CNC machines [START_REF] Ilyukhin | Nonlinear Adaptive Correction of Continuous Path Speed of the Tool for High Efficiency Robotic Machining[END_REF][START_REF] Makhanov | On the tool-path optimization of a milling robot[END_REF][START_REF] Olabi | Enhanced trajectory planning for machining with industrial six-axis robots[END_REF][START_REF] Lartigue | CNC tool path in terms of B-spline curves[END_REF] can be also applied to robotic cells. Another difficulty may arise because of robot redundancy with respect to the technological process. In fact, conventional machining process requires 5 dof only while most of industrial robots have 6 actuators. This redundancy can be used to optimize the tool path, to improve the trajectory smoothness [START_REF] Zargarbashi | The Jacobian condition number as a dexterity index in 6R machining robots[END_REF] or to reduce the joint torque in order to minimize the impact of machining forces on the robot behavior [START_REF] Zargarbashi | Posture optimization in robot-assisted machining operations[END_REF][START_REF] Vosniakos | Improving feasibility of robotic milling through robot placement optimisation[END_REF].

Another difficulty of robot application in machining is related to non-negligible compliance of robotic manipulators available on the market. For instance, in some cases the end-effector deflections due to the influence of the cutting forces may overcome 10 mm [START_REF] S.-I. Matsuoka | High-speed end milling of an articulated robot and its characteristics[END_REF].T o reduce them, robot manufacturers pay particular attention to improvement of manipulator stiffness and compensation of the compliance errors using dedicated mechanism and/or special control algorithms. To improve the manipulator stiffness, designers are obliged either to increase the link cross-section or to use advanced composites materials. It is clear that the first solution leads to increasing of manipulator moving masses and consequent reduction of dynamic properties. In contrast, utilization of composite materials essentially influences on the ⁎ Corresponding author.

E-mail address: a.klimchik@innopolis.ru (A. Klimchik).

MARK

(5 µm) [START_REF] Ramesh | Error compensation in machine toolsa review: Part I: geometric, cutting-force induced and fixture-dependent errors[END_REF]. Traditionally, they are used for processing of metal parts from parallelepiped-like crude products with high material removal rate. These operations are quite popular in automotive industry, aeronautics as well as in mold making, prototyping, etc. Contemporary CNC machines possess quite large workspace allowing essentially increase an application area. Besides, their efficiency was also proved for processing of composite materials that are utilized more and more due to perfect mass-to-strength ratio. In addition to milling, finishing and trimming operations can be also performed by CNC machines at present. Nevertheless, in spite of numerous advantages, the CNC machines remain very expensive and their workspace is limited and cannot be extended, which is crucial for aeronautic and shipbuilding. This motivates users to find an alternative solution.

One of the promising ways to overcome the above mentioned difficulties is replacing the CNC machines by industrial robots, whose cost is competitive and workspace can be easily extended (by adding extra actuated axes). An example for such an application is presented in Fig. 1. It worth mentioning that machining is fairly new application for robots. Traditionally, the market of industrial robots is shared between handling, pick and place, assembling and welding. The processing (including machining), represents insignificant part of the market, less than 5%. According to PWC study [START_REF] Mccutcheon | The new hire: How a new generation of robots is transforming manufacturing[END_REF] this shares will remain the same in nearest future, as still there are a lot of pick and place operations to be automated with the improvement in robotic vision. Nevertheless, the share of robot-based machining is continuously growing and about a quarter of new robotic cells in North America are processes oriented ones. Large part of this market share corresponds to trimming that was traditionally a high-qualified manual work, but nowadays the robots become competitive here due to increasing of their accuracy. For machining, robots are attractive due to their large and extendable workspace and competitive price that makes them a cost-effective solution for machining of large dimension parts. However, the main obstacle for robots utilization in machining is their relatively low accuracy (about 0.7 mm) and repeatability (about 0.2 mm) compared to the CNC machines, assessment of robot capability for this application have been explored by Slamani et al. [START_REF] Slamani | Assessment of the positioning performance of an industrial robot[END_REF][START_REF] Slamani | A comparative evaluation of three industrial robots using three reference measuring techniques[END_REF] with several methods. It worth mentioning that under the cutting force influence, the positioning errors can reach up to several millimeters. Nevertheless, there are a number of efficient solutions to reduce manipulator positioning errors that were discovered in research labs and progressively applied in industrial environment. The latter allows robots to compete with CNC machines in terms of accuracy, while providing essentially larger workspace. In more details, comparison of CNC machines and robots for machining is presented in Table 1.

Cutting forces in machining

The problem of the cutting forces evaluation is in the focus of mechanical engineers more than 70 years. In his pioneer work, Merchant used principle of minimum angle to develop an analytical robot price and decreases its market competitiveness. Nevertheless, both ways improve the link stiffness only, while the major manipulator elasticity is often concentrated in the actuator gears [START_REF] Dumas | Joint stiffness identification of industrial serial robots[END_REF] and can be hardly improved in practice. Another method of the compliance error reduction is based on the mechanical gravity compensators (pneumatic, hydraulic, spring-based, etc.), which produce auxiliary forces/ torques reducing impact of the link weights. However, this solution does not allow compensating the impact of the machining forces but only shifts the force-deflection relation. To overcome the problem of elastic deformations in the actuator gears, robot manufactures tends to use secondary encoders attached to the motor shaft [START_REF] Devlieg | Expanding the use of robotics in airframe assembly via accurate robot technology[END_REF] that allows to modify the actuator input in order to compensate the gear compliance. It is obvious that this approach also increases the robot price. According to our experience, the double encoders enable compensating about 65% of the compliance errors on average, but in some workspace areas the compensation level is limited by 40-50%. The main reason for this is that the robot link deformations are outside of the double encoder observability. It is clear that for the high-speed milling, where the cutting forces are high enough to cause deflection of several millimeters, such level of error compensation is not enough sufficient. In this case, it is reasonable to apply the off-line error compensation technique [START_REF] Klimchik | Compliance error compensation technique for parallel robots composed of non-perfect serial chains[END_REF][START_REF] Chen | Spatial statistical analysis and compensation of machining errors for complex surfaces[END_REF][START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF] based on the modifying the target trajectory used as the controller input. As follows for our previous research, this approach is very efficient. In particular, the off-line technique based on the simple (reduced) manipulator stiffness model allows user to compensate 85-90% of the end-effector deflections [START_REF] Dumas | Joint stiffness identification of industrial serial robots[END_REF][START_REF] Klimchik | Compliance error compensation in robotic-based milling[END_REF][START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF], while the complete stiffness model ensures the compensation level of about 95% [START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF]. However, it should be stressed that usually robot manufacturers do not provide customers with the manipulator stiffness parameters, so they must carry out dedicated experimental study to obtain the desired model [START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF][START_REF] Abele | Modeling and Identification of an Industrial Robot for Machining Applications[END_REF][START_REF] Guo | Stiffness-oriented posture optimization in robotic machining applications[END_REF][START_REF] Nubiola | Absolute calibration of an ABB IRB 1600 robot using a laser tracker[END_REF][START_REF] Dumas | Joint stiffness identification of six-revolute industrial serial robots[END_REF]. In this paper, the above mentioned problem will be also considered.

To advance robot application in machining, end-user should be provided with clear and efficient tool allowing to evaluate the final product quality expressed via the level of the end-effector deflections caused by the manipulator elasticity. These deflections can be computed both for a single work point and given force/torque or for a set of given trajectories and corresponding cutting forces [START_REF] Klimchik | Serial vs. quasi-serial manipulators: Comparison analysis of elasto-static behaviors[END_REF]. It is clear that usual approach based on different indices extracted from the Cartesian stiffness matrix [START_REF] Guo | Stiffness-oriented posture optimization in robotic machining applications[END_REF][START_REF] Nagai | A systematic approach to stiffness analysis of parallel mechanisms[END_REF][START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF] are not suitable here. In particular, commonly used performance measures based on the stiffness matrix singular values do not represent directly the machining accuracy which is the primary indicator for practicing engineers. For this reason, this paper proposes an industry oriented technique allowing to examine particular robot suitability for a give machining task and to compare several robot-based implementations.

To address the above mentioned problems, the reminder of the paper is organized as follows. Section 2 deals with the particularities of the robot based machining. Section 3 presents stiffness modeling background for industrial robots. Section 4 introduces industry-oriented performance measure adapted for machining accuracy evaluation. Section 5 contains comparison analysis of several industrial robots available on the market that are suitable for high-speed machining. Section 6 deals with experimental validation of the main theoretical results. In Section 7, the limitations of the proposed approach and perspectives are discussed. Finally, Section 8 summarizes the main contributions of the paper.

Robot-based machining

Particularities of machining with robots

Machining with robots is an intersection of two engineering fields: conventional machining and robotics. Machining sector usually prefers for these operations Cartesian CNC machines, which provide end-users with high reliability, good repeatability (2 µm) and very good precision 1D cutting force model as a function of the chip thickness and material properties [START_REF] Merchant | Mechanics of the metal cutting process. I. Orthogonal cutting and a type 2 chip[END_REF]. Later, this model has been extended to the 2D and 3D cases [START_REF] Doi | Chatter vibration of lathe tools[END_REF] and applied to orthogonal cutting such as turning [START_REF] Salje | Self-excited vibrations of systems with two degrees of freedom[END_REF] and milling [START_REF] Chérif | Generic modelling of milling forces for CAD/CAM applications[END_REF]. These results allowed preventing the chatter phenomenon leading to machine and tool damage. Further advances in this area lead to the mechanistic models of the cutting force proposed by Marttelloti [START_REF] Martellotti | An analysis of the milling process[END_REF] and Koenigsberger [START_REF] Koenigsberger | An investigation into the cutting force pulsations during milling operations[END_REF]. These models describe the machining process behavior as the function of the cutting tool geometry and some essential process parameters (such as feed rate, spindle speed, etc). Later, they were integrated in commercial software [START_REF] Weck | CAD assisted chatter-free NC tool path generation in milling[END_REF][START_REF] Shirase | Cutting force and dimensional surface error generation in peripheral milling with variable pitch helical end mills[END_REF]. The mechanistic models were able to predict deformation of the tool during machining and to adapt cutting conditions in order to guaranty high quality of the final product.

Nowadays the cutting force models allow user to take into account more details concerning the application area and cutting conditions (tool geometry, for instance). For example, in the case of the plain milling tool, the flute angle is the constant defining by the tool geometry [START_REF] Wan | A unified instantaneous cutting force model for flat end mills with variable geometries[END_REF], while for the carbide insert tool, the flute angle is a variable that varies during technological process [START_REF] Araujo | Thread milling as a manufacturing process for API threaded connection: geometrical and cutting force analysis[END_REF]. Contemporary fully mechanistic models include more physical parameters (both tool geometry and cutting conditions) [START_REF] Ali | FEM to predict the effect of[END_REF]42], and allow user to estimate cutting forces for different applications. They rely on a special cutting coefficient that characterizes the specific pressure for the cutting tool and workpiece material.

For the face milling application considered in this paper, the carbide insert tool is used and corresponding cutting force model is able to analyze the path of each teeth and to predict corresponding force [42,[START_REF] Gygax | Dynamics of single-tooth milling[END_REF] taking into account the cutting conditions as well as cutting tool geometry and orientation. The geometry of technological process and the cutting force components are presented in Fig. 2. In the tool coordinate system the instantaneous force

FF F =[ ,0, ] i n i f
i T is defined by two principal components: the normal force F n acting perpendicular to the tool face and the friction force F f acting along the face. Here, the superscript "i" indicates the nature of the force, i.e. instantaneous. In the above decomposition, the friction force F f i can be expressed via the normal force F n i using the friction coefficient k f , i.e.

p Fφ Kfa φ () = ⋅⋅ ⋅ s i n n i n z p (1)
where

K Nm m ,[/ ]

n 2 is a specific cutting coefficient, which depends on the material properties and the cutting tool. Here, φ sin takes into account variation of the cross section of the chip depending on the tool advance in the material.

For computational convenience, the instantaneous force F i is presented in the cylindrical coordinate system as

FFF F =[ , , ] i r i t i z iT 0
, where FFF ,, rtz are the radial, tangential and axial (in z-direction) components, respectively. The correspondence between these forces is defined by the rotation matrix κλ γ R (, , ) s 0 that describes the tool orientation, i.e. entering angle κ, helix angle λ s and cutting angle γ 0

κλ γ FR F =( ,,) ⋅ i s i 0 0 (2)
Using the above presentation, it is possible to project the cutting force in the workpiece frame (in the Cartesian space) as

φ FRF =( ) ⋅ C i z i 0 (3)
where the subscript "C" indicates the Cartesian space, and R z is the homogeneous rotation matrix around z-axis by the engagement angle φ. Since the cutting force varies with the engagement angle φ,i t reasonable to average it for the turn

∫ F π Fd φ = 1 2 n φ φ n i S E ( 4
)
where the angles φ S and φ E define the tool engagement interval. More details concerning these angles are given in Fig. 3. It is obvious that they depend on the technological task and cutting tool diameter. Similarly, the average cutting force in the Cartesian coordinate system can be computed as

∫ π φκ λ γd φ FR R F = 1 2 () ⋅ (, , ) ⋅ C φ φ zs i 0 0 S E (5)
Taking into account that

F F C =⋅ n with k C = [1, 0, ]
f T , one can rewrite the above expression in the final form

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ∫ π Fφ φ d φ κλ γ FR R C = 1 2 () ⋅ () ⋅ (, , ) ⋅ C φ φ nz s 0 S E ( 6
)
allowing user to present the cutting force is a function φφ κλγ Kfa F (, , , ,, ,, )

C SE sn 0 z 0 depending on a number of physical and geometric parameters (cutting conditions, tool orientation and tool geometry, material properties). The latter will be used in the following sections to compute the external force applied to the robot end-effector. In a more general case if there are several engagement zones (see Fig. 3c,f for example) it is required to integrate the cutting force for each of them φφ (, )

S j E j , which yields ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ∫ ∑ π Fφ φ d φ κλ γ FR R C = 1 2 () ⋅ () ⋅ (, , ) ⋅ C j n φ φ nz s =1 0 S j E j ( 7
)
where n is the number of engagements for the turn. It should be mentioned that the specific cutting coefficient K n depends on the material and tool properties, so it can be hardly estimated analytically.

For the conventional tools and materials its values can be found in relevant manuals. Also the specific cutting coefficient can be estimated using the orthogonal-to-oblique method based on the material database [START_REF] Wan | Mechanics and Dynamics of Multifunctional Tools[END_REF][START_REF] Kaymakci | Unified cutting force model for turning, boring, drilling and milling operations[END_REF]. However, for the particular applications with specific tool-material couple, this coefficient can be estimated from the dedicated experimental study only [START_REF] Jayaram | Estimation of the specific cutting pressures for mechanistic cutting force models[END_REF].

Once the cutting force is computed, it is possible to estimate the torque applied to the robot end-effector while machining. It can be obtained as the cross product of the tool radius vector rcosφ rsinφ r =[ , , 0 ] T and the cutting force: M rF =× 0 . Here, r is the tool radius. In more convenient form, the torque expression can be presented as

Fκ λ γ M rR C =⋅ ( × ) ⋅( ,,) ⋅ ns 0 (8)
where r (× ) is the skew-symmetric matrix based on the vector r. The above presented models have been implemented as a Matlab toolbox that integrates analytical expressions for different stages of machining process. Example of its application is presented in Fig. 4.It has been obtained for the cutting tool with three teeth (z = 3) of the radius Rm m =5 and the following machining parameters:

κ =9 0 °, γ =7 °0 , λ =4 5 °s , f mm tooth =0 . 0 8 / z , am m =5 p , K Nm m = 750 / n 2 .
Here, the specific cutting coefficient K n was estimated experimentally for the particular tool and workpiece material. As follows from the figure, the maximum cutting force/torque correspond to the complete engagement (i.e. φ = 180 °) and about 1.5 kN and 10Nm respectively. As follows from our experience, these values are non-negligible for typical industrial robots and may cause essential deflections of the end-effector (about several millimeters), which obviously affects the quality of the final product. In the following sections, the above presented force/ torque will be used as reference (i.e., a benchmark task) to estimate the manipulator compliance errors and to compare industrial robot performances for machining application.

Manipulator compliance errors and their impact on machining precision

The cutting forces considered in the previous section cause nonnegligible deflections of the robot end-effector that influence on the precision of the machining process. This section provides mathematical background for the estimation of the deflection magnitude taking into account manipulator component elasticity.

Stiffness model of an industrial manipulator

In robotics, manipulator stiffness is usually described by the Cartesian stiffness matrix K C that allows user to compute the end- effector deflections Δt as ΔtKF =⋅ C -1 [START_REF] Zargarbashi | Posture optimization in robot-assisted machining operations[END_REF] for given external force/torque F, which in the case of machining corresponds to the cutting force. This matrix is positive definite, it essentially varies throughout the robot workspace and highly depends on the robot configuration. There exist different modeling approaches that allow obtaining desired stiffness matrix: the Finite Elements Analysis (FEA), the Matrix Structural Analysis (MSA), and the Virtual Joint Modeling method (VJM). Their detailed comparison, including advantages and limitations, can be founded in [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF][START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF][START_REF] Piras | Dynamic finite-element analysis of a planar high-speed, high-precision parallel manipulator with flexible links[END_REF][START_REF] Chen | Conservative congruence transformation for joint and Cartesian stiffness matrices of robotic hands and fingers[END_REF][START_REF] Marie | Elasto-geometrical modeling and calibration of robot manipulators: application to machining and forming applications[END_REF]. From our experience, the most attractive for our application is the VJM method, which extends the conventional rigid model of manipulator by introducing several virtual joints describing elasticity of principal components [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and nonperfect parallel manipulators under internal and external loadings[END_REF]. Using the VJM-based technique, the stiffness model of a typical industrial robot is presented as a serial chain containing a number of rigid links separated by actuators and virtual springs [START_REF] Dumas | Joint stiffness identification of six-revolute industrial serial robots[END_REF][START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF] as shown in Fig. 5.

To compute the desired Cartesian stiffness matrix, it is necessary to consider simultaneously the extended geometric model tg q θ =( ,)and static equilibrium equation JF Kθ ⋅= ⋅ θ T θ . The first of them allows computing the end-effector location t for given coordinates of the actuators q and virtual joint deflections θ caused by the external loading F. The static equilibrium equation is derived using the virtual work principal and allows to find relation between the external force/ torque and the deflections of the virtual joints. It includes the matrix K θ describing the virtual joint elasticities (joint stiffness matrix) and the Jacobian of the extended geometric model with respect to the virtual joint coordinates Jg q θ θ =∂ ( , )/∂ θ . Simultaneous solution of the above mentioned equations leads to the following expression for the desired Cartesian matrix [START_REF] Yi | Geometric analysis of antagonistic stiffness in redundantly actuated parallel mechanisms[END_REF][START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF][START_REF] Klimchik | Stiffness matrix of manipulators with passive joints: computational aspects, Robot[END_REF]

KJ K J =( • • ) Cθ θ -1 θ T -1 (10)
which obviously depends on the manipulator geometry (via the Jacobian J θ ) and its elastostatic properties (via joint stiffness matrix K θ ). In more general case, when the external loading F is rather high and changes essentially the manipulator equilibrium configuration, the Cartesian stiffness matrix also depends on the Hessian H θθ (second order derivative of the extended geometric model), taking the form [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and nonperfect parallel manipulators under internal and external loadings[END_REF][START_REF] Klimchik | Static stability of manipulator configuration: influence of the external loading[END_REF]]

KJ K H J =( • ( - ) • ) Cθ θθ θ -1 θ T -1 (11)
where

H gq θ F θ =∂ ( ( , ) ⋅ )/∂ θθ 2
2 describes modification of manipulator elasticity due to applied loading.

It should be stressed that the above expressions for the Cartesian stiffness matrix (10)- [START_REF] S.-I. Matsuoka | High-speed end milling of an articulated robot and its characteristics[END_REF] were derived assuming that joint stiffness matrix K θ is constant. However, if the manipulator includes the gravity compensators (Fig. 6a,b), the equivalent joint stiffness coefficients become configuration dependent [START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF] and relevant joint stiffness matrix is presented as a sum

K KK q =+ ( ) GC θ θ 0 θ (12)
where the first term K θ 0 is constant and corresponds to the manipulator without compensator (see K θ in Eq. ( 10)); and the second term K q ( ) GC θ depends on the manipulator configuration and describes equivalent elasticity of the gravity compensator with respect to the virtual joints.

As follows from relevant study, the matrix K q ( ) GC θ is rather sparse and contains only several non-zero elements which are described by highly non-linear functions. For example, for the spring-based gravity compensator in the second actuated joint (Fig. 6a,c), the matrix K q ( )

GC θ is expressed as follows diag K a L η K q () = ( 0 , • • ,0 ,0 ,0 ,0 ) GC GC q θ (13)
where K GC is the compensator spring stiffness, a and L are the geometric parameters (see Fig. 6c), and η q is the configuration- dependent coefficient expressed as

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ η αq s s aL s αq αq = c o s (+) - • s i n (+) + c o s (+) q 0 2 2 (14)
Here, s 0 is the compensator spring preloading, q is the actuated coordinate, α and s are auxiliary geometric parameters defining compensator configuration whose meaning is clear from Fig. 6c.

Hence, knowing relevant geometric and elastic parameters of the manipulator, it is possible to compute the Cartesian stiffness matrix K C and the undesired end-effector deflections caused by the machining force F. However, in practice, robot specifications usually include the geometric parameters only, while the elastic parameters should be identified by the users from the dedicated experimental study. This problem is in the focus of the following section.

Identification of the stiffness model parameters

Methodology and basic equations

In practice, the stiffness model parameters of industrial robots are not included in datasheets and identified from the appropriate experimental data [START_REF] Joubair | Non-kinematic calibration of a six-axis serial robot using planar constraints[END_REF]. Relevant experimental setup usually requires openloop or/and closed-loop measurement system [START_REF] Nubiola | Absolute calibration of an ABB IRB 1600 robot using a laser tracker[END_REF][START_REF] Takeda | DBB-based, kinematic calibration method for in-parallel actuated mechanisms using a Fourier series[END_REF][START_REF] Ikits | Kinematic calibration using a plane constraint[END_REF] as it is presented in Fig. 7. As a rule, the identification procedure contains the following steps [START_REF] Dumas | Joint stiffness identification of industrial serial robots[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]:

(i) Selection of optimal measurement configurations ensuring the lowest impact of the measurement errors. At this step, a customized design of experiments technique is applied that is based on the appropriate performance measure (RMS of the positioning error after the compliance errors compensation [START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF][START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF][START_REF] Imoto | Optimal kinematic calibration of robots based on maximum positioning-error estimation (Theory and application to a parallel-mechanism pipe bender[END_REF], indices computed using singular values of the manipulator Jacobian [START_REF] Klimchik | Optimal pose selection for calibration of planar anthropomorphic manipulators[END_REF][START_REF] Joubair | Comparison of the efficiency of five observability indices for robot calibration[END_REF][START_REF] Sun | Observability index selection for robot calibration[END_REF][START_REF] Borm | Determination of optimal measurement configurations for robot calibration based on observability measure[END_REF], etc.). The importance of this step was shown in [START_REF] Klimchik | Optimal pose selection for calibration of planar anthropomorphic manipulators[END_REF][START_REF] Borm | Determination of optimal measurement configurations for robot calibration based on observability measure[END_REF]. (ii) Measuring reference point locations without loading. The output of this step is a set of Cartesian coordinates of the reference point and corresponding values of the joint angles obtained from the actuator encoders. (iii) Application of the external loading to the manipulator endeffector (specially designed for this experiments [START_REF] Dumas | Joint stiffness identification of six-revolute industrial serial robots[END_REF][START_REF] Klimchik | Optimal selection of measurement configurations for stiffness model calibration of anthropomorphic manipulators[END_REF]). The magnitude and direction of this loading should be either known or measured in experiments. To increase the identification accuracy, it is reasonable to apply maximum allowed force/torque but do not violate the stiffness model linearity. (iv) Measuring reference point locations with loading. Here, the output is similar to the step (ii). To avoid influence of other factors (such as reputability, heating, etc.), the measurements for unloaded and loaded manipulator should be carried on sequentially, i.e. without changing configuration between steps (ii)-(iv). (v) Repeating steps (ii)-(iv) for all measurement configurations obtained in step (i).

(vi) Estimation of the stiffness model parameters using dedicated numerical procedure and set of the experimental data (end-effecto locations with and without loading, corresponding joint coordinates and loading magnitude/direction).

Basic expression for the stiffness model parameters identification is derived from the force-deflection relation [START_REF] Zargarbashi | Posture optimization in robot-assisted machining operations[END_REF]. It is usually written in the following way [START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF][START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF]

Δ jn tA π =•; = 1 , jj (15)
where Δt j is the end-effector deflection caused by the loading, the vector π aggregates all desired parameters (compliances) to be identified, the subscript j defines the measurement experiment number, A j is so-called observation matrix corresponding to the jth experiment that is presented as

AJ J F J J F J J F =[ , ,. . . ., ] jj j T jj j T jj m jm T j 1 1 2 2 (16)
The latter is expressed via the external force F j and the columns of the Jacobians JJ J J = [ , , ...., ]

jj j j m θ1 2
of the extended geometric model for the manipulator configuration q j .

After application of the least square technique, the identification procedure reduces to the following optimization problem

∑ ΔΔ Aπ t η η Aπ t (- ) (- ) → m i n j m jj T j T j jj π =1 (17)
which should be solved with respect to the vector of the stiffness parameters π. Here η j is the matrix of the weighting coefficients that normalizes the measurement data (normalization is obligatory if position and orientation components are used simultaneously because of different units). It should be noted that the weighting matrices η j can be similar for all experiments if the noise distribution remains the same for all measurements. The above optimization problem is solved in a usual way and yields the following expression for the set of the unknown parameters

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ∑∑ Δ π Aη ηA Aη η t ˆ=• . j m j T j T j j j m j T j T j j =1 -1 =1 (18)
If the measurement noise is Gaussian (as it is observed in practice), expression [START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF] provides us with an unbiased estimate. Corresponding covariance matrix evaluating the dispersion of π ˆcan be computed as follows

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ∑∑ ∑ π A ηηA A ηηΣηηA A ηηA cov( ˆ)= j m πj T j T j πj j m πj T j T j j j T j πj j m πj T j T j πj =1 -1 =1 2 =1 -1 (19)
where the matrix Σ j 2 describes the variance of the measurement errors in the j th experiment.

Case of manipulators with gravity compensators

If robot includes the gravity compensator, the force-deflection relation becomes essentially non-linear and can be hardly presented in the form [START_REF] Chen | Spatial statistical analysis and compensation of machining errors for complex surfaces[END_REF]. Nevertheless, using specificworkaroun d [START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF],thesetoftheunknown 6 . Here, the subset ππ ( , ...) 21 22 represents the second joint compliances corresponding to different values of the joint variable q 2 . This idea allows us to present the basic equations in similar way [START_REF] Abele | Modeling and Identification of an Industrial Robot for Machining Applications[END_REF] and to keep the same expressions ( 18)-( 19) for the parameter estimates while replacing A j by B j .T h em o d i fied observation matrix B j is composed of the matrix A j elements, it is obtained by inserting in the matrix A j several zero columns that correspond to the parameters from the subset ππ ( , ...) 21 22 different from the current one.

Δ jn tB π =⋅ ; = 1 , j j ext
On the second step, using the sub-set ππ (, . . . )

ii 12

, it is possible to separate the stiffness parameters of the manipulator and the compensator parameters. For example, for the spring-based compensator, one can find separately the actuated joint stiffness K θ 0 as well as gravity compensator stiffness coefficient K GC and its preloading s 0 .

Corresponding equations can be also solved using the least square technique that yields an expression

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ∑∑ KK s K K CC C [ • ]= θ GC GC T i m i T i i m i T θi 0 0 =1 -1 =1 GC GC ( 21
)
where m GC is the number of unique angles q i in the set of measurement configurations, and the matrix C j is defined as follows

a L αq a L s a L s αq αq C =[ 1 - • cos ( -) / •( / •sin ( -) + cos( -))] i i i ii 22 (22)
Thus, this modification of the conventional identification technique developed for strictly serial manipulator allows us to find the stiffness model parameters for a more general case, for the manipulators with gravity compensators.

Identification using enhanced partial pose measurements

Conventional full pose measurement technique [START_REF] Hollerbach | Model Identification[END_REF][START_REF] Driels | Full-pose calibration of a robot manipulator using a coordinatemeasuring machine[END_REF] provides both the position p and orientation φ that above were merged in a single location vector tp φ =( , ). However, in practice conventional measurement systems (laser tracker Leica or Faro, for instance) do not provide with the orientation components φ directly, so these angles are computed using the Cartesian coordinates of several reference points [START_REF] Driels | Automated partial pose measurement system for manipulator calibration experiments[END_REF][START_REF] Goswami | Complete parameter identification of a robot from partial pose information[END_REF][START_REF] Santolaria | Laser tracker-based kinematic parameter calibration of industrial robots by improved CPA method and active retroreflector[END_REF] as shown in Fig. 8. Further, to include both p and φ in the identification equations, the full-pose measurement technique requires normalization (see Eq. ( 17)) in order to avoid the non-homogeneity problem. However, this normalization may affect the identification results and it is reasonable to exclude it. For this reason, to avoid the above mentioned difficulties, it was developed the enhanced partial pose measurement technique [START_REF] Klimchik | Geometric and elastostatic calibration of robotic manipulator using partial pose measurements[END_REF], which uses the Cartesian coordinate measurements only.

Using this idea, the identification equations can be written in a reduced form, which yields the following expression for the desired set of parameters π ext

⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ∑∑ ∑∑ Δ π BB B p ˆ= ext j m i n j ip j ip j m i n j ip j i =1 =1 () () -1 =1 =1 () TT ( 23
)
where the superscript "(p)" indicates the position components (i.e. first three lines in the case of matrix B p ( ) ), the indices "j" and "i"d e fine the manipulator configuration number and the reference point number, respectively, m and n are the number of configurations and reference points used. Essential benefit of this approach comparing to the full pose measurement technique is potential improvement of the identification accuracy due to increasing the total number of the scalar equations from m 6 to mn 3 . The principal particularity of this approach that it operates with several tool transformations p tool i , which should be also identified. It can be proved that these transformations can be computed using the same measurement data as

⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ∑∑ Δ pp D D D p [ ;. . . ]= tool tool n j m jj i j m j j G 1 =1 -1 =1 TT ( 24
)
where the matrix diag D RR =( , . . . , )

j jj n 1
is composed of n rotation matrices R j , the matrix R j defines the robot orientation that is computed from the direct kinematics,

Δ ΔΔ pp p =[ ;. . .; ] j G j G j Gn 1
is the vector of residuals between the measurements (for the unloaded manipulator) and computed positions for zero tool parameters. In the experimental section of this paper, the above presented equations will be used to identify the desired stiffness parameters of industrial robots under study.

It should be mentioned that in the case when a sophisticated manipulator CAD model is available, it is also posible to estimate stiffness coefficient using virtual experiments [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF][START_REF] Pashkevich | Accuracy improvement for stiffness modeling of parallel manipulators[END_REF]. Moreover, this approach allows user to obtaine realistic 6×6 stiffness matrices for each manipulator component separately.

Estimation of compliance errors and their compensation

The above presented stiffness model allows user to evaluate the impact of the machining force on the robot precision. On the other hand, this model can also be used for adjusting the controller input in order to compensate (or to reduce) the manipulator compliance errors. It should be mention that in the controllers of industrial robots, the manipulator motions are usually generated by means of purely kinematic control (inverse and direct kinematics) allowing obtaining actuator coordinates corresponding to the desired end-effector location t 0 without taking into account the external loading. However, in machining application this approach is insufficient. To overcome this difficulty, it is reasonable to integrate in the control scheme the manipulator stiffness model that in general case can be written as

f tF t =( | ) F -1 0 (25)
where the subscripts 'F' and '0' refer to the loaded and unloaded endeffector locations, respectively, and '|' separates arguments of the function f (). More details concerning this function are given in our previous works [START_REF] Klimchik | Compliance error compensation in robotic-based milling[END_REF][START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF]. Using this model, the problem of compliance error compensation can be presented in the following way. Let us assume that the external loading F causes undesired end-effector displacement from the location t 0 to t F . To compensate this error, it is necessary to find another reference location t*, that will be used as the controller input and for which the same loading F causes displacement to the desired location t 0 . This idea leads to the following non-linear equation

f tF t =( | *) 0 -1 (26)
that should be solved with respect to t*. As follows from our experience, the desired solution can be obtained using the simplest iterative scheme

αf tt t F t * = * +• ( - ( | *)) ′ 0 -1 ( 27
)
where α is scalar parameter ensuring iterative scheme convergence, the prime indicates modification of t* on the next iteration. In more details parameters π =[ ,ππ 12 , π 3 ,...,π] 6 can be extended up to π ext , which for a typical industrial manipulator with compensator attached to the second the iterative procedure is presented in Fig. 9. In practice, it can be implemented either on-line or off-line, as shown in Fig. 10.

Accuracy of robot-based machining and its evaluation

Industrial standards for machining accuracy

In industrial practice, there exist a number of norms to evaluate the quality of a final product. They estimate the path straightness (ISO 12780 [START_REF] Standardization | Geometrical product specifications (GPS) -Straightness -Part 1: Vocabulary and parameters of straightness[END_REF]), the surface flatness (ISO 12781 [START_REF] Standardization | Geometrical product specifications (GPS) -Flatness -Part 1: Vocabulary and parameters of straightness[END_REF]) and the path roundness (ISO 12181 [START_REF] Standardization | Geometrical product specifications (GPS) -Roundness -Part 1: Vocabulary and parameters of straightness[END_REF]), which is also called the circularity. Some details concerning their application can be found in [START_REF] Borisov | Performance evaluation of a new taut wire system for straightness measurement of machine tools[END_REF][START_REF] Raghunandan | Selection of sampling points for accurate evaluation of flatness error using coordinate measuring machine[END_REF][START_REF] Wen | Flatness error evaluation and verification based on new generation geometrical product specification (GPS)[END_REF][START_REF] Feng | Straightness error compensation for large CNC gantry type milling centers based on B-spline curves modeling[END_REF]. Robotic trimming of a straight path was studied by Slamani et al. [START_REF] Slamani | Analysis of trajectory deviation during high speed robotic trimming of carbon-fiber reinforced polymers[END_REF] who investigated the influence of the cutting parameters on the final product quality. The circularity norm was used in [START_REF] Oh | Robot accuracy evaluation using a ball-bar link system[END_REF] to estimate robot accuracy and efficiency of the proposed error compensation technique. Similar approach based on the circular tests for boring of aluminum parts for five-axis machining was applied in [START_REF] Slamani | Improvement to highspeed end mill boring accuracy by a simple compensation strategy[END_REF], where the authors compared boring with a boring bar and boring with a highspeed end mill.

From our experience, the circularity norm is the best suited one for the machining process, since it also evaluates the straightness and flatness in indirect way. In fact, the straightness violation leads to the non-uniform circle stretching/shrinkage and the flatness violation causes the circle twisting. On the other side, perfect milling of circular profiles guarantees perfect machining of straight lines and plain surfaces. For this reason, the circularity norm will be used further to evaluate the capacity of industrial robot to perform the machining tasks.

According to relevant standard (ISO 12181), the circularity calculation includes two steps: building a reference circle and estimating the deviations with respect to the reference circle. The standard defines four methods to obtain the reference circle that are called respectively the Minimum Circumscribed Circle (MCC), the Maximum Inscribed Circle (MIC), the Minimum Zone Circles (MZC) and the Least Squares Circle (LSC). The difference between these methods is illustrated in Fig. 11. In all cases, the circularity is equal to the distance between the inscribed and circumscribed circles. The principal difference is related to the center that is computed using different approaches. For example, for MIC the center point is computed for the maximum inscribed circle and it is also used for the minimum circumscribed one. In the MCC method, the center is computed for the minimum circumscribed circle and the inscribed circle is build using the same center point. In the case of LSC, the and circumscribed circles are found for the center point obtained for the least square circle. In contrast, the MZC method uses a center point for which the distance between the inscribed and circumscribed circles is minimal. It should be stressed that in the case of real measurement data all methods are competitive and provide almost the same results [START_REF] Li | The relationship between the minimum zone circle and the maximum inscribed circle and the minimum circumscribed circle[END_REF].

Usually, the reference circle is defined by a particular application. In particular, for the hole-milling it is reasonable to use the MIC reference circle that estimates the ability of the corresponding cylinder to be inserted into the hole and the circularity characterizes the clearance. In contrast, for machining of the cylinder that should enter into the hole, the MCC is more appropriate. However, MIC and MCC methods are not applicable if machining profile is asymmetrical. In this case it is prudent to use either MZC or LSC methods. In fact, the MZC evaluates true circularity; however it is the most complicated approach from numerical point of view. Thus, a reasonable alternative in engineering practice is the LSC method. Advantages of this method are confirmed by simulation results presented in Fig. 12, which show that difference between MZC and LSC is negligible in the case of conventional cutting conditions [START_REF] Calvo | Accurate evaluation of functional roundness from point coordinates[END_REF].

Evaluation of circularity for robot-based machining

Basic expressions

As it was mentioned above, the circularity is the performance measure that characterizes the difference between the radii of maximum inscribed and minimum circumscribed concentric circles obtained for the reference machining profile ρ rr =max min [START_REF] Mccutcheon | The new hire: How a new generation of robots is transforming manufacturing[END_REF] In the frame of the LSC method, the center point is obtained by solving the following optimization problem

∑ xx yy r ( ( -) + ( -) -) → min i n i i xy r =1 0 2 0 22 ,, 0 0 (29)
where xy (, ) i i are the coordinates of the machining profile provided by the measurement system, xy (, ) 0 0 is the desired center of least-square circle and ris its radius. The optimization problem ( 29) is highly non- linear one and it cannot be solved analytically. For this reason, the Newton-Raphson method is used. To improve the convergence, the start point can be estimated from the simplified optimization problem [START_REF] Kim | Geometrical tolerances: improved linear approximation of least squares evaluation of circularity by minimum variance[END_REF] ∑

xx yy (( -) + ( -) ) → min i n i i xy =1 0 2 0 2 , 0 0 (30) that yields ∑∑ ∑ x n xy n y r n xx yy = 1 ,= 1 , = 1 (- )-(- ) start i n i start i n i start i n i start i start 0 =1 0 =1 =1 0 2 0 2 (31)
After finding the center point xy (, ) 0 0 , the desired radii of the circumscribed and inscribed circles are be computed as follows is the center of corresponding least-square circle. Hence, for given set of points p i describing machining profile one can compute the circularity index ρ using expressions presented above. It is clear that these points can be obtained either experimentally or numerically, using relevant models of the manipulator (whose parameters are identified from the dedicated experiments, see Section 3.2). This paper concentrates on the second approach that requires the manipulator stiffness model only.

Modeling of the machining profile

This approach assumes that all geometric and elastic parameters of the manipulator are either given or have been already identified from the dedicated experiments. In this work, the circularity is evaluated for the benchmark profile corresponding to a circular milling task of the radius r mm = 100

. This benchmark is in a good agreement with typical industrial requirements, but it can be easily adapted to other dimensions.

Assuming that machining produces the cutting force/torque FM (, ), the manipulator end-effector deflections Δt i caused by this loading can be computed using expression

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ΔtJ qK qJ q F M =( ) • ( ) •( ) ⋅ i ii i i i θ θ -1 θ T (33)
where the subscript "i" indicatives the machining profile point number, K q ( ) i θ is the manipulator stiffness matrix corresponding to the configuration q i , Jq () i θ is the manipulator Jacobian for the same configuration. Further, taking into account particularities of the machining process, variation of the force/torque FM (, ) along the path can be taken into account via the rotation transformation of corresponding vectors (while the magnitudes remain the same).

So, for the given set of the polar angles φ { } i , it is possible to compute the set of points p { } i describing the machining profile corresponding to the desired circle:

φφ φ pp R rJq K q J q R FJqK qJ q R M = + () ⋅+ ()• ()• () ⋅() ⋅ +( ) •( ) • ( ) ⋅ ( ) ⋅ ii p ii pT ii p ii φT ii 0 θ () θ -1 θ () θ () θ -1 θ () (34)
Here, the radius-vector rdefines the profile point corresponding to φ = 0, φ R( ) i is the rotation matrix, the superscripts " p ( )" and " φ ( )" denote the position or orientation part of Jacobian matrix, i.e. col JJ J =(, ) () . Further, the obtained set of points p { } i is used to compute the desired circularity index ρ in accordance with expressions [START_REF] Mccutcheon | The new hire: How a new generation of robots is transforming manufacturing[END_REF][START_REF] Slamani | Assessment of the positioning performance of an industrial robot[END_REF][START_REF] Slamani | A comparative evaluation of three industrial robots using three reference measuring techniques[END_REF][START_REF] Merchant | Mechanics of the metal cutting process. I. Orthogonal cutting and a type 2 chip[END_REF][START_REF] Doi | Chatter vibration of lathe tools[END_REF].

pφ θ θ () θ
Using the approach presented in this Section, it is possible to evaluate circularity for the entire robot workspace (i.e. locating center of the circle in different zones) and to determine the zone in which the manipulator provides the best performance from the circularity point of view. In the following Section, the developed approach is applied to compare machining capabilities of several industrial robots used by our industrial partners.

Comparison study: assessment of robots capacities for machining

Examined industrial robots and their parameters

The technique developed in this paper has been applied to the comparison study of five industrial robot of the Kuka family (Fig. 13). They were compared with respect to the circular machining task of 100 mm radius that was placed in different workspace points. The examined robots have similar kinematics (Fig. 14), their geometric parameters are presented in Table 2. In the unloaded mode, the examined robots have comparable repeatability/accuracy in the range from 50 to 80 µm. However, their payload capacities differ essentially and vary from 100 to 500 kg as shown in Table 3. Working radii of the examined robots are vary from 2.6 to 3.9 m.

Elastostatic parameters of the examined robots are presented in Table 4. There were obtained using dedicated experimental study and identification technique developed in our previous work [START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF][START_REF] Klimchik | Geometric and elastostatic calibration of robotic manipulator using partial pose measurements[END_REF].I ti s worth mentioning that for the parameters k 2 describing equivalent compliance of the second joint, it is given the range that takes into account influence of the gravity compensator, which depends on the manipulator configuration (see Section 3 for details).

To estimate the cutting force/torque for the benchmark task it was assumed that the cutting tool radius is 5 mm and number of teeth is equal to 3. The remaining machining parameters are: κ =9 0 °, γ =7 °0 Nm Nm M = [ 0, 3, 1 0 . 5], which will be used in the following section for the circularity evaluation of robot-based machining.

Comparison analysis: robot performance in machining process

To compare accuracy of the examined robots with respect to the machining process, it was considered the circular milling task of radius 100 mm. It was assumed that the benchmark trajectory was located in the plane XOZ, which suits the majority of our industrial problems. For this benchmark task, the circularity maps were computed for all manipulators under study. Relevant results are presented in Figs. [START_REF] Chen | Spatial statistical analysis and compensation of machining errors for complex surfaces[END_REF][START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF][START_REF] Klimchik | Compliance error compensation in robotic-based milling[END_REF][START_REF] Klimchik | Design of calibration experiments for identification of manipulator elastostatic parameters[END_REF][START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF], which include two cases corresponding to different manipulator configurations (elbow-up and elbow-down). It is obvious that manipulator configuration must remain the same during continues machining process in order to avoid passing through the kinematic singularities. In addition, the figures contain optimal regions (red squares) for locating the machining tasks in square zones of different sizes (i.e. 200×200 mm, 400×400 mm and 600×600 mm, etc.). The figures also show the same zones in the middle of workspace (green squares). Summary of these results is presented in Table 5 that, for comparison purposes, also includes the circularity indices computed for the cutting force/torque corresponding to 80% of maximum loading. Some extra indicators for further comparison analysis are presented Tables 6,7, which contain accuracy for square workspace regions in the case of optimal and middle-task placements, respectively.

As follows from the obtained results, the best accuracy for the considered benchmark task is ensured by the robot Kuka KR500. In particular, for this robot the best circularity index is equal to 0.41 mm (for the optimal task placement). This advantage is achieved due to less complaint actuators that obviously increase the robot price. In general, all examined robots except KR120 ensure circularity level about 3 mm within entire workspace. It is clear that so high value of this indicator is not acceptable in practice. But it worth mentioning that by applying relevant compliance error compensation techniques [START_REF] Klimchik | Compliance error compensation technique for parallel robots composed of non-perfect serial chains[END_REF][START_REF] Klimchik | Compliance error compensation in robotic-based milling[END_REF] it is possible to improve the circularity by 80-95%, depending on the stiffness model complexity [START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF]. This allows us to achieve the circularity level of about 0.1 mm.

In practice, selection of a particular robot for the given task should be multi-objective. It is necessary to take into account both the accuracy (represented here via the circularity) and also pay attention to the price, the robot workspace size and its payload capacities. For instance, the robots KR100 and KR120 are not suitable for machining of hard materials because of the payload limitation. However, because of their lower price comparing to other examined robots, they are competitive for machining if the cutting force magnitude is less than 1 kN. For the hard materials, either KR360 or KR500 can be used depending on the required accuracy level and the workpiece dimensions. On the other hand, the robot KR360 is competitive for large- dimensional tasks only and for milling with forces higher than 2.5 kN. Otherwise, KR270 is preferable that may ensure better performance within its workspace. For extremely large-dimension parts, the robot KR120 is preferable due to its high working radius, however its accuracy is relatively low.

To simplify robot selection with regard to different objectives, Fig. 20 presents the normalized bar diagram showing four principal indicators: (i) the circularity indices, (ii) working radius, (iii) maximum payload, and (iv) relative price. Corresponding absolute values of these performance measures are presented in Tables 3,5. The figure clearly shows that increasing of the robot workspace leads to reduction of its accuracy. From the other side, the robot accuracy is well correlated with it price. These results allow practicing engineers to justify robot selection for given machining task taking into account contradicting objectives.

Integrated results of the completed comparison analysis are presented in Tables 8,9. The first of them presents machining accuracy of the examined robots (in terms of the circularity) for different magnitude of the cutting force. It is assumed here that either online or offline compensation algorithm is implemented in order to reduce the manipulator compliance errors. The second table presents the same results in a slightly different way, allowing user to select an appropriate robot from the set of the examined ones knowing the desired accuracy and expected cutting force magnitude. It is worth mentioning that here the accuracy is evaluated for 200×200 mm square zone optimally located within the robot workspace (i.e. so called minimax circularity). Hence, the developed methodology allows user to select an industrial robot taking into account particularities of the given technological task.

Experimental validation

To validate the developed methodology, several machining experiments were carried out using one of the examined robots used by our industrial partners (Kuka KR270) and available in the laboratory. A spindle PRECISE with a power of 24 kW was used to perform the tests. The robotic cell layout is typical for pre-machining tasks for the most of applications we are involved in. The experiments were targeted at milling the circular grooves of radius 60 mm for several workpiece locations. The workpiece locations were chosen according to the available fixture on the machine table, on the same height and different distances from the axis #1. The first placement is the closest one to the robot base, the third placement corresponds to the maximum working radius achieved on the machining table (the second placement is in the middle between the first and the third ones). The main difference with respect to the robot for the selected locations appears in the arm configurations (axis #2 and #3). In our experiments, the grooves were milled by a single step of the provided cutting depth, no intermediate steps were used. The robot trajectory programming was done in the Cartesian space with a fixed tool orientation. It should be mentioned that in spite of the fact that only 5 dof are required to perform the technological task (while the robot has 6 dof), the problem of the redundancy resolution and manipulator configuration optimization along the path is not important here since the milling tool axis coincides with the last axis of the robot wrist. It is worth mentioning also that here, because of some restrictions, the radius was lower comparing to the benchmark task considered in Section 5.

Nevertheless, the obtained data can be reasonably used for the validation assuming that simulation is adapted for these settings. The experimental setup is presented in Fig. 21, it includes the robot Kuka KR270 equipped with a spindle and a machining table with a workpiece attached to. In the experiments the workpiece was placed in several locations (Fig. 22a), the groove depth varied from 4 to 10 mm, and two types of materials were used (soft and hard ones). The machined workpiece are presented in Fig. 22b. To make theoretical and experimental results comparable, the simulation has been repeated for radius 60 mm and the cutting force estimated for this case.

Comparison of the theoretical and experimental results on linear and polar plots are presented in Figs. 23 and24 respectively, which contain the milling path deviations with respect to the reference circle. As follows from them the shapes of the machined grooves and estimated ones are very close to each other. For instance, for the location #1 and the cutting depth 6, 8 and 10 mm, the difference between the measured and estimated circularity indices is less than 5%. This allows us to conclude that the above presented technique for robot comparison is based on realistic performance measures.

Summarized results of the experiment study are given in Table 10, which contain the circularity indices for all machined samples. To evaluate the quality of the studied machining process and to compute the circularity, a standard 3 axis CMM machine was used. It should be noted that the entry and exit points were eliminated from the analysis in order to increase the validity of the circularity estimation. The results clearly show that augmentation of the cutting force magnitude (directly related to the cutting depth) leads to increasing of the circularity index, which can be computed by simple scaling. Also, the obtained results confirm advantages of the workpiece location #1 which is very close to the optimal task placement computed in Section 5. It is worth mentioning that the stiffness-based optimization of the task location is preferable here compared to the convention kinematic-based approaches, which does not take into account the machining accuracy.

Discussions

In spite of numerous advantages, the developed approach has some limitations that should be clarified. It is clear that the adapted performance measure (evaluating contortion of the benchmark circle) is general enough since it indirectly takes into account other desired properties such as the path straightness and surface flatness. From the other side, if the technological task is exactly known (the path and corresponding forces), it is possible to evaluate machining trajectory directly. However, the latter approach may produce the task location optimum that is not acceptable for a slightly different operation.

It is worth mentioning that the proposed methodology requires estimation of the cutting forces that highly depend on the technological tool and material properties, so it is computed rather approximately. Nevertheless, small errors in the force magnitude are not critical here since the circularity-force relation is almost linear. This fact is confirmed by experimental results showing that the groove depth increase from 6 to 8 mm leads to the circularity degradation from 0.69 to 0.89 mm. Hence, the issue of the cutting force magnitude is not essential for the robot comparison because it does not violate the robot dominance order. For this reason, the theoretical model of the cutting force used in this paper is good enough to get reasonable results for robot selection.

Another important subject to be discussed is the validity of the cutting force magnitude/direction assumption. In fact, the developed methodology assumes that the cutting force is constant with respect to the tool feed direction, and the force simply rotates while advancing along the circle. From technological point of view this assumption is equivalent to the constant groove depth condition. However, as follows from the experimental study, the cutting depth essentially varies along the trajectory (about 20%, from 0.8 to 1.3 mm depending on the workpiece location). This aspect was ignored in the circularity computing (see Section 4.2) while it is an important issue to be considered in future. The most obvious solution of this problem is to adjust the cutting force for each path-point taking into account current groove depth estimated via the relevant compliance error component. Such enhancement should obviously make simulations results closer to the experimental data.

It should be also stressed that the comparison analysis does not require the complete stiffness model. In particular, the simplified model taking into account the actuators elasticity only is sufficiently good here, since it allows predicting about 80% of the end-effector elastostatic deflections without violating the robot dominance order. In contrast, for the compliance error compensation the complete stiffness model is preferable, allowing to reach 95% error reduction level.

In the frame of the developed technique, the robots are compared considering the circularity index only. However, the compliance errors (caused by machining) emerge as the circle center shift and modification of its radius also. For instance, in machining experiments it was detected that the radii of the obtained circles differed from the desired one by 0.6-1.0 mm. Similar results are also produced by simulation. Nevertheless, from our point of view, the circularity index is the most representative for the machining accuracy evaluation because the radius modification can be easily eliminated by simple scaling of the input trajectory in the robot control system.

Finally, it is worth mentioning that the developed methodology is based on the performance indices computed using elasto-static model. At the same time, the experimental data characterize both static and dynamic properties of the robot. Nevertheless, there is good agreement between the simulation and experiments, which justifies the methodology proposed in this work.

Conclusions

The paper proposes an industry-oriented methodology allowing user to rank industrial robots with respect to the machining accuracy. The developed methodology is based on the circularity index evaluating precision of the benchmark circular trajectory completed by an examined robot. To compute this performance measure, the manipulator stiffness model is used for estimation of the compliance errors corresponding to the reference cutting force, identical for all considered robots. The cutting force is evaluated taking into account typical machining conditions and material properties. The stiffness model is obtained experimentally using dedicated technique presented in the paper. The main advantage of the accepted performance measure is its clear engineering meaning. Besides, it takes into account the roundness, path straightness and the surface flatness.

Efficiency of the developed technique was confirmed by an application example that deals with comparison study of five industrial robots of Kuka family. For each of the examined robots, the circularity was evaluated for the entire workspace and for all possible manipulator configurations. These results allowed us to get the circularity limits and to find the best task placement for each case study. Final results are presented the forms of the tables and bar diagram, which can be used for selection an appropriate robot taking into account required accuracy, payload and workspace size.

Validity of developed methodology was also proved by experimental study dealing with robot-based machining of circular grooves for several workpiece samples and different locations. The obtained circularity indices were in a good agreement with the theoretical model used in this paper. In addition, advantages of the optimal task placement computed using this performance measure were also confirmed. The experimental results also showed almost linear relation between the circularity and the cutting force magnitude, which allows user to apply the same circularity maps for different machining conditions and to adapt numerical values by simple scaling.

In future, the developed methodology will be enhanced in order to take into account variation of cutting force and cutting depth along the trajectory. Another research direction is related to investigation of machining capacities for parallel robots and their comparison with the serial ones.

Fig. 1 .

 1 Fig. 1. Example of machining process with robot.

Fig. 2 .

 2 Fig. 2. Cutting force components and machining process geometry.

Fig. 3 .

 3 Fig. 3. Tool engagement angles for different machining tasks.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Variation of forces/torques for different stages of machining process Rm m =5 , z =3 , κ = 90°, γ =7 °0 , λ =4 5 °s , fm m r e v =0 . 0 8 / z

Fig. 6 .

 6 Fig. 6. Mechanical compensators of industrial robots (a, b) and their static model (c).

Fig. 7 .

 7 Fig. 7. Experimental setup for identification of stiffness model parameters.

Fig. 8 .

 8 Fig.8. Typical measurement tool for robot calibration with three reference points.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Compliance error compensation algorithm.

Fig. 11 .

 11 Fig. 11. Estimation of the circularity: definition of ISO norms.

Fig. 12 .

 12 Fig. 12. Difference between the circularity norms.

  this case, application the above presented technique (see Section 2.2.) for completely engaged tool (φ = 180 °) yields the following cutting force NN N F = [-440 , -1370 , -635 ]

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Examined industrial robots.

Fig. 15 .

 15 Fig. 15. Circularity maps for robot KR 100 HA, mm.

Fig. 17 .

 17 Fig. 17. Circularity maps for robot KR 270-2 R2700, mm.

Fig. 18 .

 18 Fig. 18. Circularity maps for robot KR 360 L240, mm.

Fig. 19 .

 19 Fig. 19. Circularity maps for robot KR 500, mm.

Fig. 21 .

 21 Fig. 21. Experimental setup: milling circular trajectory in aluminum by Kuka KR 270 robot.

Fig. 23 .

 23 Fig. 23. Measured compliance errors along the machining trajectories.

Fig. 24 .

 24 Fig. 24. Measured compliance errors: polar plots (polar radius represents the compliance error, polar angle represents the point location on the benchmark circular trajectory).

Table 1

 1 Comparison of CNC machins and robots for machining application.

	Indicator	CNC machine	Industrial robot
	Accuracy	~0.005 mm	~0.1-1.0 mm
	Repeatability	~0.002 mm	~0.03-0.3 mm
	Workspace	Limited	Large
	Workspace extending	Impossible	Possible by adding extra actuated axis
	Kinematic architecture	Cartesian	Serial
	Number of actuated axes	3 or 5	6+
	Kinematic redundancy	Non	Yes, 1 dof at least
	Complexity of trajectory	Suitable for 3/5 axes machining	Any complex trajectory
	Relation between actuated and operational space	Linear	Non-linear
	Actuator feedback	Single encoder	Single or double encoders
	Mechanical compliance	Relatively high	Relatively Low
	Compliance error compensation	Non-required	Mechanical (Gravity compensators)
			Algorithmic (off-line and/or on-line)
	Dynamic properties	Moderate, homogeneous within the workspace	High, heterogeneous within the workspace
	Control algorithm	Continues path control	Point-to-point control
			Continues path control
	Programming language	Standardized G-code language	Manufacture Specified languages (KRL, V+, Karel, RAPID, Inform, etc.)
	Manufacturing flexibility	Single or several similar operations	Any type of operation
	Price	Competitive for 3-axis tools	Competitive for 6 dof robots
		Expensive for 5-axis tools

Table 2

 2 Principal geometric parameters of examined industrial robots, mm [ ].

	Robot	A	B	C	D	E	F
	KR100HA	350	750	1250	45	1000	210
	KR120 R3900	750	590	1350	41	1800	240
	KR270 R2700	350	675	1150	41	1200	240
	KR360 L240	500	1045	1300	55	1525	290
	KR500	500	1045	1300	55	1025	290

Table 3

 3 Examined robot specification data.

	Robot	Repeatability Workspace	Working	Payload
			volume	radius
	KR 100 HA 0.05 mm	45.95 m 3	2.6 m	100 kg
	KR120	0.06 mm	206.72 m 3	3.9 m	120 kg
	R3900
	KR270	0.06 mm	55 m 3	2.7 m	270 kg
	R2700
	KR360	0.08 mm	118 m 3	3.325 m	360 kg
	L240
	KR500	0.08 mm	68 m 3	2.825 m	500 kg

Table 4

 4 Stiffness parameters of examined robots.

	Robot	Equivalent joint compliances, μm/N
		k 1	k 2	k 3	k 4	k 5	k 6
	KR100 HA	1.92	0.28-0.55	0.56	3.31	3.83	5.42
	KR120 R3900	1.13	0.26-0.36	0.43	0.96	3.82	4.01
	KR270 R2700	0.54	0.28-0.30	0.42	2.79	3.48	2.074
	KR360 L240	0.86	0.11-0.28	0.25	2.17	1.47	2.96
	KR500	0.47	0.14-0.20	0.19	0.72	0.95	1.44

Table 5

 5 Circularity indices range of examined robots for elbow up/down configurations, mm.

	Robot	Under test loading		Under 80% of maximum
					loading
		min	max	middle	Min	max	middle
	KR 100	1.50/1.03 3.28/3.36 2.10/1.91 0.83/	1.82/	1.17/1.06
	HA				0.58	1.86
	KR120	4.00/2.54 6.66/6.57 5.08/4.05 2.67/	4.44/	3.39/2.70
	R3900				1.69	4.38
	KR270	1.36/0.84 3.05/3.13 1.93/1.64 2.04/	4.58/	2.90/2.46
	R2700				1.26	4.70
	KR360	1.69/1.02 2.85/2.81 2.29/1.80 3.34/	5.70/	4.58/3.60
	L240				2.04	5.62
	KR500	0.77/0.41 1.38/1.42 1.00/0.76 2.14/	3.83/	2.78/2.11
					1.14	3.94

Table 6

 6 Machining accuracy for optimally-placed square zones (elbow up/down configurations), mm.

	Robot	100×100 mm	200×200 mm	500×500 mm	1000×1000 mm
	KR 100 HA	1.60/1.12	1.73/1.22	2.24/1.56	-/2.40
	KR120 R3900	4.28/2.67	4.51/2.79	5.29/3.24	-/4.23
	KR270 R2700	1.47/0.90	1.60/0.97	2.06/1.23	-/1.87
	KR360 L240	1.81/1.08	1.96/1.15	2.39/1.50	-/2.01
	KR500	0.84/0.45	0.93/0.49	1.21/0.64	-/1.00

Table 7

 7 Machining accuracy for middle-placed square zones (elbow up/down configurations), mm. Evaluating robot performances with regard to principal indicators.

	Robot			100×100 mm		200×200 mm	500×500 mm	1000×1000 mm
	KR 100 HA			2.19/2.02			2.30/2.14	2.65/2.54	-/3.30
	KR120 R3900			5.21/4.10			5.35/4.21	5.81/4.57	-/5.28
	KR270 R2700			2.02/1.74			2.11/1.84	2.42/2.17	-/2.81
	KR360 L240			2.35/1.84			2.41/1.90	2.60/2.06	-/2.37
	KR500			1.05/0.81			1.11/0.86	-/1.02	-/1.33
	Fig. 20. Table 8
	Machining accuracy of examined robots under different task loading (with compliance
	error compensation), mm.
	Robots	Cutting force magnitude
		100 N	200 N	500 N	1000 N	2000 N	3000 N
	KR 100 HA	0.07	0.14	0.36	0.72	--
	KR120 R3900	0.18	0.35	0.88	1.76	--
	KR270 R2700	0.06	0.12	0.29	0.58	1.17	-
	KR360 L240	0.07	0.14	0.35	0.71	1.42	2.13
	KR500	0.03	0.06	0.14	0.28	0.57	0.85

Table 9

 9 Capability of examined robots for different machining tasks (with compliance error compensation).

	Desired	Cuttind force magnitude
	accuracy
		100 N	200 N	500N	1000N	2000N	3000N
	0.05 mm	KR 100 KR 100	KR 500
		KR 120 KR 270
		KR 270 KR 360
		KR 360 KR 500
		KR 500
	0.1 mm	KR 100 KR 100	KR 100	KR 500
		KR 120 KR 120	KR 270
		KR 270 KR 270	KR 360
		KR 360 KR 360	KR 500
		KR 500 KR 500
	0.2 mm	KR 100 KR 100	KR 100	KR 100	KR 500	KR 500
		KR 120 KR 120	KR 120	KR 270
		KR 270 KR 270	KR 270	KR 360
		KR 360 KR 360	KR 360	KR 500
		KR 500 KR 500	KR 500
	0.5 mm	KR 100 KR 100	KR 100	KR 100	KR 270	KR 360
		KR 120 KR 120	KR 120	KR 120	KR 360	KR 500
		KR 270 KR 270	KR 270	KR 270	KR 500
		KR 360 KR 360	KR 360	KR 360
		KR 500 KR 500	KR 500	KR 500
	1.0 mm	KR 100 KR 100	KR 100	KR 100	KR 270	KR 360
		KR 120 KR 120	KR 120	KR 120	KR 360	KR 500
		KR 270 KR 270	KR 270	KR 270	KR 500
		KR 360 KR 360	KR 360	KR 360
		KR 500 KR 500	KR 500	KR 500

Table 10

 10 Comparison of simulation and experimental results: circularity for different workpiece locations and cutting depth.

	Workpiece location	Cutting depth, mm		Cutting force, N (estimation)	Circularity, mm
		Desired	Measured		Estimation	Experiment
	Location #1	6.00	5.67-6.51	899.3	0.65	0.69
	Location #1	8.00	7.50-8.60	1199.1	0.86	0.89
	Location #1	10.00	10.00	1498.9	1.08	1.07
	Location #2	6.00	5.66-6.77	899.3	1.08	0.82
	Location #3	6.00	4.63-5.94	899.3	1.62	0.80

Acknowledgments

The work presented in this paper was partially funded by the project FEDER ROBOTEX, France (№ 38444).