Calibration of a Manipulator With a Regularized Parameter Identification Method

As an inverse problem, the parameter identification of manipulators is essential to in-situ calibration. Since the ill-posed inverse kinematic model is sensitive to the measurement value, even tiny errors will make the geometric model of the manipulator wrongly identified. To overcome this problem, a Regularized Parameter Identification Method (RPIM) is proposed to calibrate the geometric model of the manipulator. The inverse kinematic model of a 6 DOF manipulator is modified by a Tikhonov regularization to overcome its ill-posed problem. The regularization parameter is optimized by the improved L curve method to adjust the initial model to a well-posed one that approximates the real situation. A calibration system is designed to evaluate the effectiveness of the suggested method. The position of the selected targets is tested by using a laser tracker. The experimental result shows that the absolute position errors of the manipulator are 2.533mm without calibration, 0.472mm by RPIM, 1.445mm by Least Square Method (LSM), 1.353mm by Gradient Descent (GD), and 0.956mm by Gauss-Newton (GN). It shows that the absolute position error of RPIM is reduced by 81.331% after calibration, which is superior to other methods.


I. INTRODUCTION
The collaborative manipulators are widely applied to station docking, extravehicular maintenance, intelligent manufacturing, vehicle manufacturing, and other industrial applications [1]. As a critical component of manipulators, the manipulator can ensure the actuator locates the target successfully. However, the positioning errors of manipulators restrict the development of intelligent manufacturing on a robotic production line. Due to structural wear and performance attenuation, as is known to all, positioning accuracy can be increased after manipulators work continuously for a long time. The location error is harmful, for example, it will increase the difficulty of automatic assembly for mechanical parts, the probability of mislocation on riveting holes, the welding error for the welding manipulator, and the control difficulty of a telemedicine cooperative manipulator [2], [4], [5], [6], [7].
The associate editor coordinating the review of this manuscript and approving it for publication was Guilin Yang .
The studies on accuracy improvement are essential to manipulator maintenance [8].
The regular in-situ calibration of the manipulator can improve the positioning accuracy of the manipulator by identifying the geometric parameters of the kinematic model. To identify the geometric parameters of a kinematic model, the pose relationship of manipulator manipulators should be described first by a model. There are many kinds of models, such as the D-H model, the SD-H model, the MD-H model, the S model, the CPC/MCPC model, the POE model [9], [10], [11], [12], and other parameterized models. Among which, the D-H model is the traditional method with a wide application. However, when the two axes of a manipulator are parallel or close to parallel, it will cause a problem of singularity. To overcome this problem, the MD-H model with a simple structure is introduced by introducing the y-axis rotation [13]. In other models are hard to overcome the illposed problem of the inverse kinematic model.
Parameter identification is an important part of calibration for manipulators. There are many methods for parameter identification. The least-square method is a commonly used identification method [14]. For example, ALBERT et al. applied the LSM to the system identification of a manipulator by considering the influence of noise, environment, experimental setup, and other interference of measurement [15], [16], [17]. Meggiolaro M et al. obtained the geometric parameter errors of a manipulator by utilizing the least square method (LSM) [18]. Tang and Mooring designed a composed board with precise positioning points to decrease the calibration error of a manipulator [19], [20]. Hage et al. proposed a constraint equation by using four planes to improve the accuracy of a probe-based manipulator [21]. The maximum likelihood estimation is applied by Renders [22] and other researchers to paratha meter identification of the manipulator kinematic geometric parameters, which achieved a good result. With the continuous improvement of the computing level, the Levenberg-Marquardt (LM) method [23], the extended Kalman filter method [24], the simulated annealing algorithm [25], the parameter optimization method, and other iterative algorithms are also applied to parameter identification of the complex system.
These studies have yielded positive results in optimizing manipulator performance in spite of certain limitations. For the LSM method, although the iteration process is simple, the convergence speed is fast, and its applications are extensive, this test method has the disadvantage of a large amount of calculation. Although the Gradient Descent method has a fast convergence speed, it cannot guarantee the globalist's optimality. The Gauss-Newton is easy to implement in spite of no logic to control the step size. The Extended Kalman Filter algorithm may suffer from target loss.
However, aiming at the ill-posed problem of the calibration process, these methods are hard to solve this problem. For this reason, researchers have made many explorations and achieved good results. For example, Kostenko [26] et al. applied the Laplace operator on a regularization matrix to invert seismic coseismal slip distribution inversion [27], [28]. Huang et al. proposed a regularization matrix method with expected properties to estimate the information in the medium [9]. H.Save and his teammates solved the ill-posed problem of the gravity field by using a regularization matrix [29]. Ditmar and his followers selected the first-order derivative regularization matrix and the Kaula regularization to reflect the statistical law of the potential coefficient [30]. Li et al. applied the regularization matrix with distance weighting to the abnormal detection of hyperspectral imaging [31]. Gauthier et al. proposed a definition of a regularization matrix in the acoustic domain [32], [33]. The central idea of the regularization matrix focuses on the targeted correction of the singular values. By constructing a regularization matrix, the small singular values will be corrected and the ill-posedness of parameter identification can be improved [34].
The parameter identification of manipulator manipulators is essential to in-situ calibration. Since the ill-posed inverse kinematic model is sensitive to the measurement value, even tiny errors will make the geometric model of the manipulator wrongly identified. To overcome this problem, a Regularized Parameter Identification Method (RPIM) is proposed to calibrate the geometric model of the manipulator. It is organized as follows. In Section 2, a kinematic model of a 6 DOF manipulator is built based on the M-DH model and geometric error model. In Section 3, the kinematic model is reconstructed by a regularization approach. In Section 4, both numerical and experimental studies are carried out on a calibration system to validate the accuracy and the effectiveness of our suggested method. In Section 5, several conclusions and research expectations about our work are given.

