# Kinematic Calibration of Industrial Robots Based on Distance Information Using a Hybrid Identification Method

Authors: Guanbin Gao; Yuan Li; Fei Liu; Shichang Han
Journal: Complexity
 (2021)

---

## Abstract

To improve the positioning accuracy of industrial robots and avoid using the coordinates of the end effector, a novel kinematic calibration method based on the distance information is proposed. The kinematic model of an industrial robot is established. The relationship between the moving distance of the end effector and the kinematic parameters is analyzed. Based on the results of the analysis and the kinematic model of the robot, the error model with displacements as the reference is built, which is linearized for the convenience of the following identification. The singular value decomposition (SVD) is used to eliminate the redundant parameters of the error model. To solve the problem that traditional optimization algorithms are easily affected by data noise in high dimension identification, a novel extended Kalman filter (EKF) and regularized particle filter (RPF) hybrid identification method is presented. EKF is used in the preidentification of the linearized error model. With the preidentification results as the initial parameters, RPF is used to identify the kinematic parameters of the linearized error model. Simulations are carried out to validate the effectiveness of the proposed method, which shows that the method can identify the error of the parameters and after compensation the accuracy of the robot is improved.


---

## Body

## 1. Introduction

With high generality and industrial flexibility, robots have been widely used in modern manufacturing, which play an important role in the fields of automobile manufacturing, logistics, machinery manufacture, and so forth. Traditionally, industrial robots are used in the way of teaching programming, for example, in spot welding, handling, and other similar repetitive works. In recent years, the position accuracy of industrial robots is required to be higher and higher in the fields of mechanical manufacturing, carving, measuring, testing, and so forth.About 90% position errors are caused by the inaccuracy of the kinematic parameters in the controller [1]. The main way to improve the position accuracy of industrial robots is kinematic calibration by which the structural parameters errors due to the tolerance in machining and assembly of robots can be identified and compensated [2]. There are two levels for kinematic calibration. The first level is to identify and compensate the errors between the transducer reading of joints and the actual joint angle. The second level is to identify and compensate all the kinematic parameters of robots. Generally, the kinematic calibration of robots can be classified into four steps [3]: modelling, measurement, parameter identification, and compensation. The kinematic parameter identification of industrial robots is a high nonlinear problem. To obtain precision identified parameters, the kinematic model is required with the features of completeness, continuity, and minimality [4]. The standard Denavit–Hartenberg (DH) model [5] is the most commonly used model in kinematic modelling of serial open-chain robots. However, there will be singularity when two adjacent joints are completely or almost parallel. Hayati [6] proposed a modified DH (MDH) method, which introduces a rotation parameter to solve the singularity problem. S-model was also proposed [7] to solve the singularity problem by Stone. A Complete and Parametrically Continuous (CPC) kinematic model was proposed by [8] to provide a kinematic modelling method with completeness and continuity.In the step of measurement, a group of poses and positions of the robot are acquired for the following parameter identification usually by laser trackers [9–12] or other coordinate measuring instruments [13–15]. Nadia Schillre [16] reported a calibration method with a laser tracker to acquire data. The maximum position error of the robot is reduced from 0.4883 mm to 0.1801 mm. Albert Nubiola [17] presented a calibration case for ABB IRB 1600 using laser tracker, in which the mean position error of the robot was reduced from 0.968 mm to 0.696 mm. With the advantages of high precision, dynamic tracking, and noncontact measuring, laser trackers become the main measurement instrument in the calibration field of robots. But they are so expensive that ordinary users cannot afford them. Except laser trackers, coordinate measuring machines (CMMs) [18], measuring arms [19], ballbars [20], and so forth are also used in calibrating robots. Actually, the coordinates of the end effector of the robots are acquired by instruments in the previous researches [21–23], which are needed in the following identification. These coordinates are acquired in the measurement coordinate system of the instruments, which have to be transformed into the base coordinate system of the robot, and this transformation will bring error and increase the difficulty of identification. Unlike coordinates, distance is a relative value, which is always the same in different coordinate systems without any transformation.The third step, identification, actually is an optimization problem, in which the least square (LS) algorithm [24–27] is the most commonly used way. But LS is easily affected by the noise of data, which can lead to failure in identification. Extended Kalman filter (EKF) algorithm was more suitable for the identification of the robot, the working principle of which is similar to that of Kalman filter [28]. The expected EKF deals with nonlinear systems that retain the Taylor expansion first-order terms of the nonlinear function and ignore the higher-order terms. But EKF cannot deal with high nonlinear system and only suit Gaussian noise system [29]. Actually, the kinematic model of the robot is a highly nonlinear system.Referring to ISO-9258 [30, 31], the distance accuracy of industrial robot can be defined as the deviation between the measured distance and the command distance between the continuous moving points in workspace. To avoid using coordinate measuring instruments and the transformation between the measurement coordinate system and the base coordinate system of the robot, we propose a kinematic calibration method based on distance error modelling and the definition of distance accuracy to improve the position accuracy of the robot. A regularized Particle Filter (RPF) [32] is used to improve the ability of identification for the kinematic model of the robot in the light of the excellent performance of sample filtering method in solving highly nonlinear problem in identification. Since RPF is sensitive to its initial value, EKF is used in preidentification of the kinematic parameters of the robot. The results of the preidentification are used as the initial value of RPF.The paper is organized as follows: In Section2, we introduce the kinematic modelling and distance error modelling. In Section 3, the identification of the robot is discussed and the redundant parameters in the identification model are analyzed by numerical method based on distance error. Comparative simulations are conducted in Section 4, and some conclusions are given in Section 5.


## 2. Kinematic Modelling and Distance Error Modelling

### 2.1. Kinematic Modelling

The controller of industrial robots is essentially a semi-closed-loop system, which sends control instructions to the joint servo motors and acquires the positions of them. The positions of the joints are calculated from the positions of the servo motors, which are used to determine the poses of the industrial robots by kinematic model.An MDH method [33] is used to build the kinematic model of the industrial robot, in which the relationship between two adjacent rotation joints is described by a homogeneous transformation matrix. The kinematic model consists of the controller, which consists of four groups of parameters: link length a, link offset d, link twist angle α, and joint angle θ [33]. First, the coordinate systems of a 6-degree-of-freedom (DOF) industrial robot named ER20-C10 are established according to the MDH method, as shown in Figure 1. Then, the nominal values of kinematic parameters of the robot are obtained, as shown in Table 1.Figure 1
The kinematic model of ER20-C10 robot.Table 1
Nominal values of robot’s MDH model parameters.
No. of jointsαi−1(°)ai−1(mm)θi(°)di (mm)100θ1504290166.605−90 +θ2030−782.27θ304−90−138.826θ4761.355900θ506−900θ6125According to MDH method, the transformation matrix Ti,i−1 of two arbitrary adjacent rotation links {i−1} and {i} is as follows:(1)Ti−1,i=Rotxi−1,αi−1Transai−1,0,0Rotzi,θiTrans0,0,di=cosθi−sinθi0ai−1sinθicosαi−1cosθicosαi−1−sinαi−1−disinαi−1sinθisinαi−1cosθisinαi−1cosαi−1dicosαi−10001.For the generic industrial robot withNt-DOF, the transformation matrixes from the base frame {0} to the end effector can be written as(2)T0,E=∑i=1NtTi−1,i⋅TK,E,where TK,E=00lz1T and lz is the end effector (EE) offset along the axis of joint Nt. The first three rows of the last column of T0,E denote the position of EE, that is, P(·).


### 2.2. Distance Error Modelling

According to ISO-9258 [30, 31], between the continuous moving points during the movement of robot in workspace, the distance accuracy of industrial robot refers to the deviation between the measured distance and the command distance, illustrated as Figure 2. The distance of the two adjacent points (i and i + 1) in workspace can be expanded as(3)Δli,i+1=lti,i+1−lmi,i+1=Pti+1−Pti−Pmi+1−Pmi=Pti+1−Pmi+1−Pti−Pmi.Figure 2
The definition of distance accuracy.The position error of theith point of the acquisition data is defined as(4)Pei=Pti−Pmi=JX,where J is the Jacobean matrix and X is the error of kinematic parameters. According to equations (3) and (4), one can obtain the distance error formula shown as follows:(5)lmi,i+1⋅Δli,i+1=Pmi+1−Pmi⋅Pei+1−Pei.Equation (5) indicates that the distance error can be expressed by partial pose information. Furthermore, equation (5) can be rewritten as [34](6)Δli,i+1=Pmi+1−Pmilmi,i+1⋅Pei+1−Pei.According to equation (6), the measured distance can be calculated by two corresponding adjacent points [34]. Therefore, at every acquisition position, one has to measure twice to obtain the distance used, which is time-consuming and labor-consuming. In this paper, we proposed a modified distance model, which is more convenient and less time-consuming than the general one.The original position of the device is set asP0(xp0,yp0,zp0) in the frame of the measurement device, which is fixed. P0 is considered to be error-free, which is also set as one of the two adjacent points in equation (6). Hence, equation (6) can be rewritten as(7)Δli=Pmi−P0lmi⋅Pei,where lm(i) is the ith distance from EE to P0, which can be directly measured by the measurement system. P0 is known as the fixed point of the measurement system. The original error of the measurement system can be neglectable for the distance is a relative value. With P0 being the fixed reference point, to get the distance, one only needs to measure one time at a position. Thus, the measuring process can be simplified with proposed strategy. Equation (7) describes the relationship between the distance error model and the position error. Assuming that the deviation of the parameters is small, (7) can be linearized as the first-order term of its Taylor expansion around the nominal values of the kinematic parameters, as shown in the following equation:(8)Δli=J˜X,where J˜∈R3n×N, n is the number of calibration distances, and N is the number of parameters. The Jacobian matrix J˜ in distance error model is(9)J˜i=Pmi−P0lmi⋅Ji.To define the problem of identification with filtering algorithm, the errors for kinematic parametersX and the measured data Y are described in the state space, which can be rewritten as follows:(10)X=Δα1⋯Δα6Δa1⋯Δa6Δθ1⋯Δθ6Δd1⋯Δd6T,Xk+1=Xk,Yk=hXk+uk,where k is the number of iterations and Xk represents the parameters’ errors. Yk is the corresponding calibration distance measured by the measurement device, uk is the process noise of measurement, and h (·) is the measurement function.


## 2.1. Kinematic Modelling

The controller of industrial robots is essentially a semi-closed-loop system, which sends control instructions to the joint servo motors and acquires the positions of them. The positions of the joints are calculated from the positions of the servo motors, which are used to determine the poses of the industrial robots by kinematic model.An MDH method [33] is used to build the kinematic model of the industrial robot, in which the relationship between two adjacent rotation joints is described by a homogeneous transformation matrix. The kinematic model consists of the controller, which consists of four groups of parameters: link length a, link offset d, link twist angle α, and joint angle θ [33]. First, the coordinate systems of a 6-degree-of-freedom (DOF) industrial robot named ER20-C10 are established according to the MDH method, as shown in Figure 1. Then, the nominal values of kinematic parameters of the robot are obtained, as shown in Table 1.Figure 1
The kinematic model of ER20-C10 robot.Table 1
Nominal values of robot’s MDH model parameters.
No. of jointsαi−1(°)ai−1(mm)θi(°)di (mm)100θ1504290166.605−90 +θ2030−782.27θ304−90−138.826θ4761.355900θ506−900θ6125According to MDH method, the transformation matrix Ti,i−1 of two arbitrary adjacent rotation links {i−1} and {i} is as follows:(1)Ti−1,i=Rotxi−1,αi−1Transai−1,0,0Rotzi,θiTrans0,0,di=cosθi−sinθi0ai−1sinθicosαi−1cosθicosαi−1−sinαi−1−disinαi−1sinθisinαi−1cosθisinαi−1cosαi−1dicosαi−10001.For the generic industrial robot withNt-DOF, the transformation matrixes from the base frame {0} to the end effector can be written as(2)T0,E=∑i=1NtTi−1,i⋅TK,E,where TK,E=00lz1T and lz is the end effector (EE) offset along the axis of joint Nt. The first three rows of the last column of T0,E denote the position of EE, that is, P(·).


## 2.2. Distance Error Modelling

According to ISO-9258 [30, 31], between the continuous moving points during the movement of robot in workspace, the distance accuracy of industrial robot refers to the deviation between the measured distance and the command distance, illustrated as Figure 2. The distance of the two adjacent points (i and i + 1) in workspace can be expanded as(3)Δli,i+1=lti,i+1−lmi,i+1=Pti+1−Pti−Pmi+1−Pmi=Pti+1−Pmi+1−Pti−Pmi.Figure 2
The definition of distance accuracy.The position error of theith point of the acquisition data is defined as(4)Pei=Pti−Pmi=JX,where J is the Jacobean matrix and X is the error of kinematic parameters. According to equations (3) and (4), one can obtain the distance error formula shown as follows:(5)lmi,i+1⋅Δli,i+1=Pmi+1−Pmi⋅Pei+1−Pei.Equation (5) indicates that the distance error can be expressed by partial pose information. Furthermore, equation (5) can be rewritten as [34](6)Δli,i+1=Pmi+1−Pmilmi,i+1⋅Pei+1−Pei.According to equation (6), the measured distance can be calculated by two corresponding adjacent points [34]. Therefore, at every acquisition position, one has to measure twice to obtain the distance used, which is time-consuming and labor-consuming. In this paper, we proposed a modified distance model, which is more convenient and less time-consuming than the general one.The original position of the device is set asP0(xp0,yp0,zp0) in the frame of the measurement device, which is fixed. P0 is considered to be error-free, which is also set as one of the two adjacent points in equation (6). Hence, equation (6) can be rewritten as(7)Δli=Pmi−P0lmi⋅Pei,where lm(i) is the ith distance from EE to P0, which can be directly measured by the measurement system. P0 is known as the fixed point of the measurement system. The original error of the measurement system can be neglectable for the distance is a relative value. With P0 being the fixed reference point, to get the distance, one only needs to measure one time at a position. Thus, the measuring process can be simplified with proposed strategy. Equation (7) describes the relationship between the distance error model and the position error. Assuming that the deviation of the parameters is small, (7) can be linearized as the first-order term of its Taylor expansion around the nominal values of the kinematic parameters, as shown in the following equation:(8)Δli=J˜X,where J˜∈R3n×N, n is the number of calibration distances, and N is the number of parameters. The Jacobian matrix J˜ in distance error model is(9)J˜i=Pmi−P0lmi⋅Ji.To define the problem of identification with filtering algorithm, the errors for kinematic parametersX and the measured data Y are described in the state space, which can be rewritten as follows:(10)X=Δα1⋯Δα6Δa1⋯Δa6Δθ1⋯Δθ6Δd1⋯Δd6T,Xk+1=Xk,Yk=hXk+uk,where k is the number of iterations and Xk represents the parameters’ errors. Yk is the corresponding calibration distance measured by the measurement device, uk is the process noise of measurement, and h (·) is the measurement function.


## 3. Kinematic Parameter Identification

The effect of gradual error accumulation of links decreases the absolute position accuracy significantly for the open-chain structure of industrial robots. Kinematic parameter identification is an important way to improve the absolute position accuracy.In the kinematic parameter identification, the data of the joint angles and the corresponding distances are required. The forward kinematic solutions of the robot can be obtained by combining a group of joint angles, with which the distances from EE of the robot toP0 can be calculated. Unlike the kinematic parameter identification using laser tracker to acquire the coordinates of EE, the proposed method in this paper only needs the distances acquired by pull wire sensors or another distance measurement system.


### 3.1. Preidentification Based on EKF

The solution for function (8) can be obtained commonly using a nonlinear LS approach. Although nonlinear LS approach is a fast and efficient algorithm that can be used to solve nonlinear equations, it is very sensitive to noise. To obtain more accurate and stable solutions, EKF algorithm is used in preliminary identification.The EKF algorithm consists of two steps: prediction and update [29]. The parameters’ error X can be regarded as the state vector to be preliminarily identified. First, the prior estimate of state vector X^k+1 and the covariance matrix are predicted by the former state vector Xk and Pk.(11)X^k+1=Xk,P^k+1=Pk+Qk,where k denotes the number of iterations, Xk is N-state vector consisting of parameters’ deviation at kth iteration, X^k+1 is the prior estimation of state vector, Pk is the covariance matrix of the estimation state, and Qk is an N×N covariance matrix of the system process noise.With the prior estimate of state vector, the optimal Kalman gain can be calculated with (12), and the covariance matrix of state vector is obtained, as shown in (14). Then, the state vector can be updated by the measurement data with (13).(12)Kk+1=P^k+1J˜k+1TJ˜k+1P^k+1J˜k+1T+Rk+1−1,(13)Xk+1=X^k+1+KYk+1−J˜k+1TX^k+1,(14)Pk+1=P^k+1−Kk+1J˜k+1P^k+1,where Kk + 1 is the Kalman gain at kth iteration and Rk + 1 denotes the covariance matrix of measurement noise. Yk + 1 is the distance error measured by the measurement system. J˜ is the Jacobian matrix. Xk + 1 is the posterior estimate of state vector. The details of the EKF algorithm are illustrated in Figure 3.Figure 3
The flow chart of the EKF algorithm.


### 3.2. The RPF Algorithm

When the noise of the acquired data is non-Gaussian, it is difficult for the EKF algorithm to get a reliable identified result. The RPF algorithm does not need the model to be linear or compellent hypothesis for Gaussian noise [35], the basis of which is particle filter. The preidentification result of EKF is applied as the initial value of the RPF which can identify the state of nonlinear and non-Gaussian noise system. With a certain amount of random state samples (i.e., particles), the RPF algorithm can approximate the posterior probability density function of the identifying model. It provides a suboptimal solution with a finite number of samples. All state samples are transferred in dynamic systems through importance sampling. Subsequently, it updates the posterior distribution sequentially.Apart from the resampling step, the RPF is similar to generic sampling importance resampling (SIR) filter [36]. Considering the robot that is described by the state space model, the RPF algorithm operates Np particles to approximate the posterior density.The RPF algorithm also consists of two steps, prediction and update, similar to EKF. But it is described statistically with the state transition probability densitypxk|xk−1 and the observation likelihood probability density of the state vector pxk|y1:k.In the prediction step, the predicted state vector is estimated by the last time stateXk-1 and its measurement data Y1:k − 1.(15)pXk|Y1:k−1=∫pXk|Yk−1pXk−1|Y1:k−1dXk−1.With the latest measurement dataY1:k, the predicted state can be updated; then, its approximated discrete posterior density can be obtained:(16)pXk|Y1:k=pYk|XkpXk|Y1:k−1pYk|Y1:k−1,where pYk|Y1:k−1 is defined as the normalization constant, as shown in (20).(17)pYk|Y1:k−1=∫pYk|XkpXk|Y1:k−1dXk.With Monte Carlo method [37], the approximated discrete posterior density can be simplified, as shown in the following equation:(18)pXk|Y1:k≈∑i=1NpωkiδXk−Xki,where δ· is the Dirac delta function, Np is the number of particles, and ωki is the weight of the ith particle at kth iteration. pXk|Y1:k is the approximated discrete posterior density. According to Bayesian estimation [37], all the particles need to sample from the target posterior density with the knowledge of observation Y1:k, but this is unrealistic. Importance sampling [37] is introduced in the RPF algorithm, as well as the importance samples from a known posterior density (i.e., importance density). In the identification, the weight of the particles and the normalized weight can be defined as (19) and (20), respectively.(19)ωki=12πRexp−12ekdTR−1ekd,(20)ω˜ki=ωki∑i=1Npωki,where ekd is the distance error of the ith particle. According to Bayesian estimation [37], the minimum mean square (MMS) estimate X^ki can be described as the expectation of identified state vector combining all particles [38], as shown in the following equation:(21)X^ki≈EX^ki≈∑i=1Npω˜kiXki.In the identification model with distance error, the weight of each particle is evaluated using its distance error. The particles with lower distance error Δl have higher weight; that is,(22)Δl=lm−lr.After a number of iterations, the weight will occupy on one particle (almost close to one) and the weights of other particles are so small that they will be negligible due to the degeneracy of particles [37]. The degeneracy problem indicates that abundant computation is applied to updating particles whose weight is almost zero. The negligible particles with tiny weight make no contribution to approximation of the target posterior density. Resampling [39] should be applied to avoid this problem. Note that an effective sample size Neff is introduced to measure the degeneracy of particles [39]. If Neff is too small, resample will be triggered.(23)Neff=1∑i=1Npωki2.Resampling is an approach to prevent particles from degeneracy, which will introduce the decay of diversity of particles in turn. This problem leads to poor approximation of the posterior density. One way to solve this problem is modifying the resampling, in which the particles are resampled from the continuous approximation of the posterior density [40].The approximated continuous posterior density in the RPF [40] is(24)pX0:k|Y1:k≈∑i=1NpωkiKhX0,k−X0,ki,where(25)khx=1hnxkxh,where kh(·) is the Kernel probability density function, h > 0 is the Kernel bandwidth, and n is the dimension of the estimated state vector X. kh will be used for resampling particles on continuous approximated density. With appropriate function kh and h, the mean integrated square error (MISE) between real and regularized empirical posterior densities in (24) can be minimized [38]. The details of the RPF algorithm are illustrated in Figure 4.Figure 4
The flow chart of the RPF algorithm.It is assumed that the coordinate ofP0 is constant, which can be regarded as one of identified parameters. The offset of the adapter of the measurement device should also be considered. In general, a complete kinematic model includes 28 unknown kinematic parameters, which need to be identified in our framework. The distances used in identification should be more than 10 mm to increase the certainty of identified results. Figure 5 shows the flow chart of the identification.Figure 5
The flow chart of the proposed algorithm.


### 3.3. Determination of Redundant Parameters

Singularity analysis of the identified matrixJ is very important in identification calculation, which can lead to nonconvergence or inaccurate results. Generally speaking, the maximum number of parameters is 4r1 + 2p1 + 6 in a generic serial-link robot with r1 rotation axes and p1 translation axes [41]. To obtain reliable results with the identification algorithm, redundant parameters must be determined and excluded from the error model before identification. The condition number of the identified matrix is(26)κJ=JJ−1,where κ⋅ is the measure of sensitivity or stability of the identified matrix with regard to a small change in the input argument. In high dimension full-parameter identification, with a larger condition number, the identified matrix J is often nonfull rank. In this paper, the singular value decomposition (SVD) [42] is applied to reduce the condition number and avoid singularity problem. Considering total Jacobian matrix H = [J1, J2, …, Jn], it can be written as(27)H=U∑VT,where U∈R3n×3n, V∈R25×25, and they are orthogonal matrices; H∈R3n×N, Σ=diagσ1,σ2,σ3,…,σr, r is the rank of H, n is the number of calibration distances, N is the number of parameters, and, in this case, N is 25. Hence, the number of the parameters that are redundant can be obtained as 25 − r.After SVD calculation,r can be determined, that is, r = 20. Therefore, there are 5 redundant parameters that need to be eliminated from the process of identification. According to the property of distance error, the original coordinates of P0 are determined in the base frame. From Figure 1, it can be found that the four parameter errors in joint 1 may cause the overall change of the global coordinate system. Hence, P0 will change accordingly, which makes it identifiable. Therefore, there are 19 parameters including P0 (xp0, yp0, zp0) which need to be identified; and they can be written in the format of a vector:(28)X=Δα1⋯Δα5Δa1⋯Δa4Δθ2⋯Δθ5Δd2Δd4Δd6T.


## 3.1. Preidentification Based on EKF

The solution for function (8) can be obtained commonly using a nonlinear LS approach. Although nonlinear LS approach is a fast and efficient algorithm that can be used to solve nonlinear equations, it is very sensitive to noise. To obtain more accurate and stable solutions, EKF algorithm is used in preliminary identification.The EKF algorithm consists of two steps: prediction and update [29]. The parameters’ error X can be regarded as the state vector to be preliminarily identified. First, the prior estimate of state vector X^k+1 and the covariance matrix are predicted by the former state vector Xk and Pk.(11)X^k+1=Xk,P^k+1=Pk+Qk,where k denotes the number of iterations, Xk is N-state vector consisting of parameters’ deviation at kth iteration, X^k+1 is the prior estimation of state vector, Pk is the covariance matrix of the estimation state, and Qk is an N×N covariance matrix of the system process noise.With the prior estimate of state vector, the optimal Kalman gain can be calculated with (12), and the covariance matrix of state vector is obtained, as shown in (14). Then, the state vector can be updated by the measurement data with (13).(12)Kk+1=P^k+1J˜k+1TJ˜k+1P^k+1J˜k+1T+Rk+1−1,(13)Xk+1=X^k+1+KYk+1−J˜k+1TX^k+1,(14)Pk+1=P^k+1−Kk+1J˜k+1P^k+1,where Kk + 1 is the Kalman gain at kth iteration and Rk + 1 denotes the covariance matrix of measurement noise. Yk + 1 is the distance error measured by the measurement system. J˜ is the Jacobian matrix. Xk + 1 is the posterior estimate of state vector. The details of the EKF algorithm are illustrated in Figure 3.Figure 3
The flow chart of the EKF algorithm.


## 3.2. The RPF Algorithm

When the noise of the acquired data is non-Gaussian, it is difficult for the EKF algorithm to get a reliable identified result. The RPF algorithm does not need the model to be linear or compellent hypothesis for Gaussian noise [35], the basis of which is particle filter. The preidentification result of EKF is applied as the initial value of the RPF which can identify the state of nonlinear and non-Gaussian noise system. With a certain amount of random state samples (i.e., particles), the RPF algorithm can approximate the posterior probability density function of the identifying model. It provides a suboptimal solution with a finite number of samples. All state samples are transferred in dynamic systems through importance sampling. Subsequently, it updates the posterior distribution sequentially.Apart from the resampling step, the RPF is similar to generic sampling importance resampling (SIR) filter [36]. Considering the robot that is described by the state space model, the RPF algorithm operates Np particles to approximate the posterior density.The RPF algorithm also consists of two steps, prediction and update, similar to EKF. But it is described statistically with the state transition probability densitypxk|xk−1 and the observation likelihood probability density of the state vector pxk|y1:k.In the prediction step, the predicted state vector is estimated by the last time stateXk-1 and its measurement data Y1:k − 1.(15)pXk|Y1:k−1=∫pXk|Yk−1pXk−1|Y1:k−1dXk−1.With the latest measurement dataY1:k, the predicted state can be updated; then, its approximated discrete posterior density can be obtained:(16)pXk|Y1:k=pYk|XkpXk|Y1:k−1pYk|Y1:k−1,where pYk|Y1:k−1 is defined as the normalization constant, as shown in (20).(17)pYk|Y1:k−1=∫pYk|XkpXk|Y1:k−1dXk.With Monte Carlo method [37], the approximated discrete posterior density can be simplified, as shown in the following equation:(18)pXk|Y1:k≈∑i=1NpωkiδXk−Xki,where δ· is the Dirac delta function, Np is the number of particles, and ωki is the weight of the ith particle at kth iteration. pXk|Y1:k is the approximated discrete posterior density. According to Bayesian estimation [37], all the particles need to sample from the target posterior density with the knowledge of observation Y1:k, but this is unrealistic. Importance sampling [37] is introduced in the RPF algorithm, as well as the importance samples from a known posterior density (i.e., importance density). In the identification, the weight of the particles and the normalized weight can be defined as (19) and (20), respectively.(19)ωki=12πRexp−12ekdTR−1ekd,(20)ω˜ki=ωki∑i=1Npωki,where ekd is the distance error of the ith particle. According to Bayesian estimation [37], the minimum mean square (MMS) estimate X^ki can be described as the expectation of identified state vector combining all particles [38], as shown in the following equation:(21)X^ki≈EX^ki≈∑i=1Npω˜kiXki.In the identification model with distance error, the weight of each particle is evaluated using its distance error. The particles with lower distance error Δl have higher weight; that is,(22)Δl=lm−lr.After a number of iterations, the weight will occupy on one particle (almost close to one) and the weights of other particles are so small that they will be negligible due to the degeneracy of particles [37]. The degeneracy problem indicates that abundant computation is applied to updating particles whose weight is almost zero. The negligible particles with tiny weight make no contribution to approximation of the target posterior density. Resampling [39] should be applied to avoid this problem. Note that an effective sample size Neff is introduced to measure the degeneracy of particles [39]. If Neff is too small, resample will be triggered.(23)Neff=1∑i=1Npωki2.Resampling is an approach to prevent particles from degeneracy, which will introduce the decay of diversity of particles in turn. This problem leads to poor approximation of the posterior density. One way to solve this problem is modifying the resampling, in which the particles are resampled from the continuous approximation of the posterior density [40].The approximated continuous posterior density in the RPF [40] is(24)pX0:k|Y1:k≈∑i=1NpωkiKhX0,k−X0,ki,where(25)khx=1hnxkxh,where kh(·) is the Kernel probability density function, h > 0 is the Kernel bandwidth, and n is the dimension of the estimated state vector X. kh will be used for resampling particles on continuous approximated density. With appropriate function kh and h, the mean integrated square error (MISE) between real and regularized empirical posterior densities in (24) can be minimized [38]. The details of the RPF algorithm are illustrated in Figure 4.Figure 4
The flow chart of the RPF algorithm.It is assumed that the coordinate ofP0 is constant, which can be regarded as one of identified parameters. The offset of the adapter of the measurement device should also be considered. In general, a complete kinematic model includes 28 unknown kinematic parameters, which need to be identified in our framework. The distances used in identification should be more than 10 mm to increase the certainty of identified results. Figure 5 shows the flow chart of the identification.Figure 5
The flow chart of the proposed algorithm.


## 3.3. Determination of Redundant Parameters

Singularity analysis of the identified matrixJ is very important in identification calculation, which can lead to nonconvergence or inaccurate results. Generally speaking, the maximum number of parameters is 4r1 + 2p1 + 6 in a generic serial-link robot with r1 rotation axes and p1 translation axes [41]. To obtain reliable results with the identification algorithm, redundant parameters must be determined and excluded from the error model before identification. The condition number of the identified matrix is(26)κJ=JJ−1,where κ⋅ is the measure of sensitivity or stability of the identified matrix with regard to a small change in the input argument. In high dimension full-parameter identification, with a larger condition number, the identified matrix J is often nonfull rank. In this paper, the singular value decomposition (SVD) [42] is applied to reduce the condition number and avoid singularity problem. Considering total Jacobian matrix H = [J1, J2, …, Jn], it can be written as(27)H=U∑VT,where U∈R3n×3n, V∈R25×25, and they are orthogonal matrices; H∈R3n×N, Σ=diagσ1,σ2,σ3,…,σr, r is the rank of H, n is the number of calibration distances, N is the number of parameters, and, in this case, N is 25. Hence, the number of the parameters that are redundant can be obtained as 25 − r.After SVD calculation,r can be determined, that is, r = 20. Therefore, there are 5 redundant parameters that need to be eliminated from the process of identification. According to the property of distance error, the original coordinates of P0 are determined in the base frame. From Figure 1, it can be found that the four parameter errors in joint 1 may cause the overall change of the global coordinate system. Hence, P0 will change accordingly, which makes it identifiable. Therefore, there are 19 parameters including P0 (xp0, yp0, zp0) which need to be identified; and they can be written in the format of a vector:(28)X=Δα1⋯Δα5Δa1⋯Δa4Δθ2⋯Δθ5Δd2Δd4Δd6T.


## 4. Verification

To verify the proposed kinematic parameter identification method, a 6-DOF industrial robot ER20-C10 was used, whose maximum load and repeatability are 20 kg and 0.08 mm, respectively. The nominal values of parameters are shown in Table1. The settings of the error for each parameter are shown in Table 2. It should be noted that the error of the redundant parameters is set as zero, since these parameters are not independent. The two termination conditions of the iteration are as follows: the maximum number of iterations is 10000 and the minimum objective function is less than 10−6. Both of the covariance matrices Q and P were initialized as 10−4 × I3∗3 in EKF; and the noise of the covariance matrices of measurement R is set as 10−4 × I3∗3 in the simulation. The bandwidth h is set as 0.42. A random and bounded measurement disturbance [−0.02 mm, 0.02 mm] is added to each distance to simulate the measuring error in the actual device.Table 2
Settings of the kinematic parameters’ error.
No. of jointsΔαi−1 (°)Δai−1 (mm)Δθi(°)Δdi (mm)1000020.021.00.03−0.33−0.010.50.04040.015−0.40.060.855−0.01−0.30.0506−0.0600−0.4To verify the effectiveness of the proposed algorithm, a group of 100 poses were generated randomly for identification and verification. Considering the workspace of the robot, the six joint angles were chosen randomly from [−30°,90°], [−50°,50°], [−45°,90°], [−100°,100°], [−90°,90°], and [−360°,360°], respectively. The former 50 poses were used for identification and the other 50 poses were used for validation.To compare the effectiveness, three different algorithms (i.e., LS, EKF, and EKF + RPF) were used to identify the kinematic parameters that were compensated to the control model of the robot. The maximum distance error was used as the index to evaluate the results of compensation. Figure6 is the distance error distribution chart of the robot after compensation by LS, which shows that all the distance errors are reduced greatly, and the average error is from 0.4827 mm to 0.1120 mm by a reduction rate of 76%. Figure 7 is the distance error distribution chart of the robot after compensation by EKF, which shows that the average error is reduced from 0.4827 mm to 0.1352 mm by a reduction rate of 72%. Figure 8 is the distance error distribution chart of the robot after compensation by EKF + RPF, which shows that the average error is reduced from 0.4827 mm to 0.0780 mm by a reduction rate of 84%.Figure 6
The distance error of the robot after compensation by LS.Figure 7
The distance error of the robot after compensation by EKF.Figure 8
The distance error of the robot after compensation by EKF + RPF.Table3 shows the statistics of the maximum, mean, and RMS values of the distance errors before and after compensation by the three algorithms. By comparison, one can find that the effectiveness of RPF is the best of the three.Table 3
The distance error of robot before and after compensation by 3 different algorithms.
AlgorithmMax (mm)Mean (mm)RMS (mm)Before compensation—1.26120.48270.5791After compensationLS0.32360.11200.1376EKF0.32550.13520.1653EKF + RPF0.24040.07800.0936


## 5. Conclusion

This paper proposed a novel kinematic calibration method based on the distance information for the industrial robot, in which an EKF and RPF hybrid identification method was adopted. To simplify the data acquisition procedure and also avoid using expensive coordinate measuring instruments, the moving displacements of the end effector were used as the reference to build the kinematic parameter error identification model. Then, the error model was linearized by differentiating of the coefficient matrix. By SVD of the coefficient matrix, the redundant parameters were determined and removed from the list of parameters to be identified. The preidentification of the parameters was conducted with EKF, by which we obtained a group of initial parameters of the kinematic model. RPF was used to identify the kinematic parameters of the linearized error model with the initial parameters. Comparative simulations were conducted based on the linearized kinematic error model. The compensation results of the three algorithms are all effective in improving the accuracy of the industrial robot, but the hybrid algorithm proposed in this paper is more precise and fast in iteration calculation. In the future work, precision distance measuring device will be developed, and effectiveness of the device and the proposed algorithm will be further verified by experimental tests.



---

*Source: 8874226-2021-03-26.xml*