# Positioning Error Compensation for Industrial Robots Based on Stiffness Modelling

Authors: Yingjie Li; Guanbin Gao; Fei Liu
Journal: Complexity
 (2020)

---

## Abstract

Insufficient stiffness of industrial robots is a significant factor which affects its positioning accuracy. To improve the positioning accuracy, a novel positioning error compensation method based on the stiffness modelling is proposed in this paper. First, the positioning errors considering the end load and gravity of industrial robots due to stiffness are analyzed. Based on the results of analysis, it is found that the positioning errors can be described by two kinds of deformation errors at joints: the axial deformation error and the radial deformation error. Then, the axial deformation error is modelled by the differential relationship of kinematics equations. The model of radial deformation error is deduced through the recurrence method and rotation transformation between joints. Finally, these two models are transformed into a Cartesian coordinate system, and a positioning error compensation method based on these two models is presented. Simulations based on the finite element analysis are implemented to verify the positioning error compensation method. The results show that the suggested method can efficiently predict the positioning error according to the gravity and loads, so that the positioning accuracy of industrial robots can be improved with the proposed method.


---

## Body

## 1. Introduction

With the rapid development of intelligent manufacturing, industrial robots (noted as robots) have been widely applied in automobile manufacturing, logistics systems, mechanical processing, food packaging industries, etc. [1]. According to the report of the International Federation of Robotics (IFR), more than 72% of industrial robots are used in the low-precision occasions, e.g., sorting, palletizing, handling, spot welding, painting, and assembly of simple parts [2] for the positioning, and the trajectory accuracy of industrial robots is still relatively low [3]. Kinematic calibration is the general way to improve the absolute positioning and trajectory accuracy of industrial robots [4]. However, kinematic calibration can only compensate the static kinematic parameter errors of industrial robots, while the dynamic factors, e.g., loads of the end effector (noted as EE), speed, acceleration, gravity, and poses, also influence the positioning accuracy greatly. One of the most notable characteristics of industrial robots is the open-chain structure that provides high flexibility and large working space for industrial robots [5, 6]. However, the open-chain structure also results in a low stiffness and an error accumulation amplification which are the main reasons for the low accuracy of industrial robots. Therefore, improvement of the stiffness error has always been an important research field for industrial robots. At present, there are mainly three methods to improve the stiffness, i.e., structural strengthening method, stiffness control method, and working space compensation method.The structural strengthening is an approach to improve the absolute accuracy of robots, which enhances the integral stiffness of robots by changing the material properties of the robot’s mechanical parts and transmission parts. This approach requires the designers to consider the low positioning accuracy caused by insufficient stiffness in the design. Hence, the designers need to estimate the robot’s stiffness in advance according to its application fields. Tyapin et al. [7] proposed a method to model the physical stiffness of the driver, reducer, and transmission parts of the robots, respectively. With the corresponding weight coefficients obtained from the model, the robot’s stiffness can be further integrated and estimated. However, the stiffness model is not precise enough according to its experimental results. Liu et al. [8] presented a method of modelling the rotating joint of the robots based on the analysis of the contact relationships of the robot’s joints. Although one can calculate the stiffness of each joint of the robots accurately with this method, its calculation process is complicated and the consumption cost of calculation is considerable. By using structural strengthening, the stiffness of the robots can be improved in the stage of manufacture. However, after manufacture the stiffness of robots determined by the material properties and structure will not be changed.The stiffness control method aims to control the robots at some certain poses in which the stiffness of robots is relatively high [9, 10]. Abele et al. [11, 12] presented an adaptive machining method after measuring workpiece shape to overcome the problem of low positioning accuracy caused by the weak stiffness when machining. The results show that the machining accuracy of robots can be improved by changing the robot’s poses to adjust the torque at the joint, which can enhance the robot’s stiffness indirectly. By taking the stiffness ellipsoid as the index to assess stiffness, Guo et al. [13] improved the robot’s stiffness by optimizing its working poses with the maximum joint angle as the constraints. Combining with the redundancy, Zanchettin et al. [14] presented a method of optimizing the robot’s poses to enhance the stiffness when the robots perform drilling tasks. From the above methods, it can be found that the main idea of the stiffness control method is to strengthen the EE operation stiffness of the robots by changing or optimizing the robot’s working poses [15]. Nevertheless, the precondition of optimizing the pose to improve the positioning accuracy is that there is a contact force at the EE of robots. Therefore, when there are no contact forces, e.g., painting, prealignment, and measurement, the stiffness control method does not work.Based on the kinematic and dynamic relationship of robots, the working space compensation methods can establish a mathematical model that can describe the relationship between the positioning error and end loads along with the gravity of the robot. And it can be used to predict the positioning error in the working space to improve the positioning accuracy. Salisbury proposed the traditional stiffness model of robots based on the kinematic and static theories [16]. The stiffness modelling for robots in the Cartesian coordinates was studied by connecting the stiffness of each joint [17–19]. Abele et al. [11] provided a stiffness model without calculating the inverse of Jacobian matrix, which can simplify the calculation of error compensation. In [20], the stiffness matrix of joints was identified by measuring the EE displacement and rotation of robots. Then, the stiffness matrix in the Cartesian coordinates was derived by the stiffness mapping model from the joint to EE. Sun et al. [21] proposed a method to calculate the EE translation stiffness for serial robots. However, according to the experimental results, this method is not accurate in terms of prediction of the positioning error. Overall, working space compensation methods are simpler and more universal than structural strengthening methods and stiffness control methods in practice, while the low accuracy limits its application. The main reason is that the current methods only consider the rotary deformation around the axis of joints, which does not contain the rotary deformation around the radial direction of joints.To improve the positioning accuracy, a novel positioning error compensation model based on the stiffness for industrial robots is proposed in this paper, which describes the relationship between the positioning error and the EE load and gravity of robots. First, according to the Newton–Euler method, the driving torque for each joint is calculated. Through the deformation assumption of joints, the torque is connected with the rotation deformation at joints. Moreover, the torque is also decomposed along the axial and radial directions of joints, which is consistent with the axial deformation and radial deformation. Positioning error models including the axial and radial deformation are established by means of kinematics differential and recurrence methods, respectively. In accordance with the small deformation assumption, two kinds of positioning error models are linearized. Finally, a complete error compensation model is derived according to the error models of axial and radial deformation, which can effectively predict and compensate the positioning error and improve the positioning accuracy for industrial robots.The major contributions of this paper include the following:(1)
A novel positioning error model is proposed forn-degree-of-freedom (DOF) industrial robots based on the relationship between its kinematics and dynamics parameters, which can be applied to arbitrary multi-DOF serial robots(2)
The Newton–Euler method is introduced to calculate the balance torque of joints, which makes it convenient to calculate the positioning error caused by the EE loads and gravity of the robot(3)
The radial deformation error at joints is modelled and included into the stiffness error model, which improves the accuracy of prediction and compensation for positioning error compared to the traditional methods(4)
The positioning error model is linearized by introducing proper assumptions, which reduces the complexity of the proposed method and makes it convenient to applicationsThe rest of the paper is organized as follows: in Section2, the displacement and deformation at joints of robots are analyzed with the gravity and EE load. In Section 3, six basic assumptions are provided, and the positioning error of robots is derived. In Section 4, the finite element analysis and three-dimensional model of a general six-DOF industrial robot are used to validate the effectiveness and correctness of the proposed method. Finally, conclusions for this paper are given in Section 5.Notations. Throughout the paper, ℝ denotes the set of real numbers. ℝn is the Euclidean space of n-dimensional real vectors. ℝm×n is the space of m×n real matrix. In×n denotes the identity matrix in ℝn×n. {0} stands for the base coordinate system of the robot. {i} is the coordinate system of link i of the robot. Tji∈ℝ4×4 is the homogeneous transform from {i} to {j}. Rji∈ℝ3×3 is the rotation transform from {i} to {j}. Pji∈ℝ3 represents the description of coordinate origin of {j} in {i}. PCji∈ℝ3 denotes the centroid position of link j of robots in {i}. fi∈ℝ3 stands for forces of link i-1 acting on link i. ni∈ℝ3 stands for torque of link i-1 acting on link i. Analogously fji∈ℝ3 and nji∈ℝ3 are descriptions of fj and nj in {i}, respectively. θ˙i∈ℝ,θ¨i∈ℝ are the angular speed and angular acceleration of link i relative to {i}. Z^ji∈ℝ3 is a vector of z-axis of {j} in {i}. ωji∈ℝ3 and ω˙ji∈ℝ3 represent the angular speed and angular acceleration of link j relating to {0} in {i}. v˙ji∈ℝ3 is defined as the acceleration of {j} in {i}. v˙Cji∈ℝ3 denotes the acceleration of centroid of link j in {i}. ⋅ means the vector norm of a vector in ℝn. mi+1 means the mass of link i+1 of industrial robots. To simplify expressions, cosθ, sinθ, and 1−sinθ are replaced by cθ, sθ, and vθ , respectively. ∂·/∂x means the gradient operator to variable x. S·∈ℝ3×3 is defined as(1)S·=0−VzVyVz0−Vx−VyVx0,where Vx∈ℝ, Vy∈ℝ, and Vz∈ℝ represent the three components of a vector. Further, S˜·∈ℝ4×4 can be defined as(2)S˜·=S·3×30004×4.The sign function is defined assign·.


