Variable stiffness identification and configuration optimization of industrial robots for machining tasks

 Abstract: Industrial robots are increasingly used in machining tasks because of their high flexibility and intelligence. However, the low structural stiffness of the robot seriously affects the positional accuracy and machining quality of robot operation equipment. Studying robot stiffness characteristics and optimization methods is an effective way to improve robot stiffness performance. Accordingly, aiming at the poor accuracy of stiffness modeling caused by approximating stiffness of each joint as constant, a variable stiffness identification method is proposed based on space gridding. Then, a task-oriented axial stiffness evaluation index is proposed to realize quantitative assessment of the stiffness performance in the machining direction. Besides, by analyzing the redundant kinematic characteristics of the robot machining system, a configuration optimization method is further come up with to maximize the index. For a large number of points or trajectory processing tasks, a configuration smoothing strategy is proposed to achieve fast acquisition of optimized configurations. Finally, experiments on a KR500 robot are conducted to verify the feasibility and validity of proposed stiffness identification and configuration optimization methods.


Introduction
The application of industrial robots has significantly increased in automatic manufacturing, due to their remarkable advantages of great operation flexibility, high intelligence, low cost and low space requirement [1][2][3]. Unfortunately, the stiffness of industrial robots is only 2%  Wei Tian tw_nj@nuaa.edu.cn 1 Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 2 Beijing Institute of mechanical equipment, Beijing 100854, China to 5% of that of CNC machines, because of the inherent characteristics of their series structures [4][5][6]. Such weak stiffness of industrial robots seriously affects the robot positional accuracy and the processing quality of products, which severely limit their application in the manufacturing and assembly of high value-added products. Therefore, it is essential to deeply research the stiffness characteristics and strengthening methods to improve the operation performance of robot equipment. Many researchers have investigated robot stiffness characteristics and modeling methods. Dumas and Caro proposed a fast identification method of joint stiffness, and established the Cartesian stiffness model of KUKA KR240-2 robot through experimental results [7]. Zhou et al. identified the geometric parameters to revise the robot kinematics model and Jacobian matrix. On this basis, more precise joint stiffness matrix and the stiffness model could be obtained [8]. Klimchik et al. adopted the method of virtual axis to couple the connecting rod gravity and the external load into the stiffness modeling process, which obtained higher stiffness modeling accuracy [9]. The aforementioned stiffness identification and modeling methods are established based on assumptions of rigid connecting rods and elastic rotation axis of joint. However, the gravity center of the connecting rod，and the torsional stiffness of motor and reducer vary with the operation configuration. Besides, different sampling configurations will result in different identification results, because of the inherent nonlinearity of joint stiffness [10]. Namely, a set of joint stiffness identification result cannot satisfy the accuracy requirements of the stiffness model in the whole workspace well.
Furthermore, the stiffness performance of the robot mainly depends on three factors: (1) connecting rod structure and material properties; (2) stiffness of actuator and transmission mechanism; (3) working configuration. Comparing to first two factors, the robot configuration Variable stiffness identification and configuration optimization of industrial robots for machining tasks ·3· optimization has the advantages of high technical feasibility, good task adaptability, and no need to change the structure and control system of the robot, which shows high application value in engineering application.
Sabourin et al. developed a comprehensive optimization objective function combining kinematics performance and stiffness performance, so as to realize the configuration optimization of robot machining system [11]. Xiong et al. proposed a discrete search algorithm for the operation configuration of milling robot with optimal stiffness performance, and the optimization results significantly improved the trajectory accuracy [12]. Guo et al. took the robot stiffness ellipsoid as the stiffness evaluation index to improve the robot stiffness by configuration optimization, and obtained accuracy of higher drilling axial and countersink depth [13]. However, these studies focus on configuration optimization for processing tasks with discrete positions. When dealing with a large number of points or trajectory tasks, such as thousands of target positions in drilling tasks of aerospace components, the optimization of target points one by one seriously restricts the processing efficiency, which cannot meet the requirement of short cycle manufacturing tasks.
Motivated by these points, a variable stiffness identification method was investigated to realize accurate modeling of industrial robots, which overcomes the influence of configurations on stiffness modeling. Besides, a configuration smoothing method is raised to realize fast acquisition of optimized configurations, based on the geometric feature positions of target machining tasks.
The remainder of this paper is organized as follows. Section 2 puts forward a planning method of sampling point based on space gridding, on this basis, a variable stiffness identification method is put forward. In Section 3, a task-oriented stiffness representation model is established to estimate the stiffness performance in a given direction, and a configuration optimization strategy is come up with to maximize the performance. Aiming at the a large number of points or trajectory tasks, Section 4 proposes a robot configuration smoothing processing method to achieve fast acquisition of optimized configurations. In Section 5, the verification experiments of identification results and configuration optimization methods are designed and discussed. Finally, Section 6 concludes the paper.