II. KINEMATIC MODEL OF A 6 DOF MANIPULATO A. KINEMATIC MODEL
According to the MD-H modeling method, the conversion process between each link of the 6 DOF manipulator is designed as follows.  With the help of translation matrix Trans(·) and rotation matrix Rot(·), the relative poses are changed as [35] where, A i is the transformation matrix of the i th joint frame. 90536 VOLUME 10, 2022 Then we can get the relationship between each two adjacent joints i-1 th and i th can be written as According to Eq.(2), the end pose of the manipulator can be solved as The geometric parameter errors of the manipulator are accumulated by each joint, which include the geometric parameter errors of a i−1 , β i , α i , t i , d i . By substituting those geometric parameter errors into Eq.(1), the geometric errors are where, δ i−1 i T is the differential expression between adjacent joints of (i-1) th and i th .
The geometric parameter errors are expressed as a differential form, then the transformation matrix can be written as Then the accumulative errors of the 6 DOF manipulator between each joint can be written as The redundant parameters should be eliminated before identification. By considering the redundancy theory and the rank of matrix, the Jacobian matrix U is obtained [36].
Then the relationship between position errors and geometric parameter errors is where S is position errors vector, U is the Jacobian matrix of the M-DH parameters, and x MDH is the geometric parameter errors of the M-DH. The geometric parameter errors in Eq. (7) can be solved by the inverse kinematic model. However, the ill-posed problem of this model might lead to no solution. For this reason, the coefficient matrix should be modified by a regularization method.

III. REGULARIZATION OF THE GEOMETRIC ERROR MODEL
According to the geometric error model in Eq. (7), the model is built on geometric parameter errors According to the rotation matrix and the translation matrix, the relationship between the frames c and b can be written as where, b c R is the rotation matrix, and b p cO is the translation matrix.
Combining Eq. (7), Eq. (8), and Eq. (9), the geometric error model of the manipulator can be converted into where, p is the reference position, and p is the measured position.
In order to identify the geometric parameter errors in Eq.(10), it is necessary to overcome the problem of no solution or multiple solutions caused by the ill-posedness of the inverse model. Therefore, an appropriate regularization method is designed for the 6DOF manipulator with a highdimensional solution model.

