Calibration Method Based on Models and Least-Squares Support Vector Regression Enhancing Robot Position Accuracy

The positioning accuracy of a robot directly affects the quality of its operations. In this study, a calibration method is proposed based on combining a model with least-squares support vector regression (LSSVR) to improve robot positioning accuracy. First, a geometric error model of the robot is established. Second, singular value decomposition (SVD) and physical model analysis method are employed to process the coupling parameters in the error model to improve the accuracy and efficiency of identification. Third, as nongeometric errors hinder the construction of an accurate and complete mathematical model and affect the residual positioning errors of the robot, LSSVR is used to compensate for the residual positioning errors caused by nongeometric errors. The proposed method thus improves the accuracy and robustness of finite sample estimation. Finally, an experiment is performed on an IRB1410 robot with a parallelogram mechanism. The maximum/mean positioning errors of the robot decrease from 2.0348/1.0978 mm to 0.1659/0.0733 mm, and the effectiveness of the proposed method is verified. The proposed method has higher prediction accuracy and stability for small samples than other methods.


I. INTRODUCTION
Robots are typical electromechanical equipment with a wide range of applications in the fields of processing, assembly, and medical treatment. In actual robot application, positioning accuracy is an important performance indicator. With the development of automation and intelligent technology, higher requirements have been proposed for robot positioning accuracy, especially in the fields of pipelines, drilling, and welding. However, processing tolerances, assembly errors, flexible deformation, trajectory planning etc. result in inconsistencies between actual robot models and the nominal model, which deteriorate the robot absolute positioning accuracy [1]- [6]. One solution to this problem is to calibrate The associate editor coordinating the review of this manuscript and approving it for publication was Yangmin Li . the robot kinematics before they are utilized. Kinematics calibration can be model-based and model-free. Model-based calibration is employed to model robot error information, whereas model-free calibration is the approximation of robot kinematics.
Model-based calibration is used to model the robot error source and generally consists of four steps: modeling, measurement, identification, and compensation [7], [8]. An accurate and complete geometric error model should meet the following criteria: completeness, continuity and minimality [9]. The Denavit-Hartenberg (DH) model pro-posed by Denavit et al. is widely employed but becomes discontinuous when adjacent joints are nearly parallel [10], [11]. Researchers have evolved the DH model into the modified DH (MDH) model [10], S model [12], complete and parametrically continuous (CPC) model, and other models [13], [14]. Schröer et al. showed that the MDH model is a subset of the complete, continuous, and minimal kinematics model based on the DH model [9]. The S model violates the principle of continuity, and the CPC model has a problem with parameter redundancy [13]. The product of exponentials (POE) model proposed by He et al. [15] is also characterized by parameter redundancy [16]. Sun et al. pro-posed the more concise finite and instantaneous screw (FIS) theory [17]- [19], which satisfies completeness and continuity. However, the POE and FIS models are abstract descriptions of the transformation of joint motion based on screw theory, which are convenient for calculation and programming but unintuitive and difficult to understand [16]. Measurement equipment, such as laser trackers, are used to obtain pose information for an end-effector [20], [21]. Researchers have employed the least squares method or the Kalman filter to identify model parameters [22]. However, the influence of coupling parameters has prevented the accurate identification of model parameters [23]. Researchers have used model analysis and Jacobian matrix analysis to process coupling parameters [24], [25], but these methods are not suitable for analyzing complex multiparameter coupling models. Stable eigenvectors can be obtained from complex matrices using singular value decomposition (SVD) to determine linear coupling parameters [26], [27]. The nonlinear coupling parameters are analyzed according to the physical model of the robot, and the expressions of the coupling parameters are obtained. The expressions must be verified by a simulation model. Model-based calibration methods focus on analyzing and compensating geometric errors and neglect nongeometric errors. Due to the complexity and diversity of nongeometric errors, model-based calibration methods cannot establish a universal error model that includes all error sources [28].
The modeless calibration method approximates robot kinematics by using a fitting or approximation algorithm. Some researchers have used Fourier polynomials or other polynomials to approximate the robot end-effector positioning error, but this method is limited by complexity and low accuracy [29]. Guo et al. used the joint space grid method to predict nongeometric errors [30]. However, the spatial grid method is complex and only applicable to a small-motion space. Some researchers have used neural networks (NN) to relate the joint angle to the end-effector positioning error to reduce robot positioning error [28], [31], [32]. NNs have higher learning ability and flexibility than other methods and thus have been widely employed in the prediction of nongeometric errors. However, NNs are based on empirical risk minimization and are not suitable for small sample learning, and training results vary from person to person. The support vector machine (SVM) is a machine learning method based on finite samples and the statistical learning theory by Vapnik et al. [33]. The SVM has been widely applied to classification and regression. Unlike empirical risk minimization in NNs, the SVM follows the principle of structural risk minimization, has a strict mathematical theory and foundation, can effectively address finite sample learning, and has a robust estima-tion ability. Least-squares support vector regression (LSSVR) is an extension and improvement of the SVM [34], [35] that improves training speed and solution-finding performance while reducing the calculation complexity. LSSVR can be utilized to predict robot nongeometric error. However, the robot error source cannot be analyzed by model-free calibration or known by the user.
Model-based calibration has many advantages, including accurate representation of error source information, a short calculation time, and fast convergence. However, nongeometric errors cannot be accurately and completely modeled. These nongeometric errors affect the remaining robot positioning errors and reduce the robot positioning accuracy. Model-free calibration can be utilized to compensate the robot residual positioning error. In summary, we can use model-based and model-free calibration to compensate for geometric errors and residual positioning errors after calibration to improve robot positioning accuracy.
In this paper, a calibration method based on combining a model with LSSVR is proposed to compensate robot geometric and nongeometric errors. First, the robot calibration system is introduced, and the MDH model is chosen because of adjacent parallel joints. Second, SVD and error model analysis are employed to process the coupling parameters to improve the accuracy and efficiency of parameter identification. Third, LSSVR is applied to compensate residual positioning errors caused by nongeometric errors and to improve the stability during training with finite samples and the prediction accuracy. Finally, the correctness and effectiveness of the proposed method are verified by performing experiments on the robot calibration system and comparing the results with those of other methods.

II. KINEMATIC CALIBRATION MODEL A. FORWARD KINEMATIC MODEL FOR A ROBOT
The robot calibration system is shown in Figure 1. The system consists of an ABB IRB1410 robot with a parallelogram mechanism, measurement devices, and a spherically mounted retroreflector (SMR).
Within the forward kinematic model of the robot, the transformation matrix from the base coordinate system to the tool coordinate system is expressed as follows: The transformation matrix from the joint {i-1} coordinate system to the joint {i} coordinate system is T i−1 i . The joint angle, joint offset, link length, link twist, and Hayati parameter are denoted as θ i , d i , a i , α i , and β i , respectively.
If two neighboring joints are parallel or approximately parallel, then d i = 0, β i = 0; otherwise, β i = 0, d i = 0. In practice, it is difficult to accurately measure the flange coordinate system. T FL tool is determined as follows: where T FL SMR represents the transformation from the flange coordinate system to the SMR, and T SMR tool represents the fixed transformation from the SMR to the tool coordinate system, which can be measured accurately in the configuration. The transformation matrix from the measurement coordinate system to the SMR is expressed as follows: where T M SMR represents the position of the SMR in the measurement coordinate system, and T M b is the transfor-mation matrix from the robot base coordinate system to the measurement coordinate system. Preliminary measurements of T M b should be carried out before calibration. T FL SMR can be preliminarily determined using a computer-aided design (CAD) model. T M b , T b FL , and T FL SMR are identified simultaneously during calibration. The parameters in the calibration system are listed in Table 1, where T M b is expressed as follows:

B. ERROR MODEL OF THE QUADRILATERAL
Manufacturing and assembly errors cause errors in both the link length and the joint angle, and the parallelogram degenerates into a general quadrilateral, as shown in Figure 2. The quadrilateral mechanism is projected on plane O 2 x 2 y 2 to yield the following closed-chain constraint equation [24]: where cθ i and sθ i represent cos(θ i ) and sin(θ i ), respectively, and θ i,j represents θ i + θ j . It is easy to obtain θ 3 = θ A2 + θ P3 + θ P4 −θ 2 −π from the definition of a quadrilateral. The passive joint θ P3 and the passive joint θ P4 can be linearly related to θ 2 , θ A2 , a 2 , l A2 , l P3 , and l P4 by differentiating Equation (8): The angular error in passive joint 3 can be calculated by Equation (9) as follows: The constraints on the parallelogram l A2 = l P4 and θ 3 = −θ 2 + θ A2 can be used to simplify Equation (9) as follows:

C. ROBOT ERROR MODEL
A robot kinematics model is shown in Equation (5). The relationship between the end-effector positioning error and the geometric parameter error can be expressed as follows [36]: where P is a 3 × 1 end-effector positioning error vector; B is the error vector of the base coordinate system; q is the error vector of the robot parameters; and T is the vector error of the tool coordinate system.
. . , α 6 ] T , and β = [ β 2 ] T . J i is the corresponding identification Jacobian matrix, for which the columns are calculated using the method described in [36].
The Jacobian matrix corresponding to the error parameters of the base coordinate system is expressed as follows: where x i , y i , and z i are the direction vectors of the measurement coordinate system {i}; the coordinate system M is generated after rotating coordinate system {M } around the x M -axis about α 0 ; y M is the direction vector of coordinate system M along the y-axis; and M is generated after rotating M around the y M -axis about β 0 . The vector P b,SMR connects the origin of the base coordinate system to the origin of the end-effector coordinate system.
The Jacobian matrix corresponding to the robot error parameters is expressed as follows: where z i and y i are the direction vectors of the joint coordinate system {i} along the z-axis and y-axis, respectively; the coordinate system i is generated after rotating the joint coordinate system {i − 1} around the z i−1 -axis about θ i ; and x i is the direction vector of the coordinate system i along the x-axis. The vector P i,n connects the origin of the link {i} coordinate system to the origin of the tool coordinates {SMR}. The Jacobian matrix corresponding to the quadrilateral error parameter is calculated by using Equations (9) and (12) as follows: where J I i is the Jacobian matrix corresponding to the quadrilateral error parameter.
The Jacobian matrix corresponding to the error parameters of the flange coordinate system is expressed as follows: (16) where x FL , y FL , and z FL are the direction vectors along the x, y, and z axes, respectively, of the flange coordinate system.

III. DECOUPLING AND IDENTIFICATION OF GEOMETRIC PARAMETERS A. COUPLING PARAMETER ANALYSIS
A complete but redundant calibration model was established in Section II. As redundant parameters can lead to unstable identification, it is essential to remove these parameters to enhance the accuracy and efficiency of geometric parameter identification [23], [37]. SVD is used to determine the linear coupling parameters, and the mathematical expression of the coupling parameters is determined through the physical model analysis method.
By simplifying Equation (12) to P = J g and letting H = J T J , Equation (17) can be calculated by performing SVD on H: where U and V are both unitary matrices that satisfy U T U = I and V T V = I; = diag(σ 1 , σ 2 , . . . , σ r ) represents a diagonal matrix with a nonzero main diagonal; r is the rank of matrix H(r ≤ 36); σ 1 , σ 2 , . . . , σ r are the eigenvalues of H. Thus, there are 36 − r linear coupling parameters.
As H is a symmetrical matrix and V T = U −1 , V is a rotation matrix, and V T g is equivalent to the rotation of g. The coupling parameters of linear correlation are determined by elementary line transformation, and the mathematical expression of the linear correlation coupling parameters is determined through the physical analysis method; the mathematical is as follows: The nonlinear coupling parameters are analyzed using the physical model. The detailed analysis is as follows. The tiny VOLUME 9, 2021  The error θ 6 of θ 6 has the same effect on the end-effector position as the errors X and y t of X and y t and the mathematical expression of θ 6 is obtained according to the trigonometric function relationship among X , y t , X , and y t . The mathematical expressions for the nonlinear coupling parameters a 5 , a 5 , and θ 6 are as follows: where X = a 6 + x t . The coupling relation expression is verified by simulation. The coupling parameters θ 1 , d 1 , l P4 , a 5 , α 5 , θ 6 , x t and z t are removed from the calibration model. Table 2 is a summary of the identified parameters in the calibration model.

B. GEOMETRIC PARAMETER IDENTIFICATION
Measuring the end-effector positions of n joint configurations produces 3 × n equations to identify the robot geometric parameters. The linear error of the system is minimized by using the iterative least squares (ILS) method shown in Figure 3 to identify the model parameters: where g k is the parameter updated upon the kth iteration and g k is the parameter error updated upon the kth iteration, which is expressed by the following equation: where P k = ( P 1 k ) T ( P 2 k ) T · · · ( P n k ) T T is the end-effector positioning error matrix obtained upon the kth iteration, and J k = (J 1 k ) T (J 2 k ) T · · · (J n k ) T T is the corresponding Jacobian matrix. The iteration continues until g k is smaller than ε.

IV. NONGEOMETRIC ERROR PREDICTION
After geometric error calibration, the robot residual error is subject to the influence of the nongeometric error, and LSSVR is employed to compensate for the nongeometric error. The principle of LSSVR is to convert the quadratic programming problem to the solution of a linear equation to enhance the solving speed and accuracy [34], [35]. A training sample (x i , y i ) serves as the input and output vectors of LSSVR. A typical LSSVR model can be described as where ω is the weight vector; φ(x) is a nonlinear function; and b is the curve deviation. The optimization objective of the least-squares model is expressed as follows [38], [39]: where e i is the approximate error of the ith sample, and γ is the regularized parameter used to avoid overfitting. A Lagrange multiplier α i is introduced into Equation (24) to convert the optimization problem to the following equation: where α i is the Lagrange multiplier of each x i . The following equations are derived by setting the partial derivatives of the variables ω, b, e i in Equation (25) to zero: Eliminating e i and ω yields the following equation: where I = [1, 1, . . . , 1] T , α = [α 1 , α 2 , . . . , α n ] T , and is the kernel function for the inner product calculation in the high-dimensional space, and y = (y 1 , y 2 , . . . , y n ) T . Hence, for a new sample x, the output of LSSVR is expressed as follows: The radial basis function (RBF) is a suitable kernel function for nonlinear regression of samples that are large and small, with low and high dimensions. Training reduces the number of support vectors and parameters to be determined, and the RBF kernel function is taken as the kernel function of LSSVR. The regression function can be expressed as follows: where σ is the parameter of the kernel function. The trained LSSVR model is simplified as Lssvr(θ i ), and θ i = [θ 1 , . . . , θ 6 ]. The LSSVR model should be trained before being utilized to compensate for nongeometric errors. Three types of LSSVR are used to predict the residual positioning errors P X , P Y and P Z . The LSSVR input data consist of the joint angle θ i = [θ 1 , . . . , θ 6 ] of the robot and the output data comprising the residual positioning error after geometric error compensation. The LSSVR model is trained and then employed to predict nongeometric errors, as shown in Figure 3.
The entire calibration process of the robot system is shown in Figure 3. The initial parameter g I of the robot is shown in Table 1, and the forward kinematics model for the robot given by Equation (5) is abbreviated as f θ i , g I . Three groups of joint configurations, S 1 , S 2 and S 3 , that cover the entire measurement space are selected, and their corresponding endeffector positions, P S1 , P S2 and P S3 , are measured. S 1 and P S1 are applied to identify the geometric parameters; S 2 and P S2 are applied for LSSVR training, +S 3 and P S3 are used to evaluate and verify the precision after calibration. First, the ILS method is employed to iterate the robot geometric parameter errors until the norm g is less than the termination criterion, whereby the calibrated kinematics model f θ i , g I is obtained. Second, the residual positioning error P 2 of the sample set S 2 is calculated; S 2 and P 2 are utilized as the input and output data, respectively; and the LSSVR model is trained. The trained LSSVR model Lssvr(θ i ) is applied to predict the residual positioning error. After geometric parameter identification and residual positioning error training, the equation of the end-effector is P = f θ i , g I + Lssvr(θ i ), and the verification sets S 3 and P S3 are used to evaluate the robot positioning accuracy.
After calibration, the solution process of the inverse kinematics of the robot is shown in Figure 4. For any desired position P d in the robot's workspace, the iterative inverse solution method is utilized to calculate the joint angle θ mod , and the iteration is terminated when the termination condition e ≤ ε is satisfied. The Jacobian matrix J θ is shown in Equation (14).

V. EXPERIMENTAL VALIDATION
A. EXPERIMENT SYSTEM Figure 5 shows the calibration system consisting of an IRB1410 robot with a parallelogram mechanism, a Leica AT960 laser tracker (positioning accuracy: ±15 µm+6 µm/m), and an SMR. The equipment required for the calibration experiment is shown in Figure 5, but   the laser tracker should be placed outside the workspace of the robot, that is, the distance between the laser tracker and the robot base should be greater than 1.45 m; (2) The SMR installed at the end of the robot should be within the measurable range of the laser tracker, during the experiment, the distance between the laser tracker and the SMR can be set to 2-4 m. During mea-surement, we adjusted the relative pose between the robot and the laser tracker in Figure 5 to meet the requirements of the above two conditions. The repeated positioning precision of the IRB1410 is ±50 µm. Preliminary determination of T M b was conducted using axis fitting [40]. T FL SMR was preliminarily determined using the CAD model. The initial robot parameters are listed in Table 1.
A total of 1200 sets of joint configurations were generated within the maximum measurable space of the robot, and the corresponding end-effector positions were measured. The calibration process is as follows: identification of geometric parameters, training of LSSVR model and evaluation of compensation effect; 100 samples, N samples (N = 50, 100, 150, . . . 1000), and 100 samples covering the entire measuring volume were selected and labeled S 1 , S 2 , and S 3 , respectively; the corresponding end-effector positions were labeled P S1 , P S2 , and P S3 . S 1 was employed for geometric parameter identification, S 2 was utilized for LSSVR training, and S 3 was applied for evaluation and validation of the calibration accuracy.

B. INFLUENCE OF COUPLING ON PARAMETER IDENTIFICATION
The results obtained using Equation (12) for the model with coupling parameters and Equation (20) for the model without coupling parameters were compared in terms of the identification time and identification accuracy.
The ILS method was used to identify the geometric errors based on set S 1 , where ε was 0.0001, and the maximum identification number was 1000. Table 3 shows that the model rank was equal to the number of identified parameters after the coupling parameters had been removed. The identifica-  tion results of the redundant model and the model without redundancy are shown in Table 4 and Table 5, respectively. In particular, the values of θ 0 , θ 1 , d 5 , a 5 , α 5 , d 6 , and z 7 identified by the redundant model showed significant deviation from the nominal parameters, with unreasonable physical interpretations. The calibrated positioning errors were tested in S 3 ; the results indicated that the positioning accuracy was improved after removing the coupling parameters. An analysis of Tables 3, 4, 5, and 6 showed that the efficiency and accuracy of geometric parameter identification was enhanced after removing the coupling parameters, and the proposed decoupling method was validated. The decoupled calibration model was labeled Method 1. Calibration based on Method 1 effectively reduced the influence of geometric errors on the robot positioning accuracy, as shown in Figure 6.

C. NONGEOMETRIC ERROR PREDICTION
The end-effector position in the joint-angle set S 2 was calculated using the calibrated DH parameters and compared with the actual measured position to obtain the calibrated residual positioning errors P x , P y , and P z . The residual positioning error was predicted using three LSSVR models. The LSSVR models were trained using the joint angles in set S 2 as the input and P x , P y , and P z as the output. The trained LSSVR models were then used for nongeometric error prediction.
In this section, the LSSVR model is trained and used to predict the remaining positioning errors. The input to the LSSVR is S 2 , and the outputs are the residual positioning   errors after calibration. We choose BP neural network as the comparison model and labeled it as Method 2. Compared with other nongeometric error compensation methods, BP neural network has a better nonlinear approximation ability, and many studies have used BP neural network for nongeometric error compensation [28], [31], [41].
The number of training samples is the main factor that affects the prediction accuracy of the remaining positioning error. In this section, the influence of the number of training samples on the prediction accuracy of the LSSVR model and NN is investigated by a statistical method. Method 2 was trained five times, and the maximum and average error after each training were recorded. After the training was completed, the mean values of the maximum and average error were calculated. The influence of the sample size on  the prediction accuracy is shown in Figures 7 and 8. In the figures, the average and optimal values for training using Method 2 are plotted and compared with those obtained using the proposed method.
The statistical experimental results shown in Figures 7 and 8 show that for less than 400 training points, the error for Method 2 and the proposed method decreases as the number of samples increases. Using Method 2 for few sample data results in a large fluctuation between the average training error and the optimal training error and weak stability. For more than 400 samples, both the maximum and mean error tend to be stable, and similar optimal training results are obtained using Method 2 and the proposed method. For less than 300 samples, the prediction accuracy of the proposed method is better than that of the optimal training of Method 2. The proposed method has higher accuracy and stability for a small sample than Method 2.
In actual operation, increasing the number of training samples increases the measurement time and workload. Thus, it is necessary to determine an appropriate number of training samples to reduce the workload without reducing the robot accuracy. According to the curves in Figures 7 and 8 obtained from the statistical experimental data, when the number of samples is approximately 400, both the maximum error and the average error tend to be stable; therefore the number of training samples is set at 400. The verification set S 3 is used VOLUME 9, 2021   to evaluate the performance of Method 1, Method 2, and the proposed method. The compensated results are shown in Figure 9 and Table 7, and a part of Figure 9 is enlarged in Figure 10.
The compensated results are shown in Figure 9 and Table 7. Method 1, Method 2 and the proposed method can effectively improve the positioning accuracy of the end-effector. Method 2 and the proposed method can effectively reduce the robot nongeometric error based on the model calibration. The proposed results are almost consistent with the optimal prediction of Method 2. The proposed method can effectively reduce the robot nongeometric error based on the model calibration. The maximum and average error of the robot compensated by the proposed method are 0.1657 mm and 0.0733 mm, respectively, compared with the model-based calibration method, corresponding to an error reduction of 48.4% and 45.38%, respectively. The results indicate that the proposed method can effectively enhance robot positioning accuracy by calibrating the geometric and nongeometric errors. In addition, the proposed method overcomes the shortcomings of Method 2, which varies among persons and has poor stability when training with small samples.

VI. CONCLUSION
To improve robot positioning accuracy, a calibration method that combines a model with LSSVR is proposed to compen-sate the errors in the geometric parameters and the remaining positioning errors of a robot. The main contributions of this paper are listed below.
1. A geometric error model of the robot is established, and the coupling parameters are processed by SVD and a physical model analysis method to improve the identification accuracy and efficiency of the parameters. 2. LSSVR is used to compensate the residual positioning error caused by nongeometric errors. The effectiveness of the proposed method is verified through experiments. Compared with an NN, the proposed method has a higher prediction accuracy for small samples and improves the stability of nongeometric error prediction.