Variable Stiffness Identification Based On Space Gridding
The robot Cartesian stiffness is directly dependent on joint stiffness and robot configuration according to traditional static stiffness model. However, the joint stiffness identification and robot configuration are interrelated to each other. On the one hand, stiffness of each joint varies with robot configuration, because of the complex transmission structure, which mainly consists of motor and reducer [14]. On the other hand, different configurations lead to different structural gravity centers and kinematics parameter deformations, which effect the identification accuracy of joint stiffness [15]. Thus, variable stiffness identification is an effective method to overcome the impact of configuration on stiffness identification results.

Stiffness Identification Principle
Based on the elastic deformation assumption, the relationship of displacement of end effector (EE) and corresponding working load could be written as: where F is the wrench vector acting on the EE, and D denotes displacement vector. C is the inverse matrix of stiffness matrix K, namely robot compliance matrix. Obviously, the deformation caused by working load is nonlinear with the stiffness of each joint. Thus, when identifying the joint stiffness of in experiments, it is necessary to carry out tremendous sampling experiments. Multiple groups of static loads under different robot configurations should be applied, and the displacements at the end of the robot could be measured in turn. Simplify the right side of the Eq.(1), the relation of D and joint compliance could be shown as follows: (2) where F, A and C θ are defined as follows:

·4·
Based on the concept of generalized inverse matrix, C θ could be calculated: According to the identification method described above, the quantity and distribution of the chosen sampling points, which relate to robot configurations, directly influence the accuracy of identification results [16,17]. Hence, the fitting result of joint stiffness in the whole workspace cannot accurately reflect the performance of robots in different subspaces [18]. Therefore, it is necessary to find joint stiffness parameters suitable for different workspaces to improve stiffness characterization accuracy.

Space Gridding Principle Analysis
If the rotation angle and corresponding stiffness of each joint is determined, the robotic stiffness model can be established. Define K θ1 as the joint stiffness in a given configuration, and K θ2 is the joint stiffness in a close configuration: where E is the 2-norm of the difference between K θ1 and K θ2 . A positive number ξ approaching zero always exists which satisfies E<ξ, with the angle variation △θ closes to zero. Therefore, the corresponding joint stiffness of two close configurations can be considered highly similar to each other [19,20]. Through above analysis, the workspace of robot could be divided and the joint stiffness corresponding to each subspace could be considered as fixed values approximately. Due to the advantage of higher processing efficiency, the reachable space of robot can be divided into a set of cubic grids, which could be shown as Figure 1(a).
Eight vertices of given grid could be chosen as sampling points to envelop a single grid. Besides, the center point of the grid can be chosen as a sampling point to ensure that the stiffness parameter in the grid is closer to the real value. Thus, nine target points (Tag 1 to Tag 9) are used as sample positions to identify joint stiffness corresponding to a given grid space, as shown in Figure  1(b). In the initial sampling posture of a target position, the three axes of the tool coordinate system are consistent with the directions of the three intersecting edges of the grid, such as the sampling posture of Tag5. Besides, different sampling postures can be obtained by rotating a certain angle around any axis of the tool coordinate system. Figure  1(c) shows the generation of a new posture by rotating the y-axis of the tool coordinate system. On this basis, different sampling postures of nine sampling positions could be obtained, as shown in Figure 1(b) [21]. Moreover, the joint stiffness tends to be stable with more than 10 groups sampling configurations in a stiffness identification experiment [22][23][24]. Therefore, in each sampling position of a grid space, more than two configurations should be used to ensure the accuracy of stiffness identification results.