## 2. Analysis of Deformation at Joints

A general industrial robot is shown in Figure1. When forces are applied to the robot and the link’s gravity is considered, balance forces and torques are generated at each joint and link. According to the theories of material mechanics, if the object is subjected to forces and torques, its shape will be changed, e.g., tension, compression, shear, torsion, and bending [22]. The effect of applied forces at joints and links of the robot is more complex than the results of theoretical analysis because the joint of robots consists of many elements in reality, e.g., motors, drive shafts, gears, and reducers. Hence, the deformation of robots in reality is a combination of the above five deformations. However, since the deformation at joints is significantly larger than the deformation at links, the deformation at joints is mainly studied in this paper.Figure 1
Two types of deformation at the joint. (a) The rotary deformation at the joint. (b) The linear deformation at the joint.
(a)(b)In Figure1, two types of deformation are shown, which include the rotary deformation around the axis of joints and the linear deformation along the certain direction. The joint of the robot deflects angle Δφ around its axis and its linear deformation is Δl under effects of gravity G and the end load F.According to the Newton–Euler method [23], the relationship between the motion and the driving force (or driving torque) of industrial robots can be described through the following dynamic equations:(3)ω˙i+1i+1=Rii+1ω˙ii+Rii+1ω˙ii×θ˙i+1Z^i+1i+1+θ¨i+1Z^i+1i+1,v˙i+1i+1=Rii+1ω˙ii×Pi+1i+ω˙ii×ω˙ii×Pi+1i+v˙ii,v˙Ci+1i+1=ω˙i+1i+1×PCi+1i+1+ωi+1i+1×ωi+1i+1×PCi+1i+1+v˙i+1i+1,Fi+1i+1=mi+1v˙Ci+1i+1,Ni+1i+1=Ii+1Ci+1ω˙i+1i+1+ωi+1i+1×Ii+1Ci+1ωi+1i+1,(4)fii=Ri+1ifi+1i+1+Fii,(5)nii=Nii+Ri+1ini+1i+1+PCii×Fii+Pi+1ii×Ri+1ifi+1i+1,where Fii∈ℝ3 and Nii∈ℝ3 are the inertia force and torque of link i in {i}, respectively. Ii+1Ci+1∈ℝ3×3 denotes the inertial tensor matrix of link i+1 in Ci+1 that is the coordinate system of centroid of link i+1. When the robot stops at a position, some variables in equations (4) and (5) are constant, e.g., θ˙i+1=0, θ¨i+1=0, ωi+1i+1=0,0,0T, and ω˙i+1i+1=0,0,0Ti+1. Thus, the above dynamic equations are further simplified as(6)v˙i+1i+1=Rii+1v˙ii,v˙Ci+1i+1=v˙i+1i+1,Fi+1i+1=mi+1v˙Ci+1i+1,(7)fii=Ri+1ifi+1i+1+Fii,(8)nii=Ri+1ini+1i+1+PCiinii×Fii+Pi+1i×Ri+1ifi+1i+1.It should be noted that although the robot is at stationary state, the value ofv˙ii is not zero. According to the weak principle of equivalence, the gravity applied to the robot is equivalent to a case that the robot has an initial acceleration, which is opposite to the gravity direction. In this paper, the gravity of the robot is considered, and v˙00 is a three-dimensional vector that is opposite to the gravity.From equations (7) and (8), conclusions can be obtained that even if the robot is at stationary state, the balance force and torque exist at the joints of the robot. To display the force and torque clearly, the joint is taken out from the robot as a separate body in Figure 2, in which F and M represent the resistant force and torque, respectively. Fd and Md denote the driving force and torque, respectively. Ma and Mr stand for axial and radial components of torque M. Thus, under the influence of force F and torque M, there will be two kinds of deformation at the joint, i.e., the rotary deformation around the torque’s axis and the linear deformation along the force direction.Figure 2
The force and torque at the joint. (a) The torque and force at the joint. (b) The decomposition of torqueM.
(a)(b)From the above analysis, there are two types of deformation at the joint under influence of the gravity and end load. Although the deformation at the joint may be very small in practice, the joint, as vital but weak parts of the robot, will greatly affect the positioning accuracy of the robots. It can be observed from (a) and (b) in Figure1 that the rotary deformation and linear deformation (these two deformations are usually appearing together) lead to the positioning error of the robots. Hence, it is essential to describe the positioning error caused by joint deformation with a mathematical model and to more accurately predict the positioning error. Hence, to resolve these problems, a generalized mathematical model for industrial robots is proposed to predict the positioning error in this paper.


## 3. Positioning Error Modelling for Industrial Robots

Some assumptions should be introduced before establishing the positioning error model because they are the basis of the proposed method.Assumptions are as follows:(i)
The industrial robots only contain rotating joints but not moving joints.(ii)
The elastic deformation of the robot’s link is negligible compared to its deformation at joints.(iii)
The effect caused by the rotary deformation at joints on the positioning error of the robots is much larger than the effect caused by the linear deformation of the joint on the end positioning error.(iv)
The rotary angleΔφ∈ℝ caused by the rotary deformation at the joint is small enough so that the following equations can be regarded as meaningful within the allowable range of accuracy:(9)sinΔφ=Δφ,cosΔφ=1.(v)
There is a linear relationship between the rotary deformationΔφ at the joint and the torque n∈ℝ3 applied to the joint, as shown below:(10)Δφ=Cn,whereC∈ℝ is the flexibility coefficient of the joint. Therefore, the stiffness coefficient can be defined as K=1/C.(vi)
There are two types of rotary deformation at each joint, i.e., the rotary deformationΔθ∈ℝ around joint’s axis and the rotary deformation Δγ∈ℝ around the radial direction of the joint. According to assumption (iv), the following equations can be obtained:(11)Δθ=Cana,Δγ=Crnr,where Ca∈ℝ denotes the axial stiffness coefficient of the joint and Cr∈ℝ denotes the radial stiffness coefficient of the joint. With assumption (v), we can define Ka=1/Ca and Kr=1/Cr. na∈ℝ3 and nr∈ℝ3 stand for axial torque and radial torque at the joint.Remark 1.
In practice, general robots consist of six revolute joints [24, 25]. Thus, assumption (i) is appropriate for general-purpose robots. It is also shown that the deformation of the joint due to the insufficient stiffness of driving and transmission system accounts for 70% of the total deformations which are caused by the external load or gravity [26]. Accordingly, assumption (ii) is true in this paper. As shown in Figure 1, although the rotary deformation and linear deformation at the joint may be tiny, the effects caused by the rotary deformation on the positioning error are significant because of the magnifying effect of the link. For this reason, assumption (iii) is reasonable. In accordance with [27, 28], the joint stiffness of industrial robot is 1 N/μm. In other words, a force of 1000 N is required to generate deformation of 1 mm. However, the maximum end load of most industrial robots is less than 1000 N. When the deformation at the joint is less than 1 mm, it can be reckoned that assumption (iv) is meaningful. In the light of Hooke’s law of the material, assumption (v) is feasible. In the previous discussions, there is a resistant torque M at the joint. The vector of torque M can be further decomposed along two directions, i.e., the axial and radial direction of joint, as shown in (b) of Figure 2. Hence, it can be considered that the rotary deformation consists of rotary deformation around joint’s axis and the radial direction of the joint. Since these two kinds of deformation are different in essence, Ca and Cr are required to describe the relationship between the rotary deformation and torque applied to the joint. Thus, assumption (vi) can also be valid.
Now, it is considered that the forceFe∈ℝ3 is applied to the EE of the robot and the gravity of the robot is also included. When the robot is stationary, the driving torque nii of joint i can be obtained by equations (6) to (8). This torque can be decomposed as(12)nii=naii+nrii,
Thus, the rotary deformations around axial and radial direction at jointi are(13)Δθi=Cainaii,Δγi=Crinrii.
According to the two types of rotary deformations at the joint in equation (13), the positioning error model of axial and the radial deformation will be established, respectively, in the next sections, and the total positioning error model will be derived finally.


### 3.1. Positioning Error Model of Axial Deformation around the Joint

When the EE of the robot is at a pointPe0 in Cartesian coordinate system and its coordinates are Pe00=px0,py0,pz0T in {0}. The joint angle corresponding to Pe00 is Θ=θ1,θ2,θ3,…,θNT. Since Pe00 is a function of joint angle Θ, the differential operation of Pe00 to Θ is as follows:(14)dpx0dpy0dpz0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNdθ1dθ1⋮dθN,where N is the number of joints. According to assumption (iv), equation (14) can be written as(15)Δxa0Δya0Δza0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNΔθ1Δθ1⋮ΔθN,where Δxa0, Δya0, and Δza0 represent the three components of positioning error of the EE due to the rotary deformation around joint’s axis in {0}. Combined with equation (13), equation (15) can be further written as(16)ΔPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNna110⋯00na22⋯0⋮⋮⋱⋮00⋯naNNCa1Ca1⋮CaN.Equation (16) preliminarily indicates the relationship between the torque applied to the joint and the positioning error. However, equation (16) still cannot describe this relationship sufficiently. The main reason is that the direction of Δθi (positive or negative) is not associated with the subjected torque in Δθi=Cainaii. Therefore, Δθi=Cainaii can be rewritten as (17) by introducing a sign function as(17)Δθi=Cainaiisign−niiz,where niiz∈ℝ stands for the component z of nii. In equation (17), a negative sign is added before niiz because the driving torque of joint and the subjected torque are a pair of balance torques. Hence, on the basis of equations (15) and (17), a complete positioning error model of axial deformation around the joint is given as follows:(18)ΔPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θN×n11sign−n11z0⋯00na22sign−n22z⋯0⋮⋮⋱⋮00⋯naNNsign−nNNzCa1Ca1⋮CaN.In particular, equation (18) indicates that the positioning error ΔPa0 of the robots due to the rotary deformation around the joint is a function of the variables EE position Pe00 (or joint angle Θ=θ1,θ2,θ3,…,θNT), end load F(or na11na22,……naNN), and the gravity of the robot. Moreover, even though the end load F does not change, the positioning error ΔPa0 is also different when the robot is at various configurations.


### 3.2. Positioning Error Modelling of Radial Deformation around the Joint

#### 3.2.1. Rotation Transformation around Arbitrary Axis

For the integrity of the modelling process, the concept of rotation transformation around the arbitrary axis will be introduced in this section. It is assumed that the vectorKA=kx,ky,kzT is an identity vector in {A}. According to the right-hand rule, the rotation transformation matrix of rotating θ around axis of KA is as follows:(19)RKθ=kxkxvθ+cθkxkyvθ−kzsθkxkzvθ+kysθkxkyvθ+kzsθkykyvθ+cθkykzvθ−kxsθkxkzvθ−kysθkykzvθ+kxsθkzkzvθ+cθ=cθ⋅I3×3+vθ⋅KAKTA+sθ⋅KA.Equation (19) is also called Rodigues’s formula, but it is not linearized form concerning θ. In the light of assumption (iv), when the rotary angle Δθ is small enough, equation (19) can be further simplified as(20)RKΔθ=100010001+0−kzkykz0−kx−kykx0Δθ=I3×3+SKAΔθ.Equation (20) indicates that the rotation transformation matrix can be handled via a linearized function with respect to the variable Δθ after using assumption (iv), which is beneficial to the linearization of the positioning error model.


#### 3.2.2. Modelling of Radial Deformation

As mentioned above, it is assumed that the end point of the robot is still at pointPe0. With the influence of end loads and gravity, each joint will have a slight rotary deformation Δγi around the radial direction of nii. Thus, the coordinate of EE position can be described by (21) with the rotary deformation of each joint:(21)P1e1=Rn11Δγ1P1e0,P2e2=Rn22Δγ2P2e1,P3e3=Rn33Δγ3P3e2,⋮Piei=RniiΔγiPiei−1,⋮PNeN=RnNNΔγNPNeN−1,where N is the number of joint. Peii stands for the EE position with the rotary deformation Δγi at joint i in {i}. It should be noted that in (21), the coordinates of all points are described by homogeneous coordinates and the rotation transformation is also the homogeneous form in order to be convenient for the following coordinate transformation. Then, the variation of each joint before and after rotary deformation can be obtained based on equation (21) as follows:(22)Δ1Pr1=P1e1−P1e0,Δ2Pr2=P2e2−P2e1,Δ3Pr3=P3e3−P3e2,⋮ΔiPri=Piei−Piei−1,⋮ΔNPrN=PNeN−PNeN−1,where ΔPrii denotes the variation for the EE position at joint i before and after rotary deformation in {i}. Next, the analysis about ΔPrii will be performed. First, ΔPri0 is obtained by transforming ΔPrii into {0}:(23)ΔPri0=Ti0ΔPrii.Combining with equations (20), (21), (22), and (23), the following equation can be derived:(24)ΔPri0=Ti0S˜−n^riiPei−1i⋅Δγi,where n^rii represents an identity vector corresponding to nrii. It can be found by observing equation (21) that Pei−1i is deduced from the term of Pe0i in equation (24). Thus, Pei−1i can also be expanded to Pe0i. According to (21) and (22), the following equations are obtained:(25)Piei−1=Ti−1iPi−1ei−1,Pi−1ei−1=Δi−1Pri−1+Pi−1ei−2,Δi−1Pri−1=S~−n^i−1ri−1Pi−1ei−2⋅Δγi−1.Next, the following equation is obtained according to equation (25):(26)Pei−1i=Ti−1iS˜−n^ri−1i−1Pei−2i−1⋅Δγi−1+Pei−2i−1.Equation (26) gives a recursive relationship between EE position Pei−1i before the rotary deformation at joint i and the EE position Pei−2i−1 before the rotary deformation at joint i-1. Based on (26), equation of Pei−1i including Pe0−1i is written as follows:(27)Pei−1i=∑k=1i−1TkiS˜−n^rkkPek−1k⋅Δγk+T1iPe01.AlthoughPei−1i has been expanded to Pe0i, the above equation still contains Pe12, Pe23,..., Pei−2i−1. Substituting equation (27) into equations (24) and (28), we have(28)ΔPri0=Ti0S˜−n^rii∑k=1i−1TkiS˜−n^nkkPek−1k⋅ΔγkΔγi+T1iPe01⋅Δγi.In (28), since Δγ satisfies assumption (iv) and ΔγkΔγi is high-order terms of Δγ, then ΔγkΔγi can be removed. When ΔγkΔγi=0, equation (28) will be further simplified as(29)ΔPri0=Ti0S˜−n^riiPe0i⋅Δγi,where the last component of ΔPri0 is zero since the homogeneous coordinates are employed, i.e., ΔPri0=Δxr0,Δyr0,Δzr0,0T. Accordingly, a transformation for equation (29) should be performed to eliminate the last component of ΔPri0 that is zero. The specific transforming process is as follows:(30)Δ0Pri=Ti0S˜−n^iriPie0⋅Δγi,Δxr0Δyr0Δzr00=Ri0P0iorg01S−n^iri000Pie03×11⋅Δγi=Ri0S−n^iriPie03×10⋅Δγi,Δxr0Δyr0Δzr0=Ri0S−n^iriPie03×1⋅Δγi,P0ri3×1=Ri0S−n^iriPie03×1⋅Δγi,which shows that the EE position of robots Pe0i and the torque of the joint can be calculated by the rotary deformation Δγi, and the positioning error Pri0 caused by the rotary deformation the around radial direction of the joint can also be calculated. Then, the total positioning error due to the rotary deformation Δγi is determined as(31)Pr3×10=∑i=1NPr3×10.From equations (13), (30), and (31), the positioning error model after the rotary deformation Δγi is derived as follows:(32)Pr3×10=∑i=1NRi0S−n^riiPe0i3×1⋅Crinrii,when Yi=nriiRi0Sn^riiPe0i3×1, equation (32) will be written as(33)Pr3×10=Y1Y2⋯YNCr1Cr2⋮CrN.Similar to equation (18), equation (33) illustrates that when the force Fe is applied to the EE of the robot and its gravity is considered, the radial torque at the joint will lead to the positioning error ΔPr3×10.


### 3.3. Error Modelling including Axial and Radial Deformation

According to the analysis in Sections3.1 and 3.2, two positioning error models that correspond to the axial and radial deformation at the joint, respectively, have been acquired. Since the positioning errors ΔPa0 and ΔPr0 are in {0}, they can be composited to obtain a complete positioning error model which includes both the influence of axial and radial deformation, as shown in the following equation:(34)Δ0P=ΔPa0+ΔPr0.Combining with equations (18) and (33), the complete positioning error model is derived as follows:(35)Δ0P=JPa0ΛaJ˜Pr0Λr3×2NCaCr2N×1,where(36)JPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θN,J˜Pr0=col1col2⋯coli⋯colN,with arbitrary column(37)coli=Ri0S−n^riiPe0i3×1,Ca=Ca1,Ca2,Ca3,…,CaNT,Cr=Cr1,Cr2,Cr3,…,CrNT,Λa=na11sign−n11z0⋯00na22sign−n22z⋯0⋮⋮⋱⋮00⋯naNNsign−nNNz,Λr=nr110⋯00nr12⋯0⋮⋮⋱⋮00⋯nrNN.Equation (35) describes the effect of the end load and the gravity on the positioning error at an arbitrary position in the working space of the robot, and it can be rewritten in a simplified form as follows:(38)Δ0P=ΦF,ΘC,where ΦF,Θ is a simplified expression of JPa0ΛaJ˜Pr0Λr3×2N, and C means Ca,CrT.Remark 2.
Two conclusions can be obtained according to (35). (i) The mapping relationship between the positioning error of the robot and its loads including the end load and gravity is linear. (ii) The positioning error is affected by the load as well as the poses of robots. It should be noted that Ca and Cr stand for the overall stiffness of all components that make up the joint in (35). Hence, they cannot be used to measure the stiffness of a specific part of joints. Moreover, Ca and Cr describe the torsional stiffness and bending stiffness of joints, respectively. This is consistent with the practical situation in which the torsional stiffness is different from the bending stiffness. Finally, since Ca and Cr are the integrated stiffness, they cannot be obtained directly by measurement. However, many mathematical methods can be used to identify Ca and Cr, e.g., least square method [29, 30], genetic algorithm [31], particle swarm optimization algorithm [32], Kalman filtering algorithm [33], etc.


## 3.1. Positioning Error Model of Axial Deformation around the Joint

When the EE of the robot is at a pointPe0 in Cartesian coordinate system and its coordinates are Pe00=px0,py0,pz0T in {0}. The joint angle corresponding to Pe00 is Θ=θ1,θ2,θ3,…,θNT. Since Pe00 is a function of joint angle Θ, the differential operation of Pe00 to Θ is as follows:(14)dpx0dpy0dpz0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNdθ1dθ1⋮dθN,where N is the number of joints. According to assumption (iv), equation (14) can be written as(15)Δxa0Δya0Δza0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNΔθ1Δθ1⋮ΔθN,where Δxa0, Δya0, and Δza0 represent the three components of positioning error of the EE due to the rotary deformation around joint’s axis in {0}. Combined with equation (13), equation (15) can be further written as(16)ΔPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θNna110⋯00na22⋯0⋮⋮⋱⋮00⋯naNNCa1Ca1⋮CaN.Equation (16) preliminarily indicates the relationship between the torque applied to the joint and the positioning error. However, equation (16) still cannot describe this relationship sufficiently. The main reason is that the direction of Δθi (positive or negative) is not associated with the subjected torque in Δθi=Cainaii. Therefore, Δθi=Cainaii can be rewritten as (17) by introducing a sign function as(17)Δθi=Cainaiisign−niiz,where niiz∈ℝ stands for the component z of nii. In equation (17), a negative sign is added before niiz because the driving torque of joint and the subjected torque are a pair of balance torques. Hence, on the basis of equations (15) and (17), a complete positioning error model of axial deformation around the joint is given as follows:(18)ΔPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θN×n11sign−n11z0⋯00na22sign−n22z⋯0⋮⋮⋱⋮00⋯naNNsign−nNNzCa1Ca1⋮CaN.In particular, equation (18) indicates that the positioning error ΔPa0 of the robots due to the rotary deformation around the joint is a function of the variables EE position Pe00 (or joint angle Θ=θ1,θ2,θ3,…,θNT), end load F(or na11na22,……naNN), and the gravity of the robot. Moreover, even though the end load F does not change, the positioning error ΔPa0 is also different when the robot is at various configurations.


## 3.2. Positioning Error Modelling of Radial Deformation around the Joint

### 3.2.1. Rotation Transformation around Arbitrary Axis

For the integrity of the modelling process, the concept of rotation transformation around the arbitrary axis will be introduced in this section. It is assumed that the vectorKA=kx,ky,kzT is an identity vector in {A}. According to the right-hand rule, the rotation transformation matrix of rotating θ around axis of KA is as follows:(19)RKθ=kxkxvθ+cθkxkyvθ−kzsθkxkzvθ+kysθkxkyvθ+kzsθkykyvθ+cθkykzvθ−kxsθkxkzvθ−kysθkykzvθ+kxsθkzkzvθ+cθ=cθ⋅I3×3+vθ⋅KAKTA+sθ⋅KA.Equation (19) is also called Rodigues’s formula, but it is not linearized form concerning θ. In the light of assumption (iv), when the rotary angle Δθ is small enough, equation (19) can be further simplified as(20)RKΔθ=100010001+0−kzkykz0−kx−kykx0Δθ=I3×3+SKAΔθ.Equation (20) indicates that the rotation transformation matrix can be handled via a linearized function with respect to the variable Δθ after using assumption (iv), which is beneficial to the linearization of the positioning error model.


### 3.2.2. Modelling of Radial Deformation

As mentioned above, it is assumed that the end point of the robot is still at pointPe0. With the influence of end loads and gravity, each joint will have a slight rotary deformation Δγi around the radial direction of nii. Thus, the coordinate of EE position can be described by (21) with the rotary deformation of each joint:(21)P1e1=Rn11Δγ1P1e0,P2e2=Rn22Δγ2P2e1,P3e3=Rn33Δγ3P3e2,⋮Piei=RniiΔγiPiei−1,⋮PNeN=RnNNΔγNPNeN−1,where N is the number of joint. Peii stands for the EE position with the rotary deformation Δγi at joint i in {i}. It should be noted that in (21), the coordinates of all points are described by homogeneous coordinates and the rotation transformation is also the homogeneous form in order to be convenient for the following coordinate transformation. Then, the variation of each joint before and after rotary deformation can be obtained based on equation (21) as follows:(22)Δ1Pr1=P1e1−P1e0,Δ2Pr2=P2e2−P2e1,Δ3Pr3=P3e3−P3e2,⋮ΔiPri=Piei−Piei−1,⋮ΔNPrN=PNeN−PNeN−1,where ΔPrii denotes the variation for the EE position at joint i before and after rotary deformation in {i}. Next, the analysis about ΔPrii will be performed. First, ΔPri0 is obtained by transforming ΔPrii into {0}:(23)ΔPri0=Ti0ΔPrii.Combining with equations (20), (21), (22), and (23), the following equation can be derived:(24)ΔPri0=Ti0S˜−n^riiPei−1i⋅Δγi,where n^rii represents an identity vector corresponding to nrii. It can be found by observing equation (21) that Pei−1i is deduced from the term of Pe0i in equation (24). Thus, Pei−1i can also be expanded to Pe0i. According to (21) and (22), the following equations are obtained:(25)Piei−1=Ti−1iPi−1ei−1,Pi−1ei−1=Δi−1Pri−1+Pi−1ei−2,Δi−1Pri−1=S~−n^i−1ri−1Pi−1ei−2⋅Δγi−1.Next, the following equation is obtained according to equation (25):(26)Pei−1i=Ti−1iS˜−n^ri−1i−1Pei−2i−1⋅Δγi−1+Pei−2i−1.Equation (26) gives a recursive relationship between EE position Pei−1i before the rotary deformation at joint i and the EE position Pei−2i−1 before the rotary deformation at joint i-1. Based on (26), equation of Pei−1i including Pe0−1i is written as follows:(27)Pei−1i=∑k=1i−1TkiS˜−n^rkkPek−1k⋅Δγk+T1iPe01.AlthoughPei−1i has been expanded to Pe0i, the above equation still contains Pe12, Pe23,..., Pei−2i−1. Substituting equation (27) into equations (24) and (28), we have(28)ΔPri0=Ti0S˜−n^rii∑k=1i−1TkiS˜−n^nkkPek−1k⋅ΔγkΔγi+T1iPe01⋅Δγi.In (28), since Δγ satisfies assumption (iv) and ΔγkΔγi is high-order terms of Δγ, then ΔγkΔγi can be removed. When ΔγkΔγi=0, equation (28) will be further simplified as(29)ΔPri0=Ti0S˜−n^riiPe0i⋅Δγi,where the last component of ΔPri0 is zero since the homogeneous coordinates are employed, i.e., ΔPri0=Δxr0,Δyr0,Δzr0,0T. Accordingly, a transformation for equation (29) should be performed to eliminate the last component of ΔPri0 that is zero. The specific transforming process is as follows:(30)Δ0Pri=Ti0S˜−n^iriPie0⋅Δγi,Δxr0Δyr0Δzr00=Ri0P0iorg01S−n^iri000Pie03×11⋅Δγi=Ri0S−n^iriPie03×10⋅Δγi,Δxr0Δyr0Δzr0=Ri0S−n^iriPie03×1⋅Δγi,P0ri3×1=Ri0S−n^iriPie03×1⋅Δγi,which shows that the EE position of robots Pe0i and the torque of the joint can be calculated by the rotary deformation Δγi, and the positioning error Pri0 caused by the rotary deformation the around radial direction of the joint can also be calculated. Then, the total positioning error due to the rotary deformation Δγi is determined as(31)Pr3×10=∑i=1NPr3×10.From equations (13), (30), and (31), the positioning error model after the rotary deformation Δγi is derived as follows:(32)Pr3×10=∑i=1NRi0S−n^riiPe0i3×1⋅Crinrii,when Yi=nriiRi0Sn^riiPe0i3×1, equation (32) will be written as(33)Pr3×10=Y1Y2⋯YNCr1Cr2⋮CrN.Similar to equation (18), equation (33) illustrates that when the force Fe is applied to the EE of the robot and its gravity is considered, the radial torque at the joint will lead to the positioning error ΔPr3×10.


## 3.2.1. Rotation Transformation around Arbitrary Axis

For the integrity of the modelling process, the concept of rotation transformation around the arbitrary axis will be introduced in this section. It is assumed that the vectorKA=kx,ky,kzT is an identity vector in {A}. According to the right-hand rule, the rotation transformation matrix of rotating θ around axis of KA is as follows:(19)RKθ=kxkxvθ+cθkxkyvθ−kzsθkxkzvθ+kysθkxkyvθ+kzsθkykyvθ+cθkykzvθ−kxsθkxkzvθ−kysθkykzvθ+kxsθkzkzvθ+cθ=cθ⋅I3×3+vθ⋅KAKTA+sθ⋅KA.Equation (19) is also called Rodigues’s formula, but it is not linearized form concerning θ. In the light of assumption (iv), when the rotary angle Δθ is small enough, equation (19) can be further simplified as(20)RKΔθ=100010001+0−kzkykz0−kx−kykx0Δθ=I3×3+SKAΔθ.Equation (20) indicates that the rotation transformation matrix can be handled via a linearized function with respect to the variable Δθ after using assumption (iv), which is beneficial to the linearization of the positioning error model.


## 3.2.2. Modelling of Radial Deformation

As mentioned above, it is assumed that the end point of the robot is still at pointPe0. With the influence of end loads and gravity, each joint will have a slight rotary deformation Δγi around the radial direction of nii. Thus, the coordinate of EE position can be described by (21) with the rotary deformation of each joint:(21)P1e1=Rn11Δγ1P1e0,P2e2=Rn22Δγ2P2e1,P3e3=Rn33Δγ3P3e2,⋮Piei=RniiΔγiPiei−1,⋮PNeN=RnNNΔγNPNeN−1,where N is the number of joint. Peii stands for the EE position with the rotary deformation Δγi at joint i in {i}. It should be noted that in (21), the coordinates of all points are described by homogeneous coordinates and the rotation transformation is also the homogeneous form in order to be convenient for the following coordinate transformation. Then, the variation of each joint before and after rotary deformation can be obtained based on equation (21) as follows:(22)Δ1Pr1=P1e1−P1e0,Δ2Pr2=P2e2−P2e1,Δ3Pr3=P3e3−P3e2,⋮ΔiPri=Piei−Piei−1,⋮ΔNPrN=PNeN−PNeN−1,where ΔPrii denotes the variation for the EE position at joint i before and after rotary deformation in {i}. Next, the analysis about ΔPrii will be performed. First, ΔPri0 is obtained by transforming ΔPrii into {0}:(23)ΔPri0=Ti0ΔPrii.Combining with equations (20), (21), (22), and (23), the following equation can be derived:(24)ΔPri0=Ti0S˜−n^riiPei−1i⋅Δγi,where n^rii represents an identity vector corresponding to nrii. It can be found by observing equation (21) that Pei−1i is deduced from the term of Pe0i in equation (24). Thus, Pei−1i can also be expanded to Pe0i. According to (21) and (22), the following equations are obtained:(25)Piei−1=Ti−1iPi−1ei−1,Pi−1ei−1=Δi−1Pri−1+Pi−1ei−2,Δi−1Pri−1=S~−n^i−1ri−1Pi−1ei−2⋅Δγi−1.Next, the following equation is obtained according to equation (25):(26)Pei−1i=Ti−1iS˜−n^ri−1i−1Pei−2i−1⋅Δγi−1+Pei−2i−1.Equation (26) gives a recursive relationship between EE position Pei−1i before the rotary deformation at joint i and the EE position Pei−2i−1 before the rotary deformation at joint i-1. Based on (26), equation of Pei−1i including Pe0−1i is written as follows:(27)Pei−1i=∑k=1i−1TkiS˜−n^rkkPek−1k⋅Δγk+T1iPe01.AlthoughPei−1i has been expanded to Pe0i, the above equation still contains Pe12, Pe23,..., Pei−2i−1. Substituting equation (27) into equations (24) and (28), we have(28)ΔPri0=Ti0S˜−n^rii∑k=1i−1TkiS˜−n^nkkPek−1k⋅ΔγkΔγi+T1iPe01⋅Δγi.In (28), since Δγ satisfies assumption (iv) and ΔγkΔγi is high-order terms of Δγ, then ΔγkΔγi can be removed. When ΔγkΔγi=0, equation (28) will be further simplified as(29)ΔPri0=Ti0S˜−n^riiPe0i⋅Δγi,where the last component of ΔPri0 is zero since the homogeneous coordinates are employed, i.e., ΔPri0=Δxr0,Δyr0,Δzr0,0T. Accordingly, a transformation for equation (29) should be performed to eliminate the last component of ΔPri0 that is zero. The specific transforming process is as follows:(30)Δ0Pri=Ti0S˜−n^iriPie0⋅Δγi,Δxr0Δyr0Δzr00=Ri0P0iorg01S−n^iri000Pie03×11⋅Δγi=Ri0S−n^iriPie03×10⋅Δγi,Δxr0Δyr0Δzr0=Ri0S−n^iriPie03×1⋅Δγi,P0ri3×1=Ri0S−n^iriPie03×1⋅Δγi,which shows that the EE position of robots Pe0i and the torque of the joint can be calculated by the rotary deformation Δγi, and the positioning error Pri0 caused by the rotary deformation the around radial direction of the joint can also be calculated. Then, the total positioning error due to the rotary deformation Δγi is determined as(31)Pr3×10=∑i=1NPr3×10.From equations (13), (30), and (31), the positioning error model after the rotary deformation Δγi is derived as follows:(32)Pr3×10=∑i=1NRi0S−n^riiPe0i3×1⋅Crinrii,when Yi=nriiRi0Sn^riiPe0i3×1, equation (32) will be written as(33)Pr3×10=Y1Y2⋯YNCr1Cr2⋮CrN.Similar to equation (18), equation (33) illustrates that when the force Fe is applied to the EE of the robot and its gravity is considered, the radial torque at the joint will lead to the positioning error ΔPr3×10.


## 3.3. Error Modelling including Axial and Radial Deformation

According to the analysis in Sections3.1 and 3.2, two positioning error models that correspond to the axial and radial deformation at the joint, respectively, have been acquired. Since the positioning errors ΔPa0 and ΔPr0 are in {0}, they can be composited to obtain a complete positioning error model which includes both the influence of axial and radial deformation, as shown in the following equation:(34)Δ0P=ΔPa0+ΔPr0.Combining with equations (18) and (33), the complete positioning error model is derived as follows:(35)Δ0P=JPa0ΛaJ˜Pr0Λr3×2NCaCr2N×1,where(36)JPa0=∂px0∂θ1∂px0∂θ2⋯∂px0∂θN∂py0∂θ1∂py0∂θ2⋯∂py0∂θN∂pz0∂θ1∂pz0∂θ2⋯∂pz0∂θN,J˜Pr0=col1col2⋯coli⋯colN,with arbitrary column(37)coli=Ri0S−n^riiPe0i3×1,Ca=Ca1,Ca2,Ca3,…,CaNT,Cr=Cr1,Cr2,Cr3,…,CrNT,Λa=na11sign−n11z0⋯00na22sign−n22z⋯0⋮⋮⋱⋮00⋯naNNsign−nNNz,Λr=nr110⋯00nr12⋯0⋮⋮⋱⋮00⋯nrNN.Equation (35) describes the effect of the end load and the gravity on the positioning error at an arbitrary position in the working space of the robot, and it can be rewritten in a simplified form as follows:(38)Δ0P=ΦF,ΘC,where ΦF,Θ is a simplified expression of JPa0ΛaJ˜Pr0Λr3×2N, and C means Ca,CrT.Remark 2.
Two conclusions can be obtained according to (35). (i) The mapping relationship between the positioning error of the robot and its loads including the end load and gravity is linear. (ii) The positioning error is affected by the load as well as the poses of robots. It should be noted that Ca and Cr stand for the overall stiffness of all components that make up the joint in (35). Hence, they cannot be used to measure the stiffness of a specific part of joints. Moreover, Ca and Cr describe the torsional stiffness and bending stiffness of joints, respectively. This is consistent with the practical situation in which the torsional stiffness is different from the bending stiffness. Finally, since Ca and Cr are the integrated stiffness, they cannot be obtained directly by measurement. However, many mathematical methods can be used to identify Ca and Cr, e.g., least square method [29, 30], genetic algorithm [31], particle swarm optimization algorithm [32], Kalman filtering algorithm [33], etc.


## 4. Simulations

The proposed positioning error model will be verified by simulations. The procedure of the verification consists of six steps as shown in Figure3. It is noted that if dynamic parameters of the robot are known, step S3 can be omitted. In this paper, dynamic parameters are obtained by Computer-Aided Design (CAD) method, so step S3 is represented by a dotted box here.Figure 3
The flowchart of verification for the proposed model.


### 4.1. Kinematics Modelling

A 6-DOF general-purpose robot is used to verify the effectiveness and generality of the proposed method. First, the coordinate systems of the robot are established according to the D-H method, as shown in Figure4. To describe the EE position of the robot, the origin of {6} is set at its EE. Then, kinematics parameters of the robot are obtained, as shown in Table 1 based on the coordinate systems in Figure 4.Figure 4
Coordinate systems of the robot.Table 1
Kinematics parameters of the robot.
iαi−1 (rad)ai−1 (m)Di−1 (m)θi−1 (rad)1000.504θ12π/20.1700.101θ2 + π/2300.70−0.120θ34π/20.1400.760θ45−π/200θ5 + π6−π/200.170θ6According to the kinematics parameters in Table1, the kinematic model of the robot is established by equation (39). Moreover, the homogeneous transformation matrix Tji and the rotation transformation matrix Rji can be obtained, which will be used in the positioning error model:(39)Tii−1=cθi−sθi0ai−1sθicαi−1cθicαi−1−sαi−1−sαi−1disθisαi−1cθisαi−1cαi−1cαi−1di0001,Ti0=∏k=1iTkk−1,Tji=T0iTj0.


### 4.2. Dynamics Parameters

To acquire the driving torque,nii and dynamics parameters of the robot are established on the basis of Newton–Euler method, as shown in equations (3), (4), and (5). Generally, the dynamics parameters are obtained by identification [34] or calculation from the design parameters. The dynamics parameters of the robot can be acquired by Computer-Aided Design (CAD) method based on the kinematics model, as shown in Table 2.Table 2
The mass and centroid position of each link.
imi/kgrxii (m)ryii (m)rzii (m)131.1380.126−0.076−0.185252.1320.3900−0.025348.8130.124−0.0760.005435.08800−0.19056.67500.183061.89300−0.014


### 4.3. Simulation and Analysis

Combining with the above kinematics model and dynamics parameters, a simulation environment is constructed by using finite element method. The data set from finite element simulation is defined as practical values, which are used to identify the unknown parametersCa and Cr. Then, the identified parameters are put into the positioning error model. Finally, the predicted values according to the proposed method are compared with the results of finite element simulation.To estimate unknown parametersCa and Cr, a group of joint angles are chosen arbitrarily in the working space of the robot, i.e., Θ=44∘,−45∘,20∘,45∘,−30∘,80∘T. It should be noted that Θ needs to be converted into a radian system when calculating. A group of end loads F are used as shown in Table 3. In the light of Θ and F, the regression matrix ΦF,Θ is calculated. The positioning error Δ0P can be acquired from the finite element model as shown in Figure 5. There are 10 groups of data in Table 3.Table 3
Settings of the end load for parameters identification.
F/N0Δx (mm)0Δy (mm)0Δz (mm)(0, 0, −500)T5.6939 × 10−25.0193 × 10−2−3.7822 × 10−1(0, 0, −450)T5.5263 × 10−24.8607 × 10−2−3.6290 × 10−1(0, 0, −330)T5.1240 × 10−24.4801 × 10−2−3.2616 × 10−1(0, 0, −271)T4.9264 × 10−24.2929 × 10−2−3.0810 × 10−1(0, 0, −161)T4.5578 × 10−23.9440 × 10−2−2.7442 × 10−1(0, 0, −88)T4.3131 × 10−23.7124 × 10−2−2.5206 × 10−1(0, 0, −65)T4.2362 × 10−23.6395 × 10−2−2.4503 × 10−1(0, 0, −30)T4.1187 × 10−23.5285 × 10−2−2.3431 × 10−1(0, 0, −10)T4.0519 × 10−23.4650 × 10−2−2.2818 × 10−1(0, 0, 0)T4.0184 × 10−23.4333 × 10−2−2.2512 × 10−1Figure 5
The finite element model.Since (35) is linear with respect to parameters Ca and Cr, the least square method is used to estimate the unknown parameters Ca and Cr. The results of parameter estimation are as shown in Table 4. It can be found that the values of C^a1, C^a6, and C^r6 are zero. Nevertheless, this does not mean that the real stiffness coefficient at joints 1 and 6 is zero but means that their changes have no effects on the positioning error. In Table 4, except the case where estimated parameter C is zero, it can be also found that some identified values are negative. According to Δφ=Cn in assumption (v), when flexibility coefficient C is positive, it indicates that the direction of rotary deformation Δφ and joint torque n are the same. And when C is negative, it indicates that the direction of Δφ is opposite to the direction of n. Meanwhile, it also shows that the stiffness parameters Ca and Cr do not possess practical physical significance but merely mathematical meaning in the proposed model.Table 4
The results of parameters identification.
ParametersValue (rad/N•m)ParametersValue (rad/N•m)C^a10C^r19.9312 × 10−5C^a23.9922 × 10−5C^r2−6.7172 × 10−5C^a3−5.0087 × 10−5C^r3−2.4600 × 10−4C^a42.5951 × 10−5C^r41.3356 × 10−3C^a53.9436 × 10−4C^r51.0381 × 10−4C^a60C^r60Table 5
The data for model verification.
F/NThe measured valuesThe predicted values0Δx (mm)0Δy (mm)0Δz (mm)0Δx (mm)0Δy (mm)0Δz (mm)(0, 0, −78.80)T4.2823 × 10−23.6833 × 10−2−2.4925 × 10−14.2824 × 10−23.6832 × 10−2−2.4925 × 10−1(0, 0, −485.30)T5.6446 × 10−24.9727 × 10−2−3.7371 × 10−15.6445 × 10−24.9727 × 10−2−3.7372 × 10−1(0, 0, −478.58)T5.6221 × 10−24.9513 × 10−2−3.7166 × 10−15.6220 × 10−24.9514 × 10−2−3.7165 × 10−1(0, 0, −242.69)T4.8315 × 10−24.2031 × 10−2−2.9943 × 10−14.8316 × 10−24.2033 × 10−2−2.9942 × 10−1(0, 0, −400.14)T5.3589 × 10−24.7025 × 10−2−3.4763 × 10−15.3592 × 10−24.7026 × 10−2−3.4764 × 10−1(0, 0, −70.94)T4.2561 × 10−23.6583 × 10−2−2.4684 × 10−14.2560 × 10−23.6582 × 10−2−2.4686 × 10−1(0, 0, −210.88)T4.7250 × 10−24.1021 × 10−2−2.8969 × 10−14.7249 × 10−24.1022 × 10−2−2.8970 × 10−1(0, 0, −457.86)T5.5526 × 10−24.8856 × 10−2−3.6531 × 10−15.5524 × 10−24.8857 × 10−2−3.6532 × 10−1(0, 0, −396.10)T5.3454 × 10−24.6897 × 10−2−3.4640 × 10−15.3456 × 10−24.6897 × 10−2−3.4643 × 10−1(0, 0, −497.74)T5.6260 × 10−24.2823 × 10−2−3.7201 × 10−15.6261 × 10−24.2822 × 10−2−3.7202 × 10−1(0, 0, −327.87)T5.1170 × 10−24.9550 × 10−2−3.2551 × 10−15.1169 × 10−24.9551 × 10−2−3.2552 × 10−1(0, 0, −17.00)T4.0753 × 10−24.4733 × 10−2−2.3033 × 10−14.0752 × 10−24.4734 × 10−2−2.3032 × 10−1(0, 0, −424.56)T5.4409 × 10−23.4872 × 10−2−3.5511 × 10−15.4410 × 10−23.4873 × 10−2−3.5514 × 10−1(0, 0, −466.99)T5.5833 × 10−24.7800 × 10−2−3.6811 × 10−15.5834 × 10−24.7801 × 10−2−3.6813 × 10−1(0, 0, −339.36)T5.1554 × 10−24.9145 × 10−2−3.2903 × 10−15.1555 × 10−24.9144 × 10−2−3.2901 × 10−1(0, 0, −378.87)T5.2878 × 10−24.5097 × 10−2−3.4112 × 10−15.2876 × 10−24.5099 × 10−2−3.4113 × 10−1(0, 0, −371.56)T5.2634 × 10−24.6351 × 10−2−3.3889 × 10−15.2636 × 10−24.6353 × 10−2−3.3887 × 10−1(0, 0, −196.11)T4.6755 × 10−24.6119 × 10−2−2.8517 × 10−14.6754 × 10−24.6117 × 10−2−2.8519 × 10−1(0, 0, −327.73)T5.1164 × 10−24.0552 × 10−2−3.2546 × 10−15.1163 × 10−24.0554 × 10−2−3.2545 × 10−1(0, 0, −85.59)T4.3050 × 10−23.7048 × 10−2−2.5133 × 10−14.3051 × 10−23.7046 × 10−2−2.5132 × 10−1To measure the accuracy of the identified parameters, the index of relative error is introduced. A linear model can be expressed asAx=b and the relative error about the estimated value x^ can be defined as follows:(40)er=b−Ax^b×100%.According to (40), it can be found that limx^⟶xer=0. In other words, the value of er describes a degree of closeness between the estimated value x^ and the true value x. In addition, the value of er can also be used to measure the degree of closeness between the predicted value Ax^ and the true value b. According to (40), the relative error with respect to the identified parameters C^a and C^r is obtained as shown in Figure 6.Figure 6
The accuracy of identified parameters.It can be seen from Figure6 that the relative error between the theoretical values ΦF,ΘC^ and the measured values Δ0P is very small in the three directions, and it is between −0.0025% and 0.0025%.20 groups of different end loads are selected randomly to verify the effectiveness of the presented method, as shown in the first column of Table5. Then, the three components of the positioning error corresponding to each load can be acquired with finite element simulations as shown in the last three columns of Table 5. The predicted values ΦF,ΘC^ of the positioning error can be calculated based on these end loads and identified parameters C^a and C^r. The relative error between the predicted values ΦF,ΘC^ and the measured values Δ0P is calculated by (40). The results are as shown in Figure 7.Figure 7
The accuracy for predicted valueΔ0P.From Figure7, it can be found that the relative errors between the predicted values ΦF,ΘC^ and the measured values Δ0P are very small in all the three directions. They are all in the range of [−0.004%, 0.003%]. Compared with the relative error shown in Figure 6, the relative error in Figure 7 is larger. The main reason is that the former group of data is involved in the parameter identification, but the latter is only used to predict the positioning error of the robot. With a predicted accuracy of 99.996%, the accuracy of the model is quite high in predicting and compensating positioning error Δ0P. In practice, the predicted accuracy may reduce when the real data used to identify the parameters contains the noise.


## 4.1. Kinematics Modelling

A 6-DOF general-purpose robot is used to verify the effectiveness and generality of the proposed method. First, the coordinate systems of the robot are established according to the D-H method, as shown in Figure4. To describe the EE position of the robot, the origin of {6} is set at its EE. Then, kinematics parameters of the robot are obtained, as shown in Table 1 based on the coordinate systems in Figure 4.Figure 4
Coordinate systems of the robot.Table 1
Kinematics parameters of the robot.
iαi−1 (rad)ai−1 (m)Di−1 (m)θi−1 (rad)1000.504θ12π/20.1700.101θ2 + π/2300.70−0.120θ34π/20.1400.760θ45−π/200θ5 + π6−π/200.170θ6According to the kinematics parameters in Table1, the kinematic model of the robot is established by equation (39). Moreover, the homogeneous transformation matrix Tji and the rotation transformation matrix Rji can be obtained, which will be used in the positioning error model:(39)Tii−1=cθi−sθi0ai−1sθicαi−1cθicαi−1−sαi−1−sαi−1disθisαi−1cθisαi−1cαi−1cαi−1di0001,Ti0=∏k=1iTkk−1,Tji=T0iTj0.


## 4.2. Dynamics Parameters

To acquire the driving torque,nii and dynamics parameters of the robot are established on the basis of Newton–Euler method, as shown in equations (3), (4), and (5). Generally, the dynamics parameters are obtained by identification [34] or calculation from the design parameters. The dynamics parameters of the robot can be acquired by Computer-Aided Design (CAD) method based on the kinematics model, as shown in Table 2.Table 2
The mass and centroid position of each link.
imi/kgrxii (m)ryii (m)rzii (m)131.1380.126−0.076−0.185252.1320.3900−0.025348.8130.124−0.0760.005435.08800−0.19056.67500.183061.89300−0.014


## 4.3. Simulation and Analysis

Combining with the above kinematics model and dynamics parameters, a simulation environment is constructed by using finite element method. The data set from finite element simulation is defined as practical values, which are used to identify the unknown parametersCa and Cr. Then, the identified parameters are put into the positioning error model. Finally, the predicted values according to the proposed method are compared with the results of finite element simulation.To estimate unknown parametersCa and Cr, a group of joint angles are chosen arbitrarily in the working space of the robot, i.e., Θ=44∘,−45∘,20∘,45∘,−30∘,80∘T. It should be noted that Θ needs to be converted into a radian system when calculating. A group of end loads F are used as shown in Table 3. In the light of Θ and F, the regression matrix ΦF,Θ is calculated. The positioning error Δ0P can be acquired from the finite element model as shown in Figure 5. There are 10 groups of data in Table 3.Table 3
Settings of the end load for parameters identification.
F/N0Δx (mm)0Δy (mm)0Δz (mm)(0, 0, −500)T5.6939 × 10−25.0193 × 10−2−3.7822 × 10−1(0, 0, −450)T5.5263 × 10−24.8607 × 10−2−3.6290 × 10−1(0, 0, −330)T5.1240 × 10−24.4801 × 10−2−3.2616 × 10−1(0, 0, −271)T4.9264 × 10−24.2929 × 10−2−3.0810 × 10−1(0, 0, −161)T4.5578 × 10−23.9440 × 10−2−2.7442 × 10−1(0, 0, −88)T4.3131 × 10−23.7124 × 10−2−2.5206 × 10−1(0, 0, −65)T4.2362 × 10−23.6395 × 10−2−2.4503 × 10−1(0, 0, −30)T4.1187 × 10−23.5285 × 10−2−2.3431 × 10−1(0, 0, −10)T4.0519 × 10−23.4650 × 10−2−2.2818 × 10−1(0, 0, 0)T4.0184 × 10−23.4333 × 10−2−2.2512 × 10−1Figure 5
The finite element model.Since (35) is linear with respect to parameters Ca and Cr, the least square method is used to estimate the unknown parameters Ca and Cr. The results of parameter estimation are as shown in Table 4. It can be found that the values of C^a1, C^a6, and C^r6 are zero. Nevertheless, this does not mean that the real stiffness coefficient at joints 1 and 6 is zero but means that their changes have no effects on the positioning error. In Table 4, except the case where estimated parameter C is zero, it can be also found that some identified values are negative. According to Δφ=Cn in assumption (v), when flexibility coefficient C is positive, it indicates that the direction of rotary deformation Δφ and joint torque n are the same. And when C is negative, it indicates that the direction of Δφ is opposite to the direction of n. Meanwhile, it also shows that the stiffness parameters Ca and Cr do not possess practical physical significance but merely mathematical meaning in the proposed model.Table 4
The results of parameters identification.
ParametersValue (rad/N•m)ParametersValue (rad/N•m)C^a10C^r19.9312 × 10−5C^a23.9922 × 10−5C^r2−6.7172 × 10−5C^a3−5.0087 × 10−5C^r3−2.4600 × 10−4C^a42.5951 × 10−5C^r41.3356 × 10−3C^a53.9436 × 10−4C^r51.0381 × 10−4C^a60C^r60Table 5
The data for model verification.
F/NThe measured valuesThe predicted values0Δx (mm)0Δy (mm)0Δz (mm)0Δx (mm)0Δy (mm)0Δz (mm)(0, 0, −78.80)T4.2823 × 10−23.6833 × 10−2−2.4925 × 10−14.2824 × 10−23.6832 × 10−2−2.4925 × 10−1(0, 0, −485.30)T5.6446 × 10−24.9727 × 10−2−3.7371 × 10−15.6445 × 10−24.9727 × 10−2−3.7372 × 10−1(0, 0, −478.58)T5.6221 × 10−24.9513 × 10−2−3.7166 × 10−15.6220 × 10−24.9514 × 10−2−3.7165 × 10−1(0, 0, −242.69)T4.8315 × 10−24.2031 × 10−2−2.9943 × 10−14.8316 × 10−24.2033 × 10−2−2.9942 × 10−1(0, 0, −400.14)T5.3589 × 10−24.7025 × 10−2−3.4763 × 10−15.3592 × 10−24.7026 × 10−2−3.4764 × 10−1(0, 0, −70.94)T4.2561 × 10−23.6583 × 10−2−2.4684 × 10−14.2560 × 10−23.6582 × 10−2−2.4686 × 10−1(0, 0, −210.88)T4.7250 × 10−24.1021 × 10−2−2.8969 × 10−14.7249 × 10−24.1022 × 10−2−2.8970 × 10−1(0, 0, −457.86)T5.5526 × 10−24.8856 × 10−2−3.6531 × 10−15.5524 × 10−24.8857 × 10−2−3.6532 × 10−1(0, 0, −396.10)T5.3454 × 10−24.6897 × 10−2−3.4640 × 10−15.3456 × 10−24.6897 × 10−2−3.4643 × 10−1(0, 0, −497.74)T5.6260 × 10−24.2823 × 10−2−3.7201 × 10−15.6261 × 10−24.2822 × 10−2−3.7202 × 10−1(0, 0, −327.87)T5.1170 × 10−24.9550 × 10−2−3.2551 × 10−15.1169 × 10−24.9551 × 10−2−3.2552 × 10−1(0, 0, −17.00)T4.0753 × 10−24.4733 × 10−2−2.3033 × 10−14.0752 × 10−24.4734 × 10−2−2.3032 × 10−1(0, 0, −424.56)T5.4409 × 10−23.4872 × 10−2−3.5511 × 10−15.4410 × 10−23.4873 × 10−2−3.5514 × 10−1(0, 0, −466.99)T5.5833 × 10−24.7800 × 10−2−3.6811 × 10−15.5834 × 10−24.7801 × 10−2−3.6813 × 10−1(0, 0, −339.36)T5.1554 × 10−24.9145 × 10−2−3.2903 × 10−15.1555 × 10−24.9144 × 10−2−3.2901 × 10−1(0, 0, −378.87)T5.2878 × 10−24.5097 × 10−2−3.4112 × 10−15.2876 × 10−24.5099 × 10−2−3.4113 × 10−1(0, 0, −371.56)T5.2634 × 10−24.6351 × 10−2−3.3889 × 10−15.2636 × 10−24.6353 × 10−2−3.3887 × 10−1(0, 0, −196.11)T4.6755 × 10−24.6119 × 10−2−2.8517 × 10−14.6754 × 10−24.6117 × 10−2−2.8519 × 10−1(0, 0, −327.73)T5.1164 × 10−24.0552 × 10−2−3.2546 × 10−15.1163 × 10−24.0554 × 10−2−3.2545 × 10−1(0, 0, −85.59)T4.3050 × 10−23.7048 × 10−2−2.5133 × 10−14.3051 × 10−23.7046 × 10−2−2.5132 × 10−1To measure the accuracy of the identified parameters, the index of relative error is introduced. A linear model can be expressed asAx=b and the relative error about the estimated value x^ can be defined as follows:(40)er=b−Ax^b×100%.According to (40), it can be found that limx^⟶xer=0. In other words, the value of er describes a degree of closeness between the estimated value x^ and the true value x. In addition, the value of er can also be used to measure the degree of closeness between the predicted value Ax^ and the true value b. According to (40), the relative error with respect to the identified parameters C^a and C^r is obtained as shown in Figure 6.Figure 6
The accuracy of identified parameters.It can be seen from Figure6 that the relative error between the theoretical values ΦF,ΘC^ and the measured values Δ0P is very small in the three directions, and it is between −0.0025% and 0.0025%.20 groups of different end loads are selected randomly to verify the effectiveness of the presented method, as shown in the first column of Table5. Then, the three components of the positioning error corresponding to each load can be acquired with finite element simulations as shown in the last three columns of Table 5. The predicted values ΦF,ΘC^ of the positioning error can be calculated based on these end loads and identified parameters C^a and C^r. The relative error between the predicted values ΦF,ΘC^ and the measured values Δ0P is calculated by (40). The results are as shown in Figure 7.Figure 7
The accuracy for predicted valueΔ0P.From Figure7, it can be found that the relative errors between the predicted values ΦF,ΘC^ and the measured values Δ0P are very small in all the three directions. They are all in the range of [−0.004%, 0.003%]. Compared with the relative error shown in Figure 6, the relative error in Figure 7 is larger. The main reason is that the former group of data is involved in the parameter identification, but the latter is only used to predict the positioning error of the robot. With a predicted accuracy of 99.996%, the accuracy of the model is quite high in predicting and compensating positioning error Δ0P. In practice, the predicted accuracy may reduce when the real data used to identify the parameters contains the noise.


## 5. Conclusion

The main factors that affect the positioning accuracy of robots were analyzed considering the end loads and gravity. Based on the results of the analysis, it is found that the positioning error can be described by two parameters, i.e., the axial deformation and the radial deformation at the joint. A prediction and compensation model of positioning error was proposed based on the two kinds of deformations. The positioning error can be calculated according to the loads and gravity though the model forn-DOF industrial robots. Finite element simulation was used to verify the proposed model. The results of simulation showed that the proposed positioning error model can predict positioning errors. Future work will focus on the verification of the presented model by means of experiments and applying it to predict the positioning error under different loads to improve the positioning accuracy of industrial robots.



---

*Source: 8850751-2020-11-20.xml*