A. TIKHONOV REGULARIZATION
The key to the regularization method is to balance the weighted sum of x MDH and the residual sum of Ux MDH − S and minimize their sum. Mathematical expression where, λ is a regularization parameter to balance the model authenticity and the well-posedness. By differentiating Eq.(11), we can get In details, Then we can solve Eq. (14) as Parameter identification of error model by Eq. (14). On the one hand, the smaller of the regularization parameter, the better of the model authenticity. On the other hand, the larger the value of the λ, the better the stability of the model. Authenticity and stability should be considered at the same time. For this reason, the regularization appropriate parameters should be selected [37], [38].

B. SELECTION OF REGULARIZATION PARAMETER-L CURVE
The equation for the L curve is as follows [39] Ux λ,δ − S δ is the modulus of the residual about λ, x λ,δ is the modulus of the solution about λ. By solving Eq.(15), we can get a curve about u(λ) and v(λ).
In the vertical part, the regularization parameter is a small value, it shows poor stability but good authenticity. In the horizontal part, the regularization parameter is large, it shows good stability but poor authenticity. But, when using the L curve to select the regularization parameter, both the accuracy and the stability can be guaranteed. The curvature is used as a criterion, which expressed as The optimal regularization parameter can be obtained at the inflection point of the L curve. The regularization parameters obtained on the base of L curve are optimized in two steps by stepping method. The detailed flow chart is shown below.

IV. LASER TRACKER CALIBRATION AND ERROR COMPENSATION A. POSITION CALIBRATION
The end position of the manipulator is tested by the laser tracker. Since the laser tracker and the manipulator are not in the same frame, it is essential to make a coordinate conversion by position calibration [40]. The frame between the laser tracker and the manipulator is converted by utilizing the SVD(Singular Value Decomposition) method, which is built as where, q is the weight factor of the measurement data, generally, q i =1, p i is the i th measured position in the laser tracker frame, p i is the i th reference position in the manipulator frame, and w is the number of sampling position.
where, 6 7 T is the transformation matrix of the tool frame relative to the base frame, and 0 7 T is the transformation matrix of the tool frame relative to the manipulator end frame.
The coordinates of the target ball are      x n = p x + n x x + o x y + a x z y n = p y + n y x + o y y + a y z z n = p z + n z x + o z y + a z z (20) According to the LSM method, x, y, and z can be solved by where, D is the distance between the point (x k , y k , z k ) and (x k+1 , y k+1 , z k+1 ) by laser tracker, d is the distance between the point (x m , y m , z m ) and (x m+1 , y m+1 , z m+1 ) by MD-H,

C. ERROR COMPENSATION
The error compensation of the manipulator relies on the geometric error model and regularized parameter identification. When considering the geometric errors, the end pose of the VOLUME 10, 2022  manipulator can be rewritten as

V. NUMERICAL SIMULATION
In order to verify the suggested method, a numerical simulation is carried out. According to the number of identification geometric parameter errors, there are 9 sets of joint angles T = (t i |i = 1, 2, . . . , 6) and their corresponding reference points selected, as Table 2 shows.
Then the coordinate of those 9 points are measured by the laser tracker, as Table 3 shown. Calibrate the measured values to the manipulator frame through Eq.(17)-Eq. (21), and obtain the error between the measured position and the reference position, as shown in Table 3, the geometric error model of the manipulator is obtained by Eq. (4)-Eq. (7). Next, according to the geometric error model, through Eq. (15) and Eq. (16) and the principle of the L curve, the regularization parameters were obtained. Then, the optimization and verification are performed near the obtained regularization parameters to obtain the optimal regularization parameters.
Finally, the identified geometric parameter errors are obtained through Eq. (11), as shown in Table 4. The new geometric parameters are obtained through error compensation, and the nine-points absolute position errors were obtained, as shown in Table 5. The validity of the method is verified according to the position absolute errors.
It can be seen from the simulation results that the absolute position error of the manipulator has been improved to a certain extent after calibration. It indicates that the suggested method has a certain validity.