Variable Stiffness Identification
Combining the static stiffness identification method and space gridding principle, a variable stiffness identification method is innovatively come up with. Sampling positions could be obtained on the basis of space gridding, and the robot stiffness model, which suits the whole workspace, could be built precisely and described as: where j is the grid serial number, and () j  K is the joint stiffness corresponding to the grid space where the current robot's target position is located. The Eq.(9) could be denoted as the variable stiffness model. Compared with the traditional static stiffness model, the variable stiffness model can more accurately reflect the stiffness characteristics of different robot intervals, thereby achieving accurate evaluation of stiffness performance and precise prediction of load-induced positional error.
The steps to variable stiffness identification method could be shown as Figure 2.
According to the identification results, joint stiffness of the grid space where the target position is located can be obtained. However, when the location of a target point happens to be in some special positions as shown in Figure  3, choosing which joint stiffness is a confusing problem.
Variable stiffness identification and configuration optimization of industrial robots for machining tasks   When located in a vertex, side or plane that belongs to just one grid, the joint stiffness in this grid space could be used. In addition, a target position may also appear in the following special locations: (1) in a common vertex of two, four or eight grids; (2) in a common line of two or four grids; (3) in the contact plane of two grids. When the grid size is small enough, the aforementioned situations could be dealt with by averaging each joint stiffness in these grids to which the target point belongs.
The proposed stiffness identification method is an innovation or improvement on the static stiffness model. Therefore, although this paper takes KUKA KR500 as a case to discuss, but the method could be utilized for other kinds of robots. On the basis of accurate joint identification and modeling of robot, the evaluation and optimization methods of robot stiffness performance could be studied to realize the practical application of stiffness model in engineering.

Task-oriented Axial Stiffness Performance Index
In the robot machining process, the working load will cause the slip deformation of the robot, where the tool radial deformation affects the cutting positional accuracy

·6·
(drilling, boring, etc.) or the trajectory accuracy (milling, grinding, etc.), the tool axial deformation mainly affects the cutting surface quality (roughness, flatness, etc.). Therefore, the definition of the robot stiffness performance index should be adapted to the load distribution characteristics of the robot's specific task. Based on the different dimensions of elements, Eq.(1) could be transformed as follows: (10) where f and m indicate the force and torque matrix respectively; d and δ represent the translational and rotational displacement matrix. The posture error caused by cutting load is very small, and the on-line detection method is difficult to track the working posture of end effector, which can be ignored. Then: In addition, the main influencing factor of machining quality is the linear displacement of the end of the robot caused by the cutting load, and the change of the end posture caused by the torque has little effect on the machining accuracy and quality [25,26]. Therefore, Eq.(11) could be similarly simplified as: ,0 Apply a unit force f to the EE of robot, there is: Eq.(13) describes an ellipsoid that changes with robot operation configuration, which could be shown in Figure 4, whose directions of principal axes are eigenvectors of During the process of robot equipment operation, the direction of machining load at the end of the robot is not consistent with the main axis of the robot stiffness ellipsoid. Therefore, it is necessary to research the axial stiffness performance of robot to realize the accurate evaluation of machining performance. In Figure 4, λ t1 and λ t2 represent the two half axes of the elliptical section of the machining plane in the envelope space of the stiffness ellipsoid. λ d is the semiaxis of the elliptical section's normal vector. The square root of the three half axis length respectively represents the stiffness values of the three axes in the robot tool coordinate system [27,28].

