# Variable Stiffness Identification and Configuration Optimization of Industrial Robots for Machining Tasks – Chinese Journal of Mechanical Engineering

#### ByJiachen Jiao, Wei Tian, Lin Zhang, Bo Li, Junshan Hu, Yufei Li, Dawei Li and Jianlong Zhang

Sep 6, 2022

Figure 9 shows the experimental platform of the study. A KUKA KR500 industrial robot was used as the operation carrier in the identification and machining tasks. An ATI IP60 Omega160 force transducer, which was fixed on the flange, was used to measure the 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 selected as the side lengths of the cubes to study the effects of different grid sizes. Therefore, the calibration space was divided into 2, 16, and 128 cubic grids, respectively. Based on this, the initial configuration of the samplings could be determined. When rotating the EE at ± 10° along the y-axis of the tool coordinate system, a total of 27 sampling configurations were obtained in a cubic space.

The joint stiffness without space gridding was identified as

$$begin{gathered} {user2{K}}_{theta } = {text{diag}}(1.58 times 10^{10} ,6.12 times 10^{9} ,5.28 times 10^{9} ,; hfill \ 4.66 times 10^{8} ,2.19 times 10^{8} ,3.49 times 10^{8} ) , {rm {N cdot mm/rad.}} hfill \ end{gathered}$$

(38)

Dividing the calibration space into two symmetrical cubic grids (Figure 10(a)), the joint stiffness in two grids, namely Grids 1 and 2, were identified as

$$begin{gathered} {user2{K}}_{theta G1} = {text{diag}}(1.48 times 10^{10} ,5.57 times 10^{9} ,6.80 times 10^{9} ,; hfill \ 3.38 times 10^{8} ,1.20 times 10^{8} ,2.17 times 10^{8} ) ,{rm {N cdot mm/rad,}} hfill \ end{gathered}$$

(39)

$$begin{gathered} {user2{K}}_{theta G2} = {text{diag}}(1.50 times 10^{10} ,6.29 times 10^{9} ,4.56 times 10^{9} , hfill \ 2.79 times 10^{8} ,3.28 times 10^{8} ,1.61 times 10^{8} ) ,{rm { N cdot mm/rad.}} hfill \ end{gathered}$$

(40)

By selecting a 300 mm cubic grid as a sampling unit, the calibration space was divided into 16 grids (Figure 10(b)). The entire space was observed as four cuboid spaces: Cuboids 1 to 4, whose longest sides were parallel to the y-axis of the base coordinate system. The joint stiffness of each grid space was identified, and the fluctuation of values was observed according to these four cuboid spaces (Figure 11).

A 150 mm cubic grid was selected as a sampling unit and the calibration space was divided into 128 grid spaces (Figure 10(c)). The entire space was divided into four cuboid spaces, i.e., Cuboids 1 to 4, which could be referred to as first-order cuboid spaces. Each first-order cuboid space was further split into four second-order cuboid spaces, such as Cuboids 1.1 to 1.4, which were formed with eight grid spaces. The joint stiffness of each grid space could be calculated, and the variation regularity of the stiffness value could be investigated based on the second-order cuboid spaces (Figures 12, 13, 14, 15).

Analyzing the variation trend of the joint stiffness in Figures 11 and 12, 13, 14, 15, the stiffness of the second, third, and fifth joints largely maintained the same tendency by dividing the robot calibration space into smaller grid spaces. However, the stiffness of the first, fourth, and sixth joints exhibited different change trends, compared with the tendency with larger grid spaces. Thus, the variation trend with larger grids cannot indicate the local change in joint stiffness, which inevitably reduces the accuracy of the stiffness model of robots.

In addition, the joints whose axes were parallel to the y-direction of the base coordinate system, namely the second, third, and fifth joints, exhibited stable change trends of stiffness in the entire calibration space. However, the joints whose axes were perpendicular to the y-direction of the base coordinate system, namely the first, fourth, and sixth joints, had similar stiffness tendencies in the calibration space. We can conclude 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 the variable stiffness identification method was verified through compensation experiments of the load-induced positional error. In the calibration space, a sampling position was randomly selected in each 150 mm grid. Thus, 128 verification points were obtained. A load of 50 kg was fixed to the EE to produce position errors.

The load-induced positional errors were measured by executing control commands before and after loading. According to Eqs. (1) and (15), the load-induced positional errors were calculated as follows:

$${user2{d}} = {user2{K}}_{fd}^{ – 1} {user2{f}} = {user2{E}}left( {{user2{JK}}_{theta }^{ – 1} (j){user2{J}}^{{text{T}}} } right){user2{f}},$$

(41)

where E is a matrix composed of the first three columns of the first three rows of the objective matrix. By combining the working load and the identified joint stiffness, the load-induced positional errors were predicted, and the control commands were modified through reverse compensation. Subsequently, the compensation effect of load-induced errors was evaluated by executing the modified commands under the loaded state. Defining Ex(i), Ey(i), and Ez(i) as the load-induced positional errors of the ith position on the x-, y-, and z-axes of the base coordinate system, the load-induced absolute errors E(i) were calculated using the following equation:

$$E(i) = sqrt {left( {E_{x} (i)} right)^{2} + left( {E_{y} (i)} right)^{2} + left( {E_{z} (i)} right)^{2} } .$$

(42)