VI. EXPERIMENT AND ANALYSIS A. EXPERIMENTAL DESIGN
To illustrate the application advantages of the RPIM, the calibration system is designed, which consists of a manipulator, a measurement tool, a target ball, a cooperative target, and a laser tracker, as Fig.3 shows. The selected manipulator has a load capacity of 19 kg, a maximum range of 1024 mm, and a repetitive positioning accuracy of 0.03 mm. The Leica AT403 is used as a measurement tool, and the absolute ranging accuracy of the laser tracker is 5µ m. The target ball is fixed at the end of the manipulator by a connector.
In order to make the experimental results universal, we chose a cube measurement space with a side length of 800mm. The measuring cube is the area where the manipulator moves the most, 52 relatively uniform points are selected in this space. There are 52 scattered points tested to be measured by the laser tracker. To ensure the stability and consistency of the measurement results, each point is tested with 4 seconds pause, as well as at a temperature of around 25 • C The measured 52 points are displayed as the measured points as Fig. 4 shows. It can be seen that there is a difference  between the reference points and the measured points even though they all correspond to the same point.
Obviously, those differences can be offset by a certain operation both the rotation and the translation. The relationship between different frames can be converted by Eq. (17). The end position indicates that using the constraints of Eq.(21), the LSM is fitted to obtain P(-3.402, 15.293, 49.356), the matrix of the tool frame relative to the end frame is After that, unifying the points to the same frame, the points on different frames are calibrated by the SVD conversion to the same frame, as Fig.5 shows. In this way, the differences between the measured points and the reference points on the same scale will really show up.
It can be seen from Fig.5 that the reference points are approximately around the measured points on the same scale, indicating that the calibration results are correct. To illustrate the accuracy of the calibration method, the absolute position errors are solved by those 52 points, as Fig. 6 shows.

B. MATHEMATICAL VALIDATION AND ANALYSIS
Before uploading the calibration errors to the traditional inverse kinematic model to identify the geometric parameters VOLUME 10, 2022  of the manipulator, the ill-posed problem should be judged by its condition number, which can be expressed as where, U is the Jacobian matrix. By using the above data, the condition number is 1.157 × 10 21 , which is a great value to research the inverse kinematic model ill-posed. To improve this situation, the suggested RPIM is used for parameters identification.
To optimize the regularization parameter, the L curve method is described by Eq. (15) and Eq. (16). By using regularization parameter identification, the initial value obtained by the L curve method is 3.063 × 10 −11 . However, the regularization parameters obtained by this method may not be the optimal regularization parameters, so it needs to be further analyzed and optimized. The initial regularization parameter is selected in the range of 3.063 × 10 −20 to 3.063 × 10 −10 , the step frequency is 30. According to the principle of the minimum norm of residual sum, the residual sum (L2) norm are the smallest among the above range when λ = 3.063 × 10 −14 , as Fig.7 shows. It can be seen that the value of the inflection point is 3.06 × 10 −14 . The L2 norm has a trend of decreasing trend with a rapid rate of descent before this point while showing a slowly changing trend after this point.
In order to make the regularization parameters more advantageous, the same order of magnitude analysis is carried out on the basis of the above results, the regularization parameters in the range of 0.563 × 10 −14 to 5.063 × 10 −14 . Ten groups of data were uniformly selected for comparison. All those To verify the accuracy of those simulated points, the absolute position errors are solved for each measured point, as Fig 9 shows. It shows that the absolute position error is minimal when the regularization parameter is 3.063 × 10 −14 compared to the other cases.
As well, the absolute position errors of each point under different regularization parameters are compared. It shows from Fig. 10 that there is a trend of monotone decreasing when the regularization parameter ranges from 0.563×10 −14 to 3.063 × 10 −14 while a trend of monotone increasing when the regularization parameter is greater than 3.063 × 10 −14 . All the error curves have the same tendency, it can be fully demonstrated that the regularization parameter of 3.063 × 10 −14 is superior to other parameters.
According to the regularization parameter 3.063 × 10 −14 , the geometric parameter errors of the inverse kinematic model can be identified, as TABLE 7 shows.
To further validate the accuracy of our suggested method, the LSM, the GN, and the GD were also applied under the same experimental conditions as some comparison methods.
The absolute position errors of the simulated points are compared by NC (no calibration error), LSM, RPIM, GD, and GN. It can be seen from Figure 11 that the errors of the NC are greater than any other method. The average errors are 2.307mm for NC, 1.449mm for LSM, and 0.475mm for RPIM. In addition, the absolute position error of NC is about 6 times higher than that of RPIM, the absolute position error of LSM is about 4.2 times higher than that of RPIM, the absolute position error of GD is about 4 times higher than that of RPIM, and the absolute position error of GN is about 2 times higher than that of RPIM, respectively.
In order to verify the calibration effect of the suggested RPIM method in the actual manipulator model, there are  eighteen points randomly selected. About measured values, all those points are measured by a laser tracker and their corresponding joint angles are measured by the encoder, as TABLE 8 shows.
According to the 18 groups of measure data, the position of the simulated points is solved by RPIM. As a comparison, the corresponding measured points are also plotted on the frame, as Fig. 12 shows.
The RPIM in Fig. 12 represents the position coordinates obtained after the regularization method is calibrated. It can be seen from the figure that the simulated point is very close to the measured point. In detail, the error between each axis is that the x-axis range is -500 to 500, its mean error is 0.671mm, the y-axis range is -200 to 600, its mean error is 1.244mm, and the z-axis range is -600 to 1000, its mean error is 2.815mm. After using RPIM calibration, the x-axis error is reduced by 39.404%, the y-axis error is reduced by 88.489% and the z-axis error is reduced by 99.796%, as Fig.12 shows. In a word, our suggested method has a good effect on manipulator calibration.
It can be seen from Fig.13 that the errors along the three axes are small, the optimization effects on the z-axis are the best one than x-axis and y-axis, and the x-axis are worse than y-axis and z-axis.
As well, the NC, the LSM, the GN, the GD, and the RPIM are compared by the absolute position errors, as Fig.14 shows.
According to the data in Figure 14, the average value can be calculated as 0.468mm for RPIM, 1.435mm for LSM, and 3.199mm for NC, 1.353mm for GD, 0.956mm for GN, the errors of the RPIM, the LSM, the NC, and the GD are decreased by calibration.
According to the 18 groups and 52 groups of the experimental data, the maximum absolute position errors (MAX), standard deviation (SD) and mean absolute position errors (MAE) are also solved to evaluate the localization accuracy of the manipulator, as TABLE 9 shows.
It can be seen from TABLE 9 that the MAX of NC is almost ten times higher than that of RPIM, the MAX of LSM is almost three times higher than that of RPIM, the MAX of GD is almost four times higher than that of RPIM, and the MAX of GN is almost twice times higher than that of RPIM. The MAE of RPIM is almost one-sixth that of NC, one-fourth that of LSM, one-third that of GD, and onehalf that of GN. It shows our suggested method has a high     accuracy than traditional methods. The SD of the RPIM is the smallest one among all methods, which shows our suggested method is stable when it comes to the manipulator calibration.
The reduction percentage both of the MAE and SD relative to the absolute position before calibration is also solved    69.922% by the GD, 85.730% by the GN, and 95.705% by RPIM after calibration. In all, the kinematic model of the manipulator is successfully recognized by RPIM and validated by numerical simulation and experiment.

VII. CONCLUSION
In order to identify the geometric parameter errors of the manipulator more accurately, this paper proposes an RPIM identification method to improve the localization accuracy of the manipulator. On the one hand, when using MD-H to model, it avoids the singularity when two axes are parallel. On the other hand, ill-posed problems are raised during parameter identification, and the regularization method can solve this problem. The RPIM is an improved regularization method based on Tikhonov regularization, and the specific performance is that the improved L curve method is used for selecting regularization parameters.
In addition, the experimental part uses a laser tracker to data test, The RPIM compares with other methods in parameter identification. Experimental results show that manipulator geometric error model based on the Tikhonov regularization can effectively improve the localization accuracy of the pose. In addition, different regularization parameters have different effects on results, resulting in differences in localization accuracy. Comparing with others, RPIM has higher stability and smaller absolute position error, and better improves localization accuracy.

VIII. DECLARATION OF COMPETING INTEREST
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. XIAOFENG HU received the B.S. and M.S. degrees in instrumentation engineering from China Jiliang University, Hangzhou, China. His main research interests include the theory of instrumental accuracy, computer vision, error analysis, precision testing, and auto parts testing. VOLUME 10, 2022