Figure 4 Stiffness ellipsoid
The direction of λ d in the stiffness ellipsoid coordinate system could be defined as the unit vector [e x e y e z ], and mathematical formula of the line where it is located is expressed as follows: Thus, stiffness ellipsoid can be written according to the definition formula of ellipsoid as: where x, y and z respectively represent the coordinates of the intersection point between the straight line shown in Eq. (14) and the stiffness ellipsoid. Combining Eq. (14) and Eq. (15), λ d and k x can be calculated as below: The method also could be extended to calculate k y and k z , which could be show as follows:

Machining Plane
Variable stiffness identification and configuration optimization of industrial robots for machining tasks

·7·
where [h x h y h z ] and [l x l y l z ] are unit vectors of λ t1 and λ t2 respectively. Through these three formulas, the three-dimensional stiffness in the robot tool coordinate system can be accurately calculated, which provides a theoretical basis and evaluation standard for the precision control and quality improvement of machining tasks.

Stiffness Oriented Configuration Optimization Method
Redundant degrees of freedom (DOF) make the robot operating system have higher task adaptability. Consequently, the tasks of the robot operation system can be divided into two levels of subtasks. The first level subtask controls the EE move to the initial position and posture to meet the requirements of the operation task. As the second level subtask, configuration optimization improves the machining performance by using redundant DOFs [29]. The 6-DOF industrial robot has a functional redundant DOF in the direction of the tool axis, no matter position processing tasks (drilling, boring, etc.) or trajectory processing tasks (milling, grinding, etc.). Namely, there are infinite number of poses of EE can be used to execute machining tasks theoretically [30]. Therefore, change step of rotation angle should be set to promote the running speed of optimization method.
The feasible posture of EE is described in Cartesian coordinate system. However, the stiffness performance index of robot are evaluated in joint space. Therefore, it is necessary to use the inverse kinematics algorithm and the uniqueness principle to solve the joint angles corresponding to the current robot configuration.
The specific steps of robot configuration optimization are as follows: (1) Extract the position coordinates P (p x , p y , p z ) and initial posture N (α, β, γ) of the target point from process digital model in the off-line programming software. The current configuration the robot can be obtained by using the inverse kinematics algorithm and the uniqueness principle of the inverse solution.
(2) Use the axial stiffness evaluation index to calculate the target stiffness value in the initial configuration, denoted as k.
(3) Rotate the EE along the tool axis with θ x , where -180° ≤ θ x ≤ 180°, and defines the variable step of θ x as 10 x     . Thus, θ x could be calculated as follows: (4) Examine safety and rationality of the new configuration. If the posture of EE is unreachable or interfered, then abandon the configuration. Otherwise, calculate the current configuration  C and corresponding stiffness performance index k , and update k through the bigger one of k and k .
(5) Repeat the third step and fourth step to obtain the rotation angle with the best axial stiffness, and the corresponding optimal configuration could be determined.
The workflow of configuration optimization is shown as Figure 5.

Smooth processing strategy
In aviation manufacturing industry, any part has a large number of drilling demand. If the operation configuration of each target machining position was optimized one by one, the workload of optimization processing is large and the efficiency is low. Therefore, for processing tasks with a great quantity of discrete positions, according to the characteristics of the target product geometry or position distribution, taking the typical operation position as the ·8· optimization objective, using interpolation smoothing method to optimize the robot configuration of other target positions is an efficient processing strategy to improve the optimization efficiency. Similarly, this method can be extended to milling and other trajectory cutting tasks, which is controlled by interpolation points actually. Typical operation points could be shown in Figure 6.