The absolute positional errors before and after compensation with different grid sizes are shown in Figure 16, and Figure 17 shows the error distribution for different grids. The average absolute positional error induced by the 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 decreased 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 a 150 mm grid improved by approximately 52.54%. In conclusion, the validity of the identification results was verified, and the variable parameter error could better characterize the error model of a robot and obtain better positional accuracy.

### Configuration Optimization and Smooth Processing Experiments

The operation configuration optimization strategy and smooth processing method were verified through the simulation layout of a robot operation system in the DELMIA software. Using a workblank fixed on the tooling as the processing object, the machining path was planned in the software (Figure 18). The path from Tag 1 to Tag 2 was parallel to the y-axis of the base coordinate system, and the path from Tag 3 to Tag 4 also satisfied this scenario. According to the requirements of not using singular and joint-limit configurations, the rotation angle range of the robot EE was selected as − 90° ≤ θx ≤ 90°.

In the verification experiments, the normal stiffness of the robot motion trajectory in its work plane, defined as kv, and the tool axial stiffness kx, which were directly related to the trajectory accuracy and surface cutting quality, respectively, were considered as the evaluation indexes. The normal stiffness kv was calculated using the following equation:

$$k_{v} = sqrt {left( {k_{z} {text{cos}}theta_{x} } right)^{2} + left( {k_{y} {text{sin}}theta_{x} } right)^{2} } .$$

(43)

Thus, the fluctuations of kx and kv with different rotation angles (θx) in Tags 1 to 4 are shown in Figures 19, 20, 21, 22. Specifically, all the rotation angles corresponding to the configuration with optimal axial stiffness in Tags 1 to 4 were θx = 0°, that is, the initial robot configuration. The fluctuation of kv in Tags 1 and 4 formed a trough with θx = 0°, and two peaks were formed with θx = − 20° and θx = 10°. The optimal kv was observed at θx = − 20°. Similarly, two peaks formed in Tags 2 and 3 at θx = − 10° and θx = 20°, respectively, and the optimal kv was observed at θx =20°.

With the rotation angle farther from θx corresponding to the configuration with optimal stiffness performance, the stiffness in the target direction gradually decreased. Note that when the EE rotated away from the base coordinate system, the stiffness value decreased gradually.

Figures 23 and 24 compare kx and kv at different positions, respectively. The four positions were symmetrically distributed along the y-axis of the robot base coordinate system, and the stiffness change curves at the corresponding positions also exhibited spatial symmetry along the y-axis of the base coordinate system. In addition, the stiffness value exhibited a downward trend along the positive x-axis direction of the base coordinate system. In other words, under the same operational posture of the EE, kx or kv in Tag 1 was lower than the corresponding directional stiffness in Tag 3, and kx or kv in Tag 2 was lower than the corresponding directional stiffness in Tag 4.

To increase the trajectory accuracy in the operation process, the operation configuration with optimal kv was selected first. Based on the optimized configuration of Tags 1 to 4, the smoothing processing of the trajectory interpolation point was conducted to achieve the rapid acquisition of the optimal stiffness configuration of corresponding positions and smoothing results (Figure 25).

The smooth processing strategy simplified the configuration optimization of all target positions into several geometric feature positions, which significantly increased the configuration optimization efficiency. The optimization effect of smooth processing was evaluated by comparing it with the optimal stiffness performance obtained from the configuration optimization of the interpolation points. The interpolation positions of Tags 3 and 4 were selected, and the seven interpolation positions were defined as Tags 5 to 11. The kv obtained using the configuration optimization method and smooth processing are shown in Table 1, and Figure 26 shows a comparison between the smooth processing and configuration optimization.

As Table 1 shows, the stiffness index kv after smooth processing was highly similar to the result of configuration optimization, and the stiffness loss after fairing was less than 0.409% (in Tag 8). The stiffness of the target direction in Tags 5 and 11 after smoothing was slightly higher than that of configuration optimization. This was because the step size of the rotation angle was larger in the process of optimization, which resulted in the optimized configuration obtained by smoothing closer to the configuration with optimal stiffness.

In summary, smooth processing effectively increases the optimization efficiency of a robot operating configuration and ensures that the improvement effect of the axial stiffness through smooth processing is not significantly reduced compared with that of configuration optimization, which fully proves the effectiveness of the fairing method.

### Machining Experiment

The effects of configuration optimization and smoothing on milling quality were studied using a robot machining system. In the verification experiment, the milling task was performed in the cylinder head of an automobile engine. The planned trajectory in the machining experiment is shown in Figure 27, where Tags 1 and 2 represent the start and end positions of the trajectory, respectively.

By selecting the tool axial stiffness as the optimization objective, the rotation angles corresponding to the start and end positions of the machining trajectory were obtained using the configuration optimization method; the optimized configurations of Tags 1 and 2 are shown in Table 2, and smooth processing of the interpolation positions was conducted to obtain optimal configurations. Based on this, the trajectory control program of the robot was generated through offline programming software, and the robot could be driven to perform milling tasks. The process parameters of the robot milling experiment are listed in Table 3.

The comparison results of the milling process after configuration optimization are shown in Figure 28. The milling surface quality before configuration optimization was relatively poor, and the blades severely vibrated. The roughness of the surface was Ra 2.356, as measured using a roughness meter (Sanfeng SJ-210). After optimizing the operation configuration of the robot to increase the end operation rigidity, the obtained milling surface of the workpiece was relatively smooth, the surface roughness was Ra 0.597, and the accuracy and quality of milling were significantly improved.