# Improvement of Heavy Load Robot Positioning Accuracy by Combining a Model-Based Identification for Geometric Parameters and an Optimized Neural Network for the Compensation of Nongeometric Errors

Authors: Yuxiang Wang; Zhangwei Chen; Hongfei Zu; Xiang Zhang; Chentao Mao; Zhirong Wang
Journal: Complexity
 (2020)

---

## Abstract

The positioning accuracy of a robot is of great significance in advanced robotic manufacturing systems. This paper proposes a novel calibration method for improving robot positioning accuracy. First of all, geometric parameters are identified on the basis of the product of exponentials (POE) formula. The errors of the reduction ratio and the coupling ratio are identified at the same time. Then, joint stiffness identification is carried out by adding a load to the end-effector. Finally, residual errors caused by nongeometric parameters are compensated by a multilayer perceptron neural network (MLPNN) based on beetle swarm optimization algorithm. The calibration is implemented on a SIASUN SR210D robot manipulator. Results show that the proposed method possesses better performance in terms of faster convergence and higher precision.


---

## Body

## 1. Introduction

The positioning accuracy of robots is of great significance in advanced robotic manufacturing systems. In general, the robot shows good repeatability, but poor positioning accuracy. Many factors would result in errors in positioning, including set-up errors, manufacturing tolerance, wear and tear, transmission errors, and compliance [1]. Therefore, researchers have focused on how to improve the positioning accuracy of a robot by errors compensation according to the factors mentioned above.The calibration technique can improve the robot’s positioning accuracy through software algorithm without changing its mechanical structure or design. Generally, calibration can be classified into two types, namely geometric and nongeometric calibration. Almost 80% of robot positioning errors are caused by errors of geometric parameters, including link length errors, link twist angle errors, link offset, and joint angle offsets [2]. Therefore, most researches have concentrated on kinematic-based calibration [3–9]. Most of the researchers employed the D-H model which employs a minimum set of parameters to describe the relationship of joint coordinates for calibration. However, D-H model is not continuous when two consecutive joint axes are parallel or nearly parallel. In order to overcome these drawbacks, Hayati proposed a new model namely MDH model by adding a redundant parameter to the D-H model [10]. However, both the two methods are complicated in the process of modeling a manipulator. Brockett firstly used the POE formula to describe an open-chain robot [11]. The POE formula possesses a simpler expression to describe the relationship between the joint angles and the end-effector’s position and posture. Additionally, previous studies mostly focused on the calibration of the link lengths and the joint offsets, and little attention has been paid to the errors of the reduction ratio and the coupling ratio. The reduction ratio and coupling ratio errors mainly caused by the deviation during manufacturing and assembly process will lead to joint angle errors, further resulting in a decline in the positioning accuracy. Moreover, owing to the nonlinear trigonometric relationship between joint angles and the pose of endpoint, it is complicated to identify the transmission errors using the DH model [12]. The POE formula, however, possesses a simpler expression.However, the kinematic-based calibration has the following limitations: (1) the model only takes the geometric factors into account, but the nongeometric factors such as backlash and joint compliance, which may have significant influence on positioning accuracy, are ignored. (2) The identification process is a complex numerical procedure which may suffer from the numerical problem of ill-conditioning. (3) The implementation of the identified model is problematic due to the difficulty of modifying the controller parameters [13].In fact, the nongeometric errors are an important component of the parameters identification, especially in high-precision application scenarios. The positioning errors caused by joint stiffness account for 6%∼8% as well as other nongeometric errors such as gear backlash and temperature drift [14]. Some researchers investigated both the kinematic parameters and joint compliance [15, 16]. However, other nongeometric errors still affect the positioning accuracy significantly. Especially for a heavy load robot, not only the joint compliance but also the link compliance affects the accuracy [17, 18].Recently, for nonparametric calibration, some intelligent algorithms, such as genetic algorithm [7, 19, 20], maximum likelihood estimation [21, 22], neural network [23–29], and various hybrid algorithms [30], have emerged. Among them, neural network has been widely employed to build the relationship, especially nonlinear relationship between inputs and outputs. It can approximate any nonlinear function with arbitrary precision [31]. The neural network is utilized to build the relationship between the end-effector position and the error of the position in [25]. However, the positioning error depends on the configurations. In [26], a multilayer perceptron neural network is utilized to describe the relationship between the joint angles and the corresponding joint errors. Nguyen proposed a technique for the calibration of industrial robots by combining the geometric model-based calibration method and the ANN to identify the kinematic errors, joint compliance errors, and the nongeometric errors [29]. But, at the same time, the BP neural network has the problem that the performance relies too much on the input data and initial values of weights and biases and easily falls into the local optimum. To overcome these drawbacks, researchers have proposed many hybrid algorithms such as GA-BPNN [32, 33] and PSO-BPNN [34]. But, convergence rate of these algorithms is relatively slow, and implementation of code is complex.In this paper, an enhanced POE-based method is proposed for the identification of geometric parameters, reduction ratio, and coupling ratio. Joint stiffness identification is carried out by adding a load to the end-effector. For the nongeometric errors compensation, an optimized MLPNN based on a hybrid optimization method of beetle antennae search algorithm and particle swarm optimization is proposed. Experiments are carried out on a SIASUN SR210D robot manipulator to validate the proposed method by using a laser tracker to measure the position of the Spherical Mounted Retroreflector (SMR) located at the end flange of the manipulator.The rest of the paper is organized as follows. In Section2, the kinematic model of the SIASUN SR210D robot manipulator is presented. In Section 3, the POE-based error model including the transmission errors is established. Section 4 proposes the joint stiffness identification method. In Section 5, a hybrid beetle swarm optimization and MLPNN algorithm are presented for the nongeometric errors compensation. In Section 6, Compensation verification and analysis are presented. In Section 7, conclusions and future works are presented.