Figure 6 Several typical operation points
In addition to the special case that the workpiece surface is plane, the direction of each axis in the tool coordinate system of each machining position cannot be guaranteed to be in the same direction or in the same plane. Therefore, it is necessary to unify all target positions to the same reference coordinate system. In this paper, the base coordinate system is selected as the reference.
The rotation matrix of posture of EE can be expressed as RPY angle: Since Base i R is an identity invertible matrix, then Then a unit vector I could be got as: (24) Adjust the structure of Eq.(24) and get the expression of Rot , Rot , Rot , , , Base ix R represents the coordinate system obtained by rotating the robot base coordinate system a certain angle around the x-axis of the robot, whose RPY angle could be expressed as    (27) Therefore, the corresponding posture of the target position can be solved in proportion to the distance from the start position to the end position, so as to achieve the purpose of posture smoothing.
For point machining tasks, select the appropriate reference coordinate axis (y-axis or z-axis direction) according to the position coordinate distribution law of the point, and theoretical deflection angle of rotation matrix corresponding to any machining position could be obtained as Eq. (28)  Variable stiffness identification and configuration optimization of industrial robots for machining tasks ·9· deflection angle i   corresponding to any interpolation position in the trajectory processing operation can be expressed as: Therefore, for any intermediate position, the adjustment angle obtained from smooth processing can be written as:  Figure 8 Experimental system Figure 8 shows the experimental platform of the study. A KUKA KR500 industrial robot was used as the operation carrier in identification and machining tasks. An ATI IP60 Omega160 force transducer fixed on the flange, was used to measure working loads. An API laser tracker was used to establish coordinate systems and measure the positional errors of the robot.

Experimental Verification of Variable Stiffness Identification
In the stiffness identification experiment, a 1200 mm × 600 mm × 600 mm cuboid was planned as the calibration space, and 600 mm, 300 mm and 150 mm were chosen as side lengths respectively to study the influence of different grid sizes. Therefore, the calibration space could be divided into 2, 16 and 128 cubic grids respectively. On this foundation, the initial configuration of the samplings could be decided. Rotating the EE with ±10° along the y-axis of tool candidate system respectively, and a total of 27 sampling configurations could be got in a cubic space. The joint stiffness without space gridding could be identified as: 10 Dividing the calibration space into 2 symmetrical cubic grids as shown in Figure 9(a), joint stiffness in two grids, namely Grid 1 and Grid 2 could be identified as: Choosing 300 mm × 300 mm × 300 mm cubic grid as a sampling unit, the calibration space can be divided into 16 grids, as shown in Figure 9(b). The entire space could be seen as four cuboid spaces, including Cuboid 1 to Cuboid 4, whose longest sides are parallel to the y-axis of base coordinate system. The joint stiffness of each grid space can be identified and the fluctuation of values could be observed on the basis of these four cuboid spaces, which could be shown as Figure 10.
Choosing 150 mm × 150 mm × 150 mm cubic grid as a sampling unit and the calibration space can be divided into 128 grid spaces, which could be shown in Figure 9(c). The entire space could be divided into four cuboid spaces, including Cuboid 1 to Cuboid 4, which could be called as first-order cuboid spaces. A first-order cuboid space could be once more split into four second-order cuboid spaces, such as Cuboid1.1 to Cuboid 1.4, which are formed with eight grid spaces respectively. The joint stiffness of each ·10· grid space could be calculated and the variation regularity of stiffness value could be researched based on the second-order cuboid spaces, which could be shown as  Variable stiffness identification and configuration optimization of industrial robots for machining tasks

Figure 15
Compensation results of different size grids Analyzing the variation trend of joint stiffness in Figure  10 and Figures 11-14, the stiffness of second, third and fifth joint, largely maintain the same tendency by dividing robot calibration space into smaller grid spaces. On the other hand, the stiffness of first, fourth and sixth joint displays different change trends, compared to the tendency with bigger grid spaces. Thus, the variation trend with bigger girds cannot show clearly the local change of joint stiffness, which inevitably reduces the accuracy of stiffness model of robots.
Besides, the joints whose axes are parallel to the y direction of base coordinate system, namely second, third and fifth joints, have stable change trends of stiffness in whole calibration space. However, the joints whose axes are perpendicular to the y direction of base coordinate system, namely first, fourth and sixth joint, have similar stiffness tendencies in the calibration space. It could be concluded that the axis direction of a robot joint is related to the stiffness distribution of the joint, which can be used to optimize the machining configuration of industrial robots.
The effectiveness of variable stiffness identification method was verified by compensation experiments of load-induced positional error. In calibration space, a sampling position was randomly selected in each 150 mm grid. Thus, 128 verification points were got in total. A load of 50kg was fixed on the EE to produce position errors.
The load-induced positional errors could be measured by performing control commands before and after loading. According to Eq.(12), the load-induced positional errors could be calculated as follws: where D is a matrix composed of the first three columns of the first three rows of the objective matrix. Combining the working load and the identified joint stiffness, the load-induced positional errors were predicted and control commands were modified through reverse compensation. Then, the compensation effect of load-induced errors could be evaluated by executing modified commands under the loaded state. Defining E x (i), E y (i) and E z (i) as the load-induced positional errors of i-th position on x, y and z axis of base coordinate system, and the load-induced absolute errors E(i) could be calculated through the equation: Variable stiffness identification and configuration optimization of industrial robots for machining tasks ·13· Figure 15, and Figure 16 shows the error distribution with different grids. The average absolute positional errors induced by working load was 0.2868 mm, and the maximum error was 0.3587 mm. As the number of space grids increased from 1 to 128, the average value of load-induced positional errors after compensation reduced from 0.1201 mm to 0.0570 mm, and the maximum error decreased from 0.1610 mm to 0.1134 mm. The compensation effect with 150 mm grid could be improved by nearly 52.54%. In conclusion, the validity of identification results could be verified, and variable parameter error can better characterize the error model of the robot and obtain better positional accuracy.

Configuration Optimization and Smooth Processing Experiments
The operation configuration optimization strategy and smooth processing method are verified through simulation layout of robot operation system in DELMIA software. Taking the blank fixed on the tooling as the processing object, the machining path is planned in the software as shown in the Figure 17. The path from Tag1 to Tag2 is parallel to the y-axis of the base coordinate system, the path from the Tag3 to Tag4 also satisfies this situation. According to the requirements of keeping away from singular configurations and joint-limit configurations, the rotation angle range of the robot EE could be chosen as -90 ° ≤ θ x ≤ 90 °.

Figure 17 Simulation environment and machining trajectory
In the verification experiments, the normal stiffness of the robot motion direction in the work plane, defined as K v , and the tool axial stiffness K x , which are directly related to the trajectory accuracy and surface cutting quality separately, where K v is calculated using the equation: Thus, the fluctuation situation of K x and K v with different rotation angle θ x in Tag1 to Tag4 could be shown as Figures 18 -21. Specifically, all rotation angle corresponding to the configuration with optimal axial stiffness in Tag1 to Tag4 is θ x = 0°, that is, the initial robot configuration. The fluctuation situation of K v in Tag1 and Tag4 forms a trough with θ x = 0°, and two peaks are formed with θ x = -20° and θ x = 10°. The optimal K v could be found with θ x = -20°. Similarly, two peaks are formed in Tag2 and Tag3 with θ x = -10° and θ x = 20° respectively, and the optimal K v could be found with θ x =20°.
With the rotation angle far away from θ x corresponding to configuration with optimal stiffness performance, the stiffness in target direction gradually decreases. It is worth noting when the EE rotates away from the base coordinate system, the stiffness value decreases slower.

Figure 21
Stiffness fluctuation situation in Tag4 Figure 22 and Figure 23 compare the K x and K v in different positions. The four positions are symmetrically distributed along the y-axis of the robot base coordinate system, and the stiffness change curves at corresponding positions also show spatial symmetry along the y-axis of the base coordinate system. In addition, the stiffness value shows a downward trend along the positive x-axis direction of the base coordinate system. Namely, under the same operation posture of end effector, K x or K v in Tag1 is less than the corresponding directional stiffness in Tag3, and K x or K v in Tag2 is less than the corresponding directional stiffness in Tag4.  In order to improving the trajectory accuracy in trajectory operation process, the operation configuration with optimal K v can be selected first. On the basis of optimized configuration of Tag1 to Tag4, the smoothing processing of the trajectory interpolation point is carried out to achieve the rapid acquisition of the optimal stiffness configuration of corresponding positions, smoothing results as shown in Figure 24.

Figure 24 Smoothing results of interpolation positions
The smooth processing strategy simplifies the configuration optimization of the all target positions into several geometric feature positions, which can significantly improve the configuration optimization efficiency. The optimization effect of smooth processing could be evaluated by comparing with the optimal stiffness performance obtained from the configuration optimization of the interpolation points. Select the interpolation positions of Tag3 to Tag4, and define the 7 interpolation positions as Tag5 to Tag11. The K v of the configuration optimization method and smooth processing could be shown in Table 1, and Figure 25 shows the comparison between smoothing processing and configuration optimization. According to Table 1, it can be found that the stiffness index K v after smooth processing has a high similarity with the result of configuration optimization, and the stiffness loss after fairing is less than 0.409% (in Tag8) Variable stiffness identification and configuration optimization of industrial robots for machining tasks ·15· smoothing is even slightly higher than that of configuration optimization. The reason is that the step size of the rotation angle is larger in the process of optimization, which leads to the optimized configuration obtained by smoothing closer to the configuration with optimal stiffness.

Figure 25
Optimized posture comparison between smoothing processing and configuration optimization In summary, the smooth processing effectively improves the optimization efficiency of the robot working configuration, and ensures that the improvement effect of the axial stiffness through smooth processing is not significantly reduced, compared with that of configuration optimization one by one, which fully proves the effectiveness of the fairing method.

Machining Experiment
The effect of configuration optimization and smoothing on milling quality is studied by using robot operation equipment in Section 5. In the verification experiment, the milling task is performed in a cylinder head of automobile engine. The planned trajectory in the machining experiment could be shown as Figure 26, where Tag 1 and Tag 2 represent the start and end positions of the trajectory respectively.

Figure 26
Milling trajectory in the machining task Choosing the tool axial stiffness as optimization objective, the rotation angles corresponding to the start and end position of the machining trajectory were obtained using the configuration optimization method, the optimized configurations of Tag1 and Tag2 as shown in Table 2, and then smooth processing of the interpolation positions could be carried out to get optimal configurations. On this basis, the trajectory control program of the robot was generated through off-line programming software, and then the robot could be driven to carry out milling tasks. Process parameters of robot milling experiment could be shown as Table 3.   The comparison results of the milling process after configuration optimization could be shown in Figure 27. The milling surface quality before configuration optimization is relatively poor, and there were more serious vibration blades. The roughness of the surface is Ra 2.356 measured by roughness meter (Sanfeng SJ-210). After optimizing the operation configuration of the robot to improve the end operation rigidity, the obtained milling surface of the workpiece is relatively smooth, the surface Tag3  Tag5  Tag6 Tag7  Tag8  Tag9 Tag10 Tag11 Tag4 Configuration smooth processing Configuration optimization one by one Cylinder head Tag1  Tag2 Cutting direction ·16· roughness is Ra 0.597, the accuracy and quality of milling are significantly improved.

Conclusions
(1) A regular sampling point selection method is proposed on the basis of space gridding. On this foundation, considering joint stiffness as approximately constant values in any grid space, and a variable stiffness identification method is put forward. (2) A task-oriented axial stiffness evaluation index is raised to estimate the stiffness performance in a given direction, and a configuration optimization strategy is come up with to maximize the performance by utilizing the redundant DOF of robot equipment. (3) Aiming at the a large number of points or trajectory tasks, a robot configuration smoothing processing method is come up with to achieve fast acquisition of optimized configurations, which effectively improves the efficiency of attitude optimization. (4) Experiment of the stiffness identification was completed on KUKA KR500 system. The difference of joint stiffness in different spaces is verified by analyzing the identification results, and stiffness variation regularity of each joint is clearly got in calibration space. (5) The experiments of configuration optimization and smooth processing methods were tested in simulation environment. The experimental results showed that smooth processing strategy improves the optimization efficiency while the stiffness loss is very small. According to the machining results in a cylinder head of automobile engine, the milling quality has been improved obviously after the configuration optimization, and the validity of these methods are verified.