## 2. Kinematic Model of the SIASUN SR210D Robot Manipulator

In this section, a kinematic model of the SIASUN SR210D robot manipulator is described in detail. The structure of the 6R manipulator is in Figure1. L1,L2,L3,L4,andL5 are the link lengths and ω1,ω2,ω3,ω4,ω5,andω6 are the joint angular velocities, respectively. The nominal value of robot structure parameters mentioned above is listed in Table 1.Figure 1
Structure of 6R SIASUN SR210D manipulator.Table 1
Nominal value of robot structure parameters.
L
i, i = 1–5 (mm)
312
1075
235
1282
260
dθi, i = 1–6 (°)
0
0
0
0
0
0Then, the screw of each joint is(1)ξ1=000001ξ2=00−L10−10ξ3=L20−L10−10ξ4=0L2+L30100ξ5=L2+L30−L1−L40−10ξ6=0L2+L40100.The forward kinematics can be expressed as mentioned below:(2)PcT=eξ^1θ1⋯eξ^nθnPcT0,where θi is the angle of the ith joint and PcT0 is the end-effector position in the base frame corresponding to θi=0.


## 3. Geometric Parameters and Transmission Errors Identification

### 3.1. Model of Geometric Parameters Errors

The position of the tool center point (TCP) in the measure frame is shown below:(3)P=M⋅PcT=M⋅gθ⋅PcT0,where P is the position of the TCP, M is a homogenous matrix that represents the transformation from the measuring device frame to the base frame, and PCT represents the position of TCP expressed in the base frame corresponding to θ=θ1,θ2,…,θn, gθ=eξ^1ϑ1eξ^2ϑ2⋯eξ^nϑn.Then, we have the deviation of the TCP position by taking the differential on both sides of equation (3):(4)dP=dM⋅g⋅Pcr0+M⋅dg⋅Pcr0+M⋅g⋅dPcr0,=dM⋅M−1⋅P+M⋅dg⋅g−1⋅Pcr+M⋅g⋅dPcr0,where(5)g⋅dPcT0=Rgbg00dPcT00=RgdPcT00.Taking a close observation of (4), we notice that dM × M−1 ∈ se(3), dg × g−1 ∈ se(3), and the element s belonging to lie algebra satisfies(6)sq∗=ω^Vo0q1=ω^q+υ0=I3−q^00υω.d
P can be expressed as(7)dP=I−P^00dM⋅M−1V+MI−P^cT00dg⋅g−1V+MRgdPcT00.Since rigid body motion can be realized by a rotation and a translation about an axis,M can be expressed as(8)M=eξ^ϑ0∈SE3.Take the differential on both sides:(9)dM⋅M−1=deξ^0ϑ0⋅e−ξ^0ϑ0=θ0∫01eξ^0ϑ0sdξ^0e−ξ^0ϑ0sds.The matrixg is(10)g=Rgbg01∈SE3.Take the differential on both sides:(11)dg⋅g-1=deξ^1θ1e-ξ^1θ1+eξ^1θ1deξ^2θ2e-ξ^2θ2e-ξ^1θ1+⋯+eξ^1θ1eξ^2θ2…eξ^n-1θn-1deξ^nθne-ξ^nθn…e-ξ^2θ2e-ξ^1θ1,=deξ^1θ1e-ξ^1θ1+Adeξ^1θ1deξ^2θ2e-ξ^2θ2+⋯+Adeξ1θ1eξ^2θ2⋯eξ^n-1θn-1deξ^nθne-ξ^nθn,where deξ^iθie−ξ^iθi is(12)deξ^iθie−ξ^iθi=∫01eξ^iθisdξ^iθie−ξ^iθisdS+∫01eξ^iθisξ^,=θi∫01Adeξ^iθisdξ^ids+ξ^idθi.


### 3.2. Model of Transmission Errors

Transmission errors mainly include the reduction ratio error and coupling ratio error. The errors are caused by the process of manufacturing and assembly. Both the errors lead to positioning inaccuracy.Considering the transmission errors, the nominal value of theith joint angle can be given by(13)θin=θin′+∑j=1ncjinθjn′,ki=rinria=θim/θin′θim/θia′,hji=cjia−cjin,where rin and ria are the nominal and identified reduction ratios of the ith joint. cjin and cjia are the nominal and identified coupling ratios. θin′ and θia′ are the nominal and actual output angles of the ith joint reducer. θim is the ith joint motor rotation angle. ki and hji are the errors of the reduction ratio and the coupling ratio, namely the transmission errors.Then, theith joint angle θi can be expressed as(14)θi=θi0+θia′+∑j=1ncjiaθja∗=θi0+kiθin−ki∑j=1ncjinθjn′+∑j=1ncjinθja∗+∑j=1nhjiθja∗,=θi0+kiθin+∑j=1nkj−kicijnθjn+∑j=1nhijkjθjn,where θi0 is the zero positioning error of the joint angle θi, namely the ith joint offset and θj is the jth joint angle which is coupled with θi.Furthermore, the deviation of theith joint angle θi can be deduced as(15)dθi=dθi0+θin−∑j=1ncjinθjndki+∑j=1nkjθjndhjt.


### 3.3. Linearization of the Error Model

Linearization of the error model is a premise for parameters identification. The error model can be obtained by combining the two error models mentioned above as(16)dP=θ0I−P^O0∫01Adeξ^0θ0sdsdξ0+MI−P^cTO0dg⋅g−1V+MRgdPcr00,where(17)dg⋅g−1V=θ1∫01Adeξ^1θ1sdsdξ1+ξ1dθ1+Adeξ^1θ1θ2∫01Adeξ^2θ2Sdsdξ2+ξ2dθ2+⋯+Adeξ^1θ1eξ^2θ2…eξ^n−1θn−1θn∫01AdeξnθnSdsdξn+ξndθn,where eξ^iθi can be expressed as a homogeneous matrix through Taylor expansion and the adjoint transformation associated with eξ^iθi can be expressed as a 6 × 6 matrix. Equation (16) can be deduced as(18)dP=Ax,where(19)dP=Pa−Pn,where Pa is the actual TCP position measured by the laser tracker and Pn is the nominal TCP position obtained by the forward kinematics. The parameters to be identified can be expressed as(20)x=dξ0T,dξ1T…,dξnT…,dθn0…,dkn,dh21…,dhji,dPcTTT.A is the corresponding coefficient matrix which can be expressed as(21)A=Q0Q1…QnT1…TnK1⋯KnH21…HjiMRg.With the measurement of the TCP position at different configurations, we can get the following equation set:(22)dP˜=dP1⋯dPn=A1⋯Anx=A˜x,where n is the number of measurements.The least-squares solution forx is(23)x=A˜T⋅A˜−1⋅A˜T⋅dP˜.


## 3.1. Model of Geometric Parameters Errors

The position of the tool center point (TCP) in the measure frame is shown below:(3)P=M⋅PcT=M⋅gθ⋅PcT0,where P is the position of the TCP, M is a homogenous matrix that represents the transformation from the measuring device frame to the base frame, and PCT represents the position of TCP expressed in the base frame corresponding to θ=θ1,θ2,…,θn, gθ=eξ^1ϑ1eξ^2ϑ2⋯eξ^nϑn.Then, we have the deviation of the TCP position by taking the differential on both sides of equation (3):(4)dP=dM⋅g⋅Pcr0+M⋅dg⋅Pcr0+M⋅g⋅dPcr0,=dM⋅M−1⋅P+M⋅dg⋅g−1⋅Pcr+M⋅g⋅dPcr0,where(5)g⋅dPcT0=Rgbg00dPcT00=RgdPcT00.Taking a close observation of (4), we notice that dM × M−1 ∈ se(3), dg × g−1 ∈ se(3), and the element s belonging to lie algebra satisfies(6)sq∗=ω^Vo0q1=ω^q+υ0=I3−q^00υω.d
P can be expressed as(7)dP=I−P^00dM⋅M−1V+MI−P^cT00dg⋅g−1V+MRgdPcT00.Since rigid body motion can be realized by a rotation and a translation about an axis,M can be expressed as(8)M=eξ^ϑ0∈SE3.Take the differential on both sides:(9)dM⋅M−1=deξ^0ϑ0⋅e−ξ^0ϑ0=θ0∫01eξ^0ϑ0sdξ^0e−ξ^0ϑ0sds.The matrixg is(10)g=Rgbg01∈SE3.Take the differential on both sides:(11)dg⋅g-1=deξ^1θ1e-ξ^1θ1+eξ^1θ1deξ^2θ2e-ξ^2θ2e-ξ^1θ1+⋯+eξ^1θ1eξ^2θ2…eξ^n-1θn-1deξ^nθne-ξ^nθn…e-ξ^2θ2e-ξ^1θ1,=deξ^1θ1e-ξ^1θ1+Adeξ^1θ1deξ^2θ2e-ξ^2θ2+⋯+Adeξ1θ1eξ^2θ2⋯eξ^n-1θn-1deξ^nθne-ξ^nθn,where deξ^iθie−ξ^iθi is(12)deξ^iθie−ξ^iθi=∫01eξ^iθisdξ^iθie−ξ^iθisdS+∫01eξ^iθisξ^,=θi∫01Adeξ^iθisdξ^ids+ξ^idθi.


## 3.2. Model of Transmission Errors

Transmission errors mainly include the reduction ratio error and coupling ratio error. The errors are caused by the process of manufacturing and assembly. Both the errors lead to positioning inaccuracy.Considering the transmission errors, the nominal value of theith joint angle can be given by(13)θin=θin′+∑j=1ncjinθjn′,ki=rinria=θim/θin′θim/θia′,hji=cjia−cjin,where rin and ria are the nominal and identified reduction ratios of the ith joint. cjin and cjia are the nominal and identified coupling ratios. θin′ and θia′ are the nominal and actual output angles of the ith joint reducer. θim is the ith joint motor rotation angle. ki and hji are the errors of the reduction ratio and the coupling ratio, namely the transmission errors.Then, theith joint angle θi can be expressed as(14)θi=θi0+θia′+∑j=1ncjiaθja∗=θi0+kiθin−ki∑j=1ncjinθjn′+∑j=1ncjinθja∗+∑j=1nhjiθja∗,=θi0+kiθin+∑j=1nkj−kicijnθjn+∑j=1nhijkjθjn,where θi0 is the zero positioning error of the joint angle θi, namely the ith joint offset and θj is the jth joint angle which is coupled with θi.Furthermore, the deviation of theith joint angle θi can be deduced as(15)dθi=dθi0+θin−∑j=1ncjinθjndki+∑j=1nkjθjndhjt.


## 3.3. Linearization of the Error Model

Linearization of the error model is a premise for parameters identification. The error model can be obtained by combining the two error models mentioned above as(16)dP=θ0I−P^O0∫01Adeξ^0θ0sdsdξ0+MI−P^cTO0dg⋅g−1V+MRgdPcr00,where(17)dg⋅g−1V=θ1∫01Adeξ^1θ1sdsdξ1+ξ1dθ1+Adeξ^1θ1θ2∫01Adeξ^2θ2Sdsdξ2+ξ2dθ2+⋯+Adeξ^1θ1eξ^2θ2…eξ^n−1θn−1θn∫01AdeξnθnSdsdξn+ξndθn,where eξ^iθi can be expressed as a homogeneous matrix through Taylor expansion and the adjoint transformation associated with eξ^iθi can be expressed as a 6 × 6 matrix. Equation (16) can be deduced as(18)dP=Ax,where(19)dP=Pa−Pn,where Pa is the actual TCP position measured by the laser tracker and Pn is the nominal TCP position obtained by the forward kinematics. The parameters to be identified can be expressed as(20)x=dξ0T,dξ1T…,dξnT…,dθn0…,dkn,dh21…,dhji,dPcTTT.A is the corresponding coefficient matrix which can be expressed as(21)A=Q0Q1…QnT1…TnK1⋯KnH21…HjiMRg.With the measurement of the TCP position at different configurations, we can get the following equation set:(22)dP˜=dP1⋯dPn=A1⋯Anx=A˜x,where n is the number of measurements.The least-squares solution forx is(23)x=A˜T⋅A˜−1⋅A˜T⋅dP˜.


## 4. Joint Stiffness Identification

The stiffness of the manipulator arm is much greater than that of the joint. For this reason, the joint stiffness identification is investigated in this section.For a joint, the dynamic model can be described as(24)Mθθ¨+Cθ,θ˙+τgθ+τfθ˙=τm,where Mθ is the inertial matrix, Cθ,θ˙ is the Coriolis and centrifugal terms, τgθ relates to the gravitational torque, and τfθ˙ is the friction torque in the joint. τm is the input torque from the joint motor.On the assumption that the angular deformation of the reducer is proportional to the input torque, the relationship between the input torque and the deformation of the reducer is as follows.(25)τmi=ki∗δθki,where ki is the stiffness coefficient of the ith joint and δθki is the angular deformation of the ith joint. The joint torque τm is obtained as(26)τm=CeϕI,where I is the motor current that can be read directly from the controller of the robot. Ce is the motor potential constant, and ϕ is the magnetic flux.For a manipulator, the stiffness matrix of the robot joint can be obtained as follows by ignoring the interaction of each joint.(27)K=k1⋯0⋮⋱⋮0⋯kn,where n is the number of degrees of freedom.Substitute equation (25) into equation (24), and the relationship between joint angle deformation and motor control current can be given by(28)δθki=Ceϕki∗Ii.The flexibility matrix of the robot joint is as follows:(29)C=Ce1ϕ1k1⋯0⋮⋱⋮0⋯Cmϕnkn.According to the differential error relationship, the joint angle deviation vectorx=δθk1⋯,δθknT can be obtained.Since only the gravity load is loaded during the experiment, the direction of gravity is consistent with the direction of the first axis of the robot and the moment of gravity decomposition to the first axis is 0, and the first axis is not identified. Similarly, the flexibility coefficient of the sixth axis is not identified.Substitute the flexibility matrix into the equation below to improve the positioning accuracy by compensating the joint angles.(30)θ∗=θn+IC,where θ∗ is the modified joint angle value of the robot. θn is the nominal value of the joint angle, and the vector I is the current value of each joint motor.


## 5. MLPNN Based on a Hybrid Optimization Method of Beetle Antennae Search Algorithm and Particle Swarm Optimization for Nongeometric Errors Compensation

The positioning error is still large after the geometric parameters, transmission errors, and joint stiffness compensation. The residual positioning errors are caused by nongeometric errors which is difficult to model. In order to overcome the drawbacks, in this section, a MLPNN based on a new hybrid optimization method is proposed. This new model will be presented in details.


### 5.1. Beetle Antennae Search Algorithm (BAS)

The beetle antennae search algorithm (BAS) is a metaheuristic algorithm that is inspired by the searching behavior of longhorn beetles [35, 36]. A vector xt at tth time instant (t = 1, 2, …) denotes the position of the beetle. At position x, fx is defined as the fitness function to represent the concentration of odour. Normally, we can use two rules to simplify the algorithm, including search and detection behavior. In an unknown environment, the beetle searches randomly. A normalized random unit vector b⟶ is utilized to model the searching behavior as follows:(31)b⟶=randsm,1randsm,1,where rands(.) represents a random function and m is the number of dimensions of the position. xr and x1 are defined to imitate the searching behavior of both the right and left antennae:(32)xr=xt+dtb⟶,xl=xt+dtb⟶,where d is the sensing length of antennae.An iterative form is proposed to model the detecting behavior:(33)xt+1=xt+δtb⟶signfxr−fx1,where δ is the step size and sign(.) represents a sign function.The antennae lengthd and step size δ are updated according to the following formula.(34)dt=0.95dt−1+0.01,δt=0.95δt−1.Compared with other heuristic algorithms such as PSO, GA, and ABC, the BAS algorithm possesses better performance in terms of faster speed and simpler implementation. A second order Michalewicz function shown below was used to validate the algorithm.(35)fx=sinx1sinx12π20+sinx2sin2x22π20,where the minimized value satisfies fx∗ = −1.8013, locating in x∗ = (2.20, 1.57). The initial values of d and δ are 2 and 0.5. Taking the absolute error less than 0.0001 as the end evolution, the time these algorithms mentioned above takes is shown in Table 2.Table 2
Time comparisons of different algorithms in a two-dimensional situation.
Algorithms
BAS
PSO
GA
ABC
T (s)
0.087843
0.215134
0.262317
0.132563The results showed that the search speed of BAS algorithm was tremendous. The BAS algorithm can fit well for the real-time task.


### 5.2. Beetle Swarm Optimization (BSO)

The performance of beetle antennae search algorithm relies heavily on the initial position, and it will easily fall into the local optimum when dealing with high-dimensional problems [37–39]. In order to overcome these drawbacks, further improvements are made by combining the antennae search algorithm and particle swarm optimization. That is the beetle swarm optimization algorithm.In beetle swarm optimization, each particle in the standard PSO is defined as a beetle. Each beetle represents a potential solution to the optimization problem. In the iterative process, the update of the position relies not only on the individual and global historical best solution but also on the way of beetle antennae search. The formula can be expressed as follows [34]:(36)Xisk+1=Xisk+λVisk+1−λξisk,where i is the ith beetle, s is the sth dimension, k is the current number of iterations, Vis is the speed of the beetle, ξis is the increase in movement, and λ is a positive constant. Vis is expressed as:(37)Visk+1=ωVisk+c1rand1pisk−Xisk+c2rand2pgsk,where w is the inertia weight, c1 and c2 are constants, and rand1 and rand2 are two random functions in (0, 1). Pi and Pg are the individual and group extremities.The inertia weight is decreasing in the process as follows:(38)ω=ωmax−ωmax−ωminK∗k,where ωmax and ωmin represent the maximum and minimum values of ω and k and K are the current number of iterations and the maximum number of iterations.ξ defines the increase in movement and is expressed as(39)ξisk+1=δk∗Visk∗signfXrsk−fXrsk,where δ is the step size. The position of the left and right antenna can be expressed as(40)Xrsk+1=Xrsk+ViSk∗d2,Xlsk+1=Xlsk+Visk∗d2.


### 5.3. MLPNN

Artificial neural network and, in particular, the multilayer perceptron neural networks are widely used in many application areas over the years. Typically, the multilayer perceptron neural network contains three kinds of layers, which are the input layer, hidden layer, and output layer. It is claimed that a three-layered feed-forward neural network can approximate any nonlinear function with arbitrary accuracy. The number of neurons in input and output layers is normally selected according to the actual needs. However, the best number of neurons in the hidden layer can just be determined by trial-and-error methods.Generally, each neuron can be denoted by the following equation:(41)yj=f∑i=1nωijxi+bj,where xi and yj are the input and output values of the jth neuron, respectively; ωij denotes the connection weight from the ith neuron in the previous layer to the jth neuron in the latter layer; bj represents bias value in the jth neuron; and f is normally a sigmoid function as follows:(42)fx=11+e−x.Its derivation is(43)f′x=fx1−fx.The error of thekth neuron in the output layer is(44)ek=Yk−0k,where Yk and OK are the desired and actual value of the ith neuron in the output layer. The total error of the output layer is(45)E=12∑k=1mek2,where m is the number of neurons in the output layer.For nongeometric errors compensation, the joint angles are selected as the inputs and the positioning errors in each axis are selected as the outputs. By using the Jacobian matrix, the positioning errors are mapped to the errors of joint angles. The compensation process is shown in Figure2.Figure 2
NN-based nongeometric errors compensation.


### 5.4. The MLPNN Optimized by BSO

The validity and accuracy of MLFNN may be reduced if the weights and biases are improperly selected. The BSO algorithm is used to optimize the network parameters of MLFNN in this process. In this section, the BSO-MLFNN method is proposed to overcome the shortcomings of low accuracy. The process is shown in Figure3.Figure 3
Algorithm flow chart of BSO-MLPNN.


## 5.1. Beetle Antennae Search Algorithm (BAS)

The beetle antennae search algorithm (BAS) is a metaheuristic algorithm that is inspired by the searching behavior of longhorn beetles [35, 36]. A vector xt at tth time instant (t = 1, 2, …) denotes the position of the beetle. At position x, fx is defined as the fitness function to represent the concentration of odour. Normally, we can use two rules to simplify the algorithm, including search and detection behavior. In an unknown environment, the beetle searches randomly. A normalized random unit vector b⟶ is utilized to model the searching behavior as follows:(31)b⟶=randsm,1randsm,1,where rands(.) represents a random function and m is the number of dimensions of the position. xr and x1 are defined to imitate the searching behavior of both the right and left antennae:(32)xr=xt+dtb⟶,xl=xt+dtb⟶,where d is the sensing length of antennae.An iterative form is proposed to model the detecting behavior:(33)xt+1=xt+δtb⟶signfxr−fx1,where δ is the step size and sign(.) represents a sign function.The antennae lengthd and step size δ are updated according to the following formula.(34)dt=0.95dt−1+0.01,δt=0.95δt−1.Compared with other heuristic algorithms such as PSO, GA, and ABC, the BAS algorithm possesses better performance in terms of faster speed and simpler implementation. A second order Michalewicz function shown below was used to validate the algorithm.(35)fx=sinx1sinx12π20+sinx2sin2x22π20,where the minimized value satisfies fx∗ = −1.8013, locating in x∗ = (2.20, 1.57). The initial values of d and δ are 2 and 0.5. Taking the absolute error less than 0.0001 as the end evolution, the time these algorithms mentioned above takes is shown in Table 2.Table 2
Time comparisons of different algorithms in a two-dimensional situation.
Algorithms
BAS
PSO
GA
ABC
T (s)
0.087843
0.215134
0.262317
0.132563The results showed that the search speed of BAS algorithm was tremendous. The BAS algorithm can fit well for the real-time task.


## 5.2. Beetle Swarm Optimization (BSO)

The performance of beetle antennae search algorithm relies heavily on the initial position, and it will easily fall into the local optimum when dealing with high-dimensional problems [37–39]. In order to overcome these drawbacks, further improvements are made by combining the antennae search algorithm and particle swarm optimization. That is the beetle swarm optimization algorithm.In beetle swarm optimization, each particle in the standard PSO is defined as a beetle. Each beetle represents a potential solution to the optimization problem. In the iterative process, the update of the position relies not only on the individual and global historical best solution but also on the way of beetle antennae search. The formula can be expressed as follows [34]:(36)Xisk+1=Xisk+λVisk+1−λξisk,where i is the ith beetle, s is the sth dimension, k is the current number of iterations, Vis is the speed of the beetle, ξis is the increase in movement, and λ is a positive constant. Vis is expressed as:(37)Visk+1=ωVisk+c1rand1pisk−Xisk+c2rand2pgsk,where w is the inertia weight, c1 and c2 are constants, and rand1 and rand2 are two random functions in (0, 1). Pi and Pg are the individual and group extremities.The inertia weight is decreasing in the process as follows:(38)ω=ωmax−ωmax−ωminK∗k,where ωmax and ωmin represent the maximum and minimum values of ω and k and K are the current number of iterations and the maximum number of iterations.ξ defines the increase in movement and is expressed as(39)ξisk+1=δk∗Visk∗signfXrsk−fXrsk,where δ is the step size. The position of the left and right antenna can be expressed as(40)Xrsk+1=Xrsk+ViSk∗d2,Xlsk+1=Xlsk+Visk∗d2.


## 5.3. MLPNN

Artificial neural network and, in particular, the multilayer perceptron neural networks are widely used in many application areas over the years. Typically, the multilayer perceptron neural network contains three kinds of layers, which are the input layer, hidden layer, and output layer. It is claimed that a three-layered feed-forward neural network can approximate any nonlinear function with arbitrary accuracy. The number of neurons in input and output layers is normally selected according to the actual needs. However, the best number of neurons in the hidden layer can just be determined by trial-and-error methods.Generally, each neuron can be denoted by the following equation:(41)yj=f∑i=1nωijxi+bj,where xi and yj are the input and output values of the jth neuron, respectively; ωij denotes the connection weight from the ith neuron in the previous layer to the jth neuron in the latter layer; bj represents bias value in the jth neuron; and f is normally a sigmoid function as follows:(42)fx=11+e−x.Its derivation is(43)f′x=fx1−fx.The error of thekth neuron in the output layer is(44)ek=Yk−0k,where Yk and OK are the desired and actual value of the ith neuron in the output layer. The total error of the output layer is(45)E=12∑k=1mek2,where m is the number of neurons in the output layer.For nongeometric errors compensation, the joint angles are selected as the inputs and the positioning errors in each axis are selected as the outputs. By using the Jacobian matrix, the positioning errors are mapped to the errors of joint angles. The compensation process is shown in Figure2.Figure 2
NN-based nongeometric errors compensation.


## 5.4. The MLPNN Optimized by BSO

The validity and accuracy of MLFNN may be reduced if the weights and biases are improperly selected. The BSO algorithm is used to optimize the network parameters of MLFNN in this process. In this section, the BSO-MLFNN method is proposed to overcome the shortcomings of low accuracy. The process is shown in Figure3.Figure 3
Algorithm flow chart of BSO-MLPNN.


## 6. Compensation Verification and Analysis

As shown in Figure4, a spherical mounted retroreflector (SMR) was mounted at the end flange of the manipulator, and the position of the SMR, namely, the tool center point (TCP) position, was measured by a laser tracker (FARO Vantage) with the positioning accuracy of 10 μm + 2.5 μm/m. Experiments were performed on a heavy load SIASUN SR210D robot.Figure 4
Experimental scene.The main steps of the experiment were as follows: (1) a set of 50 robot joint angle values in the controller and the corresponding positions without any load were obtained in the measuring coordinate system. Geometric parameters and transmission errors were identified using the improved POE model. (2) Identify the joint stiffness by measuring 50 points with a 210 kg load after geometric parameters and transmission errors compensation. (3) Compensate the nongeometric errors based on the proposed BSO-MLPNN by measuring 200 points after the two steps mentioned above.The measurement points should be selected evenly in the workspace, so that the measuring configurations should go through all controllable DoFs and close to the boundary of the workspace, where the maximum errors are most prominent.


### 6.1. Geometric Parameters and Transmission Errors Identification

According to equation (22), the actual geometric parameters, transmission errors and coupling ratio can be identified. The nominal and identified value of geometric parameters and reduction ratio are shown in Tables 3 and 4.Table 3
Nominal value and identified value of robot geometric parameters.
Nominal parameters
L
i (mm), i = 1–5
312
1075
235
1282
260
0
θ
i
0 (°), i = 1–6
0
0
0
0
0
0
Identified parameters
L
i
′ (mm), i = 1–5
313.1011
1075.1209
234.2737
1281.3124
260
0
θ
i
0
′ (°), i = 1–6
0
0.0044
0.0479
0.4524
−0.1892
0Table 4
Nominal and identified value of reduction ratio errors.
Nominal reduction ratio errors
k
i, i = 1–6
1
1
1
1
1
1
Identified reduction ratio errors
k
i
′, i = 1–6
1.0002
1.0013
1.0004
0.9998
0.9996
1.0003The nominal coupling ratio error is 0 and the identified coupling ratio error is−2.5342×10−4. The positioning accuracy of the 50 points before and after calibration is shown in Figure 5. The average positioning accuracy of measurement points is enhanced from 1.1288 mm to 0.2898 mm. For validation, another set of 50 points is selected. The validation results show that the average positioning accuracy is 0.3008 mm which is 0.011 mm larger than 0.2898 mm. Details are presented in Table 5.Figure 5
Positioning accuracy before and after geometric calibration.Table 5
Absolute positioning accuracy of the SIASUN SR210D robot.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before calibration
1.1288
0.4709
2.3640
After robot link geometry and transmission errors compensation of measurement points
0.2898
0.1464
0.8226
After robot link geometry and transmission errors compensation of validation points
0.3008
0.1532
0.8436


### 6.2. Stiffness Identification and Compensation

According to equation (28), the joint stiffness is identified. The stiffness results are shown in Table 6. The positioning accuracy of the 50 points before and after calibration is shown in Figure 6. The average positioning accuracy is enhanced from 7.0261 mm to 0.9847 mm. Details are presented in Table 7.Table 6
Flexibility factors for stiffness calibration identification.
c
ii, i = 1–6
0
1.89e − 5
−2.49e − 5
−5.35e − 5
−8.26e − 5
0Figure 6
Positioning accuracy after the calibration of stiffness parameters under the full load condition.Table 7
Absolute positioning accuracy of the SIASUN SR210D manipulator.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before joint stiffness compensation
7.0261
1.2351
10.0898
After joint stiffness compensation
0.9847
0.5340
2.4176


### 6.3. Nongeometric Errors Compensation Based on BSO-MLPNN

Figure6 shows that the positioning errors are still large after the geometric parameters, transmission errors, and the joint stiffness compensation. These errors caused by nongeometric sources are compensated by the BSO-MLPNN method. For this purpose, 200 evenly distributed points in the workspace and corresponding joint angles are collected for training the BSO-MLPNN. In the training process, joint angles are selected as the inputs and the residual positioning errors are selected as the outputs. The residual positioning errors in each axis are shown in Figure 7.Figure 7
Residual positioning errors in each axis.The positioning accuracy before and after nongeometric calibration at the measurement points is shown in Figure8. The results show that, after the nongeometric errors compensation at the measurement points, the average positioning accuracy is enhanced from 0.9262 (mm) to 0.5907 (mm). Details are shown in Table 8.Figure 8
Absolute positioning accuracy before and after nongeometric compensation of measurement points.Table 8
Absolute positioning accuracy of the SIASUN SR210D manipulator of measurement points.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before BSO-MLPNN calibration
0.9262
0.1632
1.7435
After BSO-MLPNN calibration
0.5907
0.0760
0.8834Normally, the positioning accuracy of the measurement points is higher than that of the validation points. For validation, another set of 200 points are selected. The positioning errors are shown in Figure9. The validation results show that the average positioning accuracy is enhanced from 0.9370 mm to 0.6522 mm. Details are shown in Table 9.Figure 9
Absolute positioning accuracy before and after nongeometric errors compensation of validation points.Table 9
Absolute positioning accuracy of the SIASUN SR210D manipulator of validation points.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before BSO-MLPNN calibration
0.9370
0.1871
1.7633
After BSO-MLPNN calibration
0.6522
0.0887
0.9834Therefore, the proposed three stages calibration method including geometric parameters, transmission errors, joint stiffness, and nongeometric errors compensation is effective.


## 6.1. Geometric Parameters and Transmission Errors Identification

According to equation (22), the actual geometric parameters, transmission errors and coupling ratio can be identified. The nominal and identified value of geometric parameters and reduction ratio are shown in Tables 3 and 4.Table 3
Nominal value and identified value of robot geometric parameters.
Nominal parameters
L
i (mm), i = 1–5
312
1075
235
1282
260
0
θ
i
0 (°), i = 1–6
0
0
0
0
0
0
Identified parameters
L
i
′ (mm), i = 1–5
313.1011
1075.1209
234.2737
1281.3124
260
0
θ
i
0
′ (°), i = 1–6
0
0.0044
0.0479
0.4524
−0.1892
0Table 4
Nominal and identified value of reduction ratio errors.
Nominal reduction ratio errors
k
i, i = 1–6
1
1
1
1
1
1
Identified reduction ratio errors
k
i
′, i = 1–6
1.0002
1.0013
1.0004
0.9998
0.9996
1.0003The nominal coupling ratio error is 0 and the identified coupling ratio error is−2.5342×10−4. The positioning accuracy of the 50 points before and after calibration is shown in Figure 5. The average positioning accuracy of measurement points is enhanced from 1.1288 mm to 0.2898 mm. For validation, another set of 50 points is selected. The validation results show that the average positioning accuracy is 0.3008 mm which is 0.011 mm larger than 0.2898 mm. Details are presented in Table 5.Figure 5
Positioning accuracy before and after geometric calibration.Table 5
Absolute positioning accuracy of the SIASUN SR210D robot.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before calibration
1.1288
0.4709
2.3640
After robot link geometry and transmission errors compensation of measurement points
0.2898
0.1464
0.8226
After robot link geometry and transmission errors compensation of validation points
0.3008
0.1532
0.8436


## 6.2. Stiffness Identification and Compensation

According to equation (28), the joint stiffness is identified. The stiffness results are shown in Table 6. The positioning accuracy of the 50 points before and after calibration is shown in Figure 6. The average positioning accuracy is enhanced from 7.0261 mm to 0.9847 mm. Details are presented in Table 7.Table 6
Flexibility factors for stiffness calibration identification.
c
ii, i = 1–6
0
1.89e − 5
−2.49e − 5
−5.35e − 5
−8.26e − 5
0Figure 6
Positioning accuracy after the calibration of stiffness parameters under the full load condition.Table 7
Absolute positioning accuracy of the SIASUN SR210D manipulator.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before joint stiffness compensation
7.0261
1.2351
10.0898
After joint stiffness compensation
0.9847
0.5340
2.4176


## 6.3. Nongeometric Errors Compensation Based on BSO-MLPNN

Figure6 shows that the positioning errors are still large after the geometric parameters, transmission errors, and the joint stiffness compensation. These errors caused by nongeometric sources are compensated by the BSO-MLPNN method. For this purpose, 200 evenly distributed points in the workspace and corresponding joint angles are collected for training the BSO-MLPNN. In the training process, joint angles are selected as the inputs and the residual positioning errors are selected as the outputs. The residual positioning errors in each axis are shown in Figure 7.Figure 7
Residual positioning errors in each axis.The positioning accuracy before and after nongeometric calibration at the measurement points is shown in Figure8. The results show that, after the nongeometric errors compensation at the measurement points, the average positioning accuracy is enhanced from 0.9262 (mm) to 0.5907 (mm). Details are shown in Table 8.Figure 8
Absolute positioning accuracy before and after nongeometric compensation of measurement points.Table 8
Absolute positioning accuracy of the SIASUN SR210D manipulator of measurement points.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before BSO-MLPNN calibration
0.9262
0.1632
1.7435
After BSO-MLPNN calibration
0.5907
0.0760
0.8834Normally, the positioning accuracy of the measurement points is higher than that of the validation points. For validation, another set of 200 points are selected. The positioning errors are shown in Figure9. The validation results show that the average positioning accuracy is enhanced from 0.9370 mm to 0.6522 mm. Details are shown in Table 9.Figure 9
Absolute positioning accuracy before and after nongeometric errors compensation of validation points.Table 9
Absolute positioning accuracy of the SIASUN SR210D manipulator of validation points.
Mean (mm)
Standard deviation (mm)
Maximum (mm)
Before BSO-MLPNN calibration
0.9370
0.1871
1.7633
After BSO-MLPNN calibration
0.6522
0.0887
0.9834Therefore, the proposed three stages calibration method including geometric parameters, transmission errors, joint stiffness, and nongeometric errors compensation is effective.


## 7. Conclusion and Future Works

This paper proposed a three-stage calibration method for enhancing robot positioning accuracy. The geometric parameters and transmission errors are identified in the first stage. Then, joint stiffness compensation is carried out. For the nongeometric errors compensation, a hybrid algorithm, namely, BSO-MLPNN, is introduced.Experiments were performed on the SIASUN SR210D robot manipulator. The average positioning accuracy was enhanced from 7.0261 mm to 0.5907 mm. Results demonstrated the effectiveness and correctness of the proposed method. In addition, the validation results showed that the manipulator after calibration has the same level of positioning accuracy in the entire workspace.In future, nongeometric factors will be modeled in order to obtain more accurate knowledge of error sources and accelerate the convergence rate.



---

*Source: 5896813-2020-01-16.xml*