# Uncertainty Estimation of Robot Geometric Parameters and End-Effecter Position Based on New Generation GPS

Authors: Xiulan Wen; Shun He; GuiFang Qiao; Dongxia Wang; Aiguo Song; ChuanShuai Kang; Zhongyan Lv
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The robot end-effecter positioning accuracy can be improved by the calibration of robot geometric parameters errors. According to the requirements of new generation geometrical product specification (GPS), the calibration uncertainty should be given when the calibration results are given. In this paper, the modified Denavit-Hartenberg method (MDH) of six-joint series robot is established and the joint movement trajectory method is applied to calibrate the robot geometric parameters. The uncertainty contributors significant are analyzed and the calibration uncertainty of robot geometric parameters is estimated based on the guide to the expression of uncertainty in measurement (GUM). In order to overcome the limitations of GUM for highly nonlinear model and reduce computational cost based on Monte Carlo Simulation (MCS) error estimation, an adaptive MCS (AMCS) is proposed to estimate the uncertainty distribution of robot end-effector position. Simulation and practical example are illustrated and the experiments results confirm that not only can the proposed method evaluate the calibration uncertainty of geometric parameters, but also the uncertainty distribution of end-effecter positions in the whole robot workspace can be estimated by AMCS in which the number of MCS trials can be selected adaptively and the quality of the numerical results can be controlled directly. The proposed method not only can evaluate the uncertainty of six-joint series robot geometric parameters and end-effecter position rapidly and accurately, but also can be popularized to the estimation of calibration uncertainty of other kinds of robot geometric parameters.


---

## Body

## 1. Introduction

The main indexes of industrial robot positioning performance are repetitive positioning accuracy and absolute positioning accuracy [1]. Now most industrial manufacturers adopt the ISO 9283 norm [2], but the only upfront information regarding the positioning performance of an industrial robot continues to be a single measure specified as something like “positioning performance according to ISO 9283,” which actually refers to the average unidirectional position repeatability and accuracy at five poses obtained from thirty cycles [3, 4]. If the robot end-effector poses are manually taught, repetitive positioning accuracy is all that matters. The demand for industrial robots with higher absolute positioning accuracy has been growing continuously, especially in the fields of aerospace, medical treatment, and robot-assisted measurement, where positions are defined in a virtual space with respect to an absolute or relative coordinate system. Based on investigation of the error contribution from various sources, the robot geometric errors are responsible for about 90% of the total positioning error [5]. Therefore, many scholars have devoted to the study on geometric parameters calibration to improve absolute positioning accuracy [4–9]. The absolute accuracy of an ABB IRB1600 industrial robot was improved by using a 29-parameter calibration model which takes into account all possible geometry errors. The least squares optimization technique was developed to find the 29 error parameters that best fit the measures acquired with a laser tracker [4]. In order to improve the absolute positioning accuracy of an industrial robot with a parallelogram mechanism, Guo et al. proposed a multilevel calibration technique [6]. Joubair and Bonev developed a geometric calibration method based on sphere and distance constraints to improve the accuracy of a six-axis serial industrial robot in a specific target workspace [7]. Möller et al. presented an approach for increasing the absolute positioning accuracy of an industrial milling robot with help of a stereo camera system [8]. Nubiola presented the successful use of this measurement system for absolute robot calibration and the mean absolute positioning error was improved from 0.873 mm to 0.479 mm [9]. Wu et al. dealt with geometric calibration of industrial robots and focused on the reduction of the measurement noise impact by means of proper selection of the manipulator configurations in calibration experiments. The geometric calibration results for a KUKAKR-270 industrial robot showed that manipulator positioning accuracy was 5.5 times better compared to the noncalibrated robot [10].In terms of metrology based on new generation GPS, the measurement results must be given with a quantitative evaluation of their accuracy with measurement uncertainty [11]. Different approaches [12–17] have been proposed to quantify the uncertainty caused by measurement, modeling of complex physical systems or fuzzy systems. Among them, GUM [18, 19] is a highly practical and easily implement approach for cases where the measurement model is linear and the probability distribution of the output is normal distribution. It provides an analytical methodology for assessing and reporting the measurement uncertainty and can be applied to most fields of physical measurements. However, when the physical model is complex or the probability distribution of the output is significantly asymmetric, the GUM method to evaluate the measurement uncertainty will result in the inaccurate results. Supplement 1 to the GUM (GUM S1) is recommended by the Joint Committee for Guides in Metrology (JCGM) [20] which is implemented by Monte Carlo Simulation. ISO/TS15530 states that the most effective method for estimating uncertainty is computer simulation, and more specifically, a Monte Carlo simulation [21].Calibrations are measurements which determine quantity values. Because calibration devices, calibration standards, reagents, and tools are not perfect, environmental conditions, processes, procedures, and people are also imperfect and variable, each calibration has uncertainty associated with it [19]. According to new generation GPS, the calibration results must be given with a quantitative evaluation of their accuracy. Therefore, the study on the evaluation of calibration uncertainty of robot geometric parameters has drawn the attention of scholars and several researches have been reported [22–24]. Santolaria and Ginés presented Monte Carlo method to evaluate the uncertainty for a robot arm calibrated [22]. Li and Qu used GUM to evaluate the calibration uncertainty of geometric parameters [23]. Olarra and Axinte proposed a calibration method which enabled completely automatic identification of the kinematics of a walking hexapod robotic machine tool and proposed an analytical methodology to estimate the uncertainty of the identified geometric parameters [24].To estimate the calibration uncertainty of robot geometric parameters and end-effector position, the remainder of this paper is organized as follows. In Section2, the model of robot geometric parameters is established. The calibration procedure of geometric parameters and the calculation method of the end-effector position are given. In Section 3, The GUM method is used to calculate the calibration uncertainty of robot geometric parameters. AMCS for the uncertainty estimation of robot end-effector position is developed and the flowchart of AMCS is given. In Section 4, simulation and practical example are given. The uncertainty contributors significant are analyzed. And the calibration uncertainties of geometric parameters and end-effector position are estimated. Finally, conclusions are reached in Section 5.


## 2. Geometric Parameters Calibration and End-Effector Position Calculation

### 2.1. Geometric Parameters Calibration

#### 2.1.1. Robot Geometric Model

A number of different approaches including Denavit-Hartenberg (DH) method, S model method, Complete and Parametrically Continuous (CPC) model method, zero reference model method, and Product of Exponentials (POE) formula method exist for developing the geometric model of a robot manipulator [25]. The most popular model has been the DH model which is based on homogeneous transformation matrices [3]. This procedure consists of establishing coordinate systems on each joint axis. Each coordinate system is then related to the next through a specific set of coefficients in the homogeneous transformation matrices. To be suitable for robot calibration, the robot geometric model must satisfy three basic requirements: model completeness, parameter minimality, and model continuity. The conventional DH is not enough to achieve a complete model for geometric parameters calibration. In addition, when the adjacent joint axes are parallel, little changes at this state may cause dramatic changes at the common perpendicular between the two axes. The changes may result in robot singularity, which lead to discontinuity of the geometric model. Therefore, a modified DH (MDH) method was proposed, the core of which was to add an angle parameter to the DH model to describe the position relationship between adjacent parallel axes, to cancel the incompleteness problem [25–27]. In this paper, MDH model is adopted to establish the geometric parameters model of six-joint series robot for avoiding the singularity of conventional DH model.


#### 2.1.2. Acquisition of Joint Axis Space Position

For obtaining the spatial position of the robot joint axis, each axis is controlled to move separately and the target ball position is measured by the laser tracker. The spatial plane in which the trajectory is located is fitted.Supposing that the coordinates of measured points areP
j
(
x
j
,
y
j
,
z
j
)
,
j
=
1,2
,
⋯
n,n is the number of measured points. The space plane equation is expressed as(1)
A
x
+
B
y
+
C
z
+
1
=
0
X
=
A
B
C
T is the normal vector of the plane, which is also the direction vector n of the robot joint axis.Equation (1) can be rewritten as a matrix form:(2)
S
X
=
Lwhere S
=
x
1
y
1
z
1
⋮
⋮
⋮
x
n
y
n
z
n, L
=
-
1
⋮
-
1,To solve (2) by using the least square method, we have(3)
X
=
S
T
S
-
1
S
T
LThe trajectory of axis is a space circle whose center and radius are solved by using the space projection method [28] and its objective function is expressed as(4)
F
x
0
,
y
0
,
z
0
,
R
,
λ
=
∑
j
=
1
n
x
j
′
-
x
0
2
+
y
j
′
-
y
0
2
+
z
j
′
-
z
0
2
-
R
2
+
λ
A
x
0
+
B
y
0
+
C
z
0
+
1where P
j
′
(
x
j
′
,
y
j
′
,
z
j
′
) is the projection point of P
j
(
x
j
,
y
j
,
z
j
) on the fitting plane, (x
0,y
0,z
0) andRare the coordinates of sphere center and sphere radius, respectively.λ is a spatial constraint coefficient. (x
0,y
0,z
0) andR can be obtained by using the extremum method. Therefore, the spatial position of joint axis can be determined by the plane normal vector X and sphere center coordinates (x
0,y
0,z
0).


#### 2.1.3. Establishment of Coordinate System

The coordinate system{
T
i
} (i = 1, 2, …,H,His the number of joints) of nonparallel joints is established based on the axis equations of jointi (J
i) and jointi+1 (J
i+1) and shown in Figure 1. The axis direction vector n
i
+
1 ofJ
i+1 is set to the axisz
i. The common perpendicular of the axes ofJ
i andJ
i+1 is set to the axisx
i. The intersection of the axesx
i andz
i is the originO
i. The axisy
i is obtained by the right hand rule: y
i
=
z
i
×
x
i.C
i is the center of axis trajectory whenJ
i rotates separately.Figure 1
Coordinate system establishment for nonparallel joints.For parallel joints the joint coordinate system{
T
i
-
1
} and the axis direction vector n
i
+
1 ofJ
i+1 are shown in Figure 2. The intersection of the auxiliary plane A
x
+
B
y
+
C
z
+
1
=
0 and the axis ofJ
i+1 is the originO
i of the coordinate system {
T
i
}. The auxiliary vector h is the unit one from O
i
-
1 toO
i. The vector n
i
+
1 is set to the axisz
i. The common perpendicular of vector h and the axisz
i is set toy
i
.The axisx
i is obtained by the right-hand rule: x
i
=
z
i
×
y
i.Figure 2
Coordinate system establishment for parallel joints.


#### 2.1.4. Calibration Procedure

The MDH method is used to establish the geometric model of six-joint series robot, in which the parameters are link lengtha
i, offset distanced
i, joint angle θ
i, and joint torsion angle α
i and β
i. α
i expresses the rotation angle from the axis z
i to z
i
+
1 around the axis x
i and β
i is the rotation angle from the axis z
i to z
i
+
1 around the axis y
i. The actual values of geometric parametersa
i,d
i, θ
i, α
i, and β
i can be calibrated according to the conversion relation between the coordinate systems of robot adjacent joints and their geometry definition. The calibration flow is as follows.Step 1.
Install the target ball of laser tracker at the robot end-effector and set the robot at the initial zero position.Step 2.
Control the first joint to rotate in a measurable range and measure the target ball position by using laser tracker. Control this joint to return to the initial zero position after the motion of axis ends. Then, control the next joint to rotate separately until the last one. The target ball position of each joint is measured and the data are recorded.Step 3.
Fit the trajectory circleC
i of each jointiby Spatial Analyzer (SA) software provided by Leica laser tracker and calculate the normal vector X
i of the plane in which the space circle lies.Step 4.
Establish the coordinate system of each joint according to Section2.1.3. Calculate and obtain the calibrated geometric parameters of each joint according to the geometry definition.


### 2.2. Calculation of Robot End-Effector Position

Based on the MDH method, the translation matrixT
i
i
-
1 from the i
-
1 coordinate toicoordinate can be calculated as [26](5)
i
i
-
1
T
=
R
o
t
z
,
θ
i
∗
T
r
a
n
s
z
,
d
i
∗
T
r
a
n
s
x
,
a
i
∗
R
o
t
x
,
α
i
∗
R
o
t
y
,
β
i
=
c
θ
i
c
β
i
-
s
θ
i
s
α
i
s
β
i
-
s
θ
i
c
α
i
c
θ
i
s
β
i
+
s
θ
i
s
α
i
c
β
i
a
i
c
θ
i
s
θ
i
c
β
i
+
c
θ
i
s
α
i
s
β
i
c
θ
i
c
α
i
s
θ
i
s
β
i
-
c
θ
i
s
α
i
c
β
i
a
i
s
θ
i
-
c
α
i
s
β
i
s
α
i
c
α
i
c
β
i
d
i
0
0
0
1wherecands stand for cos and sin, respectively. When the adjacent joints are parallel or near parallel, d
i
=
0, or else β
i
=
0.The transformation matrix from the base coordinate to the end effector can be expressed as(6)
T
N
0
=
∏
i
=
1
H
T
i
i
-
1
=
n
x
o
x
a
x
P
x
n
y
o
x
a
x
P
y
n
z
o
z
a
z
P
z
0
0
0
1From (5) and (6) it is shown that the end-effector position P
=
P
x
P
y
P
z
T can be calculated by the geometric parametersa
i,d
i, θ
i, α
i, and β
i
(
i
=
1,2
,
…
,
H
) and be expressed as(7)
P
=
f
a
i
,
α
i
,
d
i
,
β
i
,
θ
iIt can be seen that the functional relationship between the robot geometric parameters and its end-effector position is highly nonlinear.


## 2.1. Geometric Parameters Calibration

### 2.1.1. Robot Geometric Model

A number of different approaches including Denavit-Hartenberg (DH) method, S model method, Complete and Parametrically Continuous (CPC) model method, zero reference model method, and Product of Exponentials (POE) formula method exist for developing the geometric model of a robot manipulator [25]. The most popular model has been the DH model which is based on homogeneous transformation matrices [3]. This procedure consists of establishing coordinate systems on each joint axis. Each coordinate system is then related to the next through a specific set of coefficients in the homogeneous transformation matrices. To be suitable for robot calibration, the robot geometric model must satisfy three basic requirements: model completeness, parameter minimality, and model continuity. The conventional DH is not enough to achieve a complete model for geometric parameters calibration. In addition, when the adjacent joint axes are parallel, little changes at this state may cause dramatic changes at the common perpendicular between the two axes. The changes may result in robot singularity, which lead to discontinuity of the geometric model. Therefore, a modified DH (MDH) method was proposed, the core of which was to add an angle parameter to the DH model to describe the position relationship between adjacent parallel axes, to cancel the incompleteness problem [25–27]. In this paper, MDH model is adopted to establish the geometric parameters model of six-joint series robot for avoiding the singularity of conventional DH model.


### 2.1.2. Acquisition of Joint Axis Space Position

For obtaining the spatial position of the robot joint axis, each axis is controlled to move separately and the target ball position is measured by the laser tracker. The spatial plane in which the trajectory is located is fitted.Supposing that the coordinates of measured points areP
j
(
x
j
,
y
j
,
z
j
)
,
j
=
1,2
,
⋯
n,n is the number of measured points. The space plane equation is expressed as(1)
A
x
+
B
y
+
C
z
+
1
=
0
X
=
A
B
C
T is the normal vector of the plane, which is also the direction vector n of the robot joint axis.Equation (1) can be rewritten as a matrix form:(2)
S
X
=
Lwhere S
=
x
1
y
1
z
1
⋮
⋮
⋮
x
n
y
n
z
n, L
=
-
1
⋮
-
1,To solve (2) by using the least square method, we have(3)
X
=
S
T
S
-
1
S
T
LThe trajectory of axis is a space circle whose center and radius are solved by using the space projection method [28] and its objective function is expressed as(4)
F
x
0
,
y
0
,
z
0
,
R
,
λ
=
∑
j
=
1
n
x
j
′
-
x
0
2
+
y
j
′
-
y
0
2
+
z
j
′
-
z
0
2
-
R
2
+
λ
A
x
0
+
B
y
0
+
C
z
0
+
1where P
j
′
(
x
j
′
,
y
j
′
,
z
j
′
) is the projection point of P
j
(
x
j
,
y
j
,
z
j
) on the fitting plane, (x
0,y
0,z
0) andRare the coordinates of sphere center and sphere radius, respectively.λ is a spatial constraint coefficient. (x
0,y
0,z
0) andR can be obtained by using the extremum method. Therefore, the spatial position of joint axis can be determined by the plane normal vector X and sphere center coordinates (x
0,y
0,z
0).


### 2.1.3. Establishment of Coordinate System

The coordinate system{
T
i
} (i = 1, 2, …,H,His the number of joints) of nonparallel joints is established based on the axis equations of jointi (J
i) and jointi+1 (J
i+1) and shown in Figure 1. The axis direction vector n
i
+
1 ofJ
i+1 is set to the axisz
i. The common perpendicular of the axes ofJ
i andJ
i+1 is set to the axisx
i. The intersection of the axesx
i andz
i is the originO
i. The axisy
i is obtained by the right hand rule: y
i
=
z
i
×
x
i.C
i is the center of axis trajectory whenJ
i rotates separately.Figure 1
Coordinate system establishment for nonparallel joints.For parallel joints the joint coordinate system{
T
i
-
1
} and the axis direction vector n
i
+
1 ofJ
i+1 are shown in Figure 2. The intersection of the auxiliary plane A
x
+
B
y
+
C
z
+
1
=
0 and the axis ofJ
i+1 is the originO
i of the coordinate system {
T
i
}. The auxiliary vector h is the unit one from O
i
-
1 toO
i. The vector n
i
+
1 is set to the axisz
i. The common perpendicular of vector h and the axisz
i is set toy
i
.The axisx
i is obtained by the right-hand rule: x
i
=
z
i
×
y
i.Figure 2
Coordinate system establishment for parallel joints.


### 2.1.4. Calibration Procedure

The MDH method is used to establish the geometric model of six-joint series robot, in which the parameters are link lengtha
i, offset distanced
i, joint angle θ
i, and joint torsion angle α
i and β
i. α
i expresses the rotation angle from the axis z
i to z
i
+
1 around the axis x
i and β
i is the rotation angle from the axis z
i to z
i
+
1 around the axis y
i. The actual values of geometric parametersa
i,d
i, θ
i, α
i, and β
i can be calibrated according to the conversion relation between the coordinate systems of robot adjacent joints and their geometry definition. The calibration flow is as follows.Step 1.
Install the target ball of laser tracker at the robot end-effector and set the robot at the initial zero position.Step 2.
Control the first joint to rotate in a measurable range and measure the target ball position by using laser tracker. Control this joint to return to the initial zero position after the motion of axis ends. Then, control the next joint to rotate separately until the last one. The target ball position of each joint is measured and the data are recorded.Step 3.
Fit the trajectory circleC
i of each jointiby Spatial Analyzer (SA) software provided by Leica laser tracker and calculate the normal vector X
i of the plane in which the space circle lies.Step 4.
Establish the coordinate system of each joint according to Section2.1.3. Calculate and obtain the calibrated geometric parameters of each joint according to the geometry definition.


## 2.1.1. Robot Geometric Model

A number of different approaches including Denavit-Hartenberg (DH) method, S model method, Complete and Parametrically Continuous (CPC) model method, zero reference model method, and Product of Exponentials (POE) formula method exist for developing the geometric model of a robot manipulator [25]. The most popular model has been the DH model which is based on homogeneous transformation matrices [3]. This procedure consists of establishing coordinate systems on each joint axis. Each coordinate system is then related to the next through a specific set of coefficients in the homogeneous transformation matrices. To be suitable for robot calibration, the robot geometric model must satisfy three basic requirements: model completeness, parameter minimality, and model continuity. The conventional DH is not enough to achieve a complete model for geometric parameters calibration. In addition, when the adjacent joint axes are parallel, little changes at this state may cause dramatic changes at the common perpendicular between the two axes. The changes may result in robot singularity, which lead to discontinuity of the geometric model. Therefore, a modified DH (MDH) method was proposed, the core of which was to add an angle parameter to the DH model to describe the position relationship between adjacent parallel axes, to cancel the incompleteness problem [25–27]. In this paper, MDH model is adopted to establish the geometric parameters model of six-joint series robot for avoiding the singularity of conventional DH model.


## 2.1.2. Acquisition of Joint Axis Space Position

For obtaining the spatial position of the robot joint axis, each axis is controlled to move separately and the target ball position is measured by the laser tracker. The spatial plane in which the trajectory is located is fitted.Supposing that the coordinates of measured points areP
j
(
x
j
,
y
j
,
z
j
)
,
j
=
1,2
,
⋯
n,n is the number of measured points. The space plane equation is expressed as(1)
A
x
+
B
y
+
C
z
+
1
=
0
X
=
A
B
C
T is the normal vector of the plane, which is also the direction vector n of the robot joint axis.Equation (1) can be rewritten as a matrix form:(2)
S
X
=
Lwhere S
=
x
1
y
1
z
1
⋮
⋮
⋮
x
n
y
n
z
n, L
=
-
1
⋮
-
1,To solve (2) by using the least square method, we have(3)
X
=
S
T
S
-
1
S
T
LThe trajectory of axis is a space circle whose center and radius are solved by using the space projection method [28] and its objective function is expressed as(4)
F
x
0
,
y
0
,
z
0
,
R
,
λ
=
∑
j
=
1
n
x
j
′
-
x
0
2
+
y
j
′
-
y
0
2
+
z
j
′
-
z
0
2
-
R
2
+
λ
A
x
0
+
B
y
0
+
C
z
0
+
1where P
j
′
(
x
j
′
,
y
j
′
,
z
j
′
) is the projection point of P
j
(
x
j
,
y
j
,
z
j
) on the fitting plane, (x
0,y
0,z
0) andRare the coordinates of sphere center and sphere radius, respectively.λ is a spatial constraint coefficient. (x
0,y
0,z
0) andR can be obtained by using the extremum method. Therefore, the spatial position of joint axis can be determined by the plane normal vector X and sphere center coordinates (x
0,y
0,z
0).


## 2.1.3. Establishment of Coordinate System

The coordinate system{
T
i
} (i = 1, 2, …,H,His the number of joints) of nonparallel joints is established based on the axis equations of jointi (J
i) and jointi+1 (J
i+1) and shown in Figure 1. The axis direction vector n
i
+
1 ofJ
i+1 is set to the axisz
i. The common perpendicular of the axes ofJ
i andJ
i+1 is set to the axisx
i. The intersection of the axesx
i andz
i is the originO
i. The axisy
i is obtained by the right hand rule: y
i
=
z
i
×
x
i.C
i is the center of axis trajectory whenJ
i rotates separately.Figure 1
Coordinate system establishment for nonparallel joints.For parallel joints the joint coordinate system{
T
i
-
1
} and the axis direction vector n
i
+
1 ofJ
i+1 are shown in Figure 2. The intersection of the auxiliary plane A
x
+
B
y
+
C
z
+
1
=
0 and the axis ofJ
i+1 is the originO
i of the coordinate system {
T
i
}. The auxiliary vector h is the unit one from O
i
-
1 toO
i. The vector n
i
+
1 is set to the axisz
i. The common perpendicular of vector h and the axisz
i is set toy
i
.The axisx
i is obtained by the right-hand rule: x
i
=
z
i
×
y
i.Figure 2
Coordinate system establishment for parallel joints.


## 2.1.4. Calibration Procedure

The MDH method is used to establish the geometric model of six-joint series robot, in which the parameters are link lengtha
i, offset distanced
i, joint angle θ
i, and joint torsion angle α
i and β
i. α
i expresses the rotation angle from the axis z
i to z
i
+
1 around the axis x
i and β
i is the rotation angle from the axis z
i to z
i
+
1 around the axis y
i. The actual values of geometric parametersa
i,d
i, θ
i, α
i, and β
i can be calibrated according to the conversion relation between the coordinate systems of robot adjacent joints and their geometry definition. The calibration flow is as follows.Step 1.
Install the target ball of laser tracker at the robot end-effector and set the robot at the initial zero position.Step 2.
Control the first joint to rotate in a measurable range and measure the target ball position by using laser tracker. Control this joint to return to the initial zero position after the motion of axis ends. Then, control the next joint to rotate separately until the last one. The target ball position of each joint is measured and the data are recorded.Step 3.
Fit the trajectory circleC
i of each jointiby Spatial Analyzer (SA) software provided by Leica laser tracker and calculate the normal vector X
i of the plane in which the space circle lies.Step 4.
Establish the coordinate system of each joint according to Section2.1.3. Calculate and obtain the calibrated geometric parameters of each joint according to the geometry definition.


## 2.2. Calculation of Robot End-Effector Position

Based on the MDH method, the translation matrixT
i
i
-
1 from the i
-
1 coordinate toicoordinate can be calculated as [26](5)
i
i
-
1
T
=
R
o
t
z
,
θ
i
∗
T
r
a
n
s
z
,
d
i
∗
T
r
a
n
s
x
,
a
i
∗
R
o
t
x
,
α
i
∗
R
o
t
y
,
β
i
=
c
θ
i
c
β
i
-
s
θ
i
s
α
i
s
β
i
-
s
θ
i
c
α
i
c
θ
i
s
β
i
+
s
θ
i
s
α
i
c
β
i
a
i
c
θ
i
s
θ
i
c
β
i
+
c
θ
i
s
α
i
s
β
i
c
θ
i
c
α
i
s
θ
i
s
β
i
-
c
θ
i
s
α
i
c
β
i
a
i
s
θ
i
-
c
α
i
s
β
i
s
α
i
c
α
i
c
β
i
d
i
0
0
0
1wherecands stand for cos and sin, respectively. When the adjacent joints are parallel or near parallel, d
i
=
0, or else β
i
=
0.The transformation matrix from the base coordinate to the end effector can be expressed as(6)
T
N
0
=
∏
i
=
1
H
T
i
i
-
1
=
n
x
o
x
a
x
P
x
n
y
o
x
a
x
P
y
n
z
o
z
a
z
P
z
0
0
0
1From (5) and (6) it is shown that the end-effector position P
=
P
x
P
y
P
z
T can be calculated by the geometric parametersa
i,d
i, θ
i, α
i, and β
i
(
i
=
1,2
,
…
,
H
) and be expressed as(7)
P
=
f
a
i
,
α
i
,
d
i
,
β
i
,
θ
iIt can be seen that the functional relationship between the robot geometric parameters and its end-effector position is highly nonlinear.


## 3. Uncertainty Estimation

Measurement uncertainty is a nonnegative parameter characterizing the dispersion of the quantity values being attributed to a measurand based on the information used. GUM is the internationally accepted master document for evaluating uncertainty and it requires the use of a first-order Taylor series expansion for propagating uncertainties. When the mathematical model of measurand is highly nonlinear the use of this linear approximation may be inadequate, MCS is recommended. Furthermore, in order to ensure the stability of required results, the adaptive MCS allows adjusting the number of MCS trials [20].


### 3.1. Uncertainty Estimation of Geometric Parameters

According to the GUM method, the uncertainty of the measurement results generally includes several components, which can be classified into Type-A and Type-B uncertainties based on the numerical evaluation methods. Type-A uncertainty is evaluated by statistical approach and characterized by the experimental standard deviation. Type-B uncertainty is obtained by nonstatistical approach and characterized by empirical or standard deviation estimated from the assumed probability distribution. For the situation with sufficient and independent test data, Type-A uncertainty is more objective than other techniques for uncertainty quantification [19]. A series of measurements are conducted and the standard deviation of tests is obtained; then Type-A uncertainty can be represented as(8)
u
A
=
u
b
-
=
∑
i
=
1
N
b
i
-
b
-
2
N
N
-
1whereNis the test times,b
i represents the result of theith test, and b
- represents the mean-value of b
i
(
i
=
1,2
,
…
,
N
). In our experiments, Type-A uncertainty is used to calculate the calibration uncertainty of robot geometric parameters [23].


### 3.2. Uncertainty Estimation of Robot End-Effector Position

#### 3.2.1. MCS

Typically, a measurandY is not measured directly and it is determined bym other quantities X
1
,
X
2
,
…
,
X
m through a functional relationshipf:(9)
Y
=
f
X
1
,
X
2
,
…
,
X
mIf the functional relationship betweenY and its input quantities X
i
(
i
=
1
,
…
,
m
) is highly nonlinear and a first-order Taylor expansion of the relationship is not an acceptable approximation, in such cases, MCS is required to evaluate the uncertainty [20]. The heart of MCS is repeatedly sampling from probability density function (PDF) for X
i
(
i
=
1
,
…
,
m
) and evaluating model (9) in each case.MCS can be summarized as follows.(1) Set the numberMof MCS to be conducted.(2) GenerateMvectors by sampling from the assigned PDFs, as realizations of the set ofm input quantitiesX
i.(3) Form the corresponding model value ofY and yieldM model values for each such vector.(4) Sort theseM model values into strictly increasing order and use the sorted model values to provide the distribution functionG.(5) UseG to form an estimatey ofY and the standard uncertaintyu(y) associated withy.(6) UseG to form an appropriate coverage interval forYfor a stipulated coverage probabilityp.


#### 3.2.2. Uncertainty Estimation of End-Effector Position Based on AMCS

It can be seen that the numberM of MCS needs to be selected in advance and there is no direct control over the quality of the numerical results provided by MCS. It cannot guarantee that any specific preassigned numberM will suffice. As the Monte Carlo trials increases, the average quality of the system is improved by decreasing the statistical uncertainty of the quantities of interest. The statistical reliability of the measurement uncertainty evaluated by MCS is of utmost importance and AMCS is an effective method to ensure the calculation reliability.From (5), (6), and (7) we can see that the relationship of the end-effector position P
(
P
x
,
P
y
,
P
z
) and the geometric parametersa
i,d
i, θ
i, α
i, and β
i
(
i
=
1,2
,
…
,
H
) is highly nonlinear. Considering that a fixed numberMcannot guarantee the reliability of the calibration uncertainty for robot end-effector position, AMCS is proposed to estimate the uncertainty of robot end-effector position. A basic implementation of AMCS involves carrying out an increasing number of MCS trials until the important results such as expected valueP, standard uncertaintyu(P), and lower and higher coverage interval endpoints P
l
o
w and P
h
i
g
h have stabilized if twice their standard deviations are equal or smaller than the associated numerical toleranceδ, which relies on the number n
d
i
g of significant digits of the results [17, 20]. The AMCS flowchart for uncertainty estimation of robot end-effector position is shown in Figure 3.Figure 3
AMCS for uncertainty estimation of end-effector positions.


## 3.1. Uncertainty Estimation of Geometric Parameters

According to the GUM method, the uncertainty of the measurement results generally includes several components, which can be classified into Type-A and Type-B uncertainties based on the numerical evaluation methods. Type-A uncertainty is evaluated by statistical approach and characterized by the experimental standard deviation. Type-B uncertainty is obtained by nonstatistical approach and characterized by empirical or standard deviation estimated from the assumed probability distribution. For the situation with sufficient and independent test data, Type-A uncertainty is more objective than other techniques for uncertainty quantification [19]. A series of measurements are conducted and the standard deviation of tests is obtained; then Type-A uncertainty can be represented as(8)
u
A
=
u
b
-
=
∑
i
=
1
N
b
i
-
b
-
2
N
N
-
1whereNis the test times,b
i represents the result of theith test, and b
- represents the mean-value of b
i
(
i
=
1,2
,
…
,
N
). In our experiments, Type-A uncertainty is used to calculate the calibration uncertainty of robot geometric parameters [23].


## 3.2. Uncertainty Estimation of Robot End-Effector Position

### 3.2.1. MCS

Typically, a measurandY is not measured directly and it is determined bym other quantities X
1
,
X
2
,
…
,
X
m through a functional relationshipf:(9)
Y
=
f
X
1
,
X
2
,
…
,
X
mIf the functional relationship betweenY and its input quantities X
i
(
i
=
1
,
…
,
m
) is highly nonlinear and a first-order Taylor expansion of the relationship is not an acceptable approximation, in such cases, MCS is required to evaluate the uncertainty [20]. The heart of MCS is repeatedly sampling from probability density function (PDF) for X
i
(
i
=
1
,
…
,
m
) and evaluating model (9) in each case.MCS can be summarized as follows.(1) Set the numberMof MCS to be conducted.(2) GenerateMvectors by sampling from the assigned PDFs, as realizations of the set ofm input quantitiesX
i.(3) Form the corresponding model value ofY and yieldM model values for each such vector.(4) Sort theseM model values into strictly increasing order and use the sorted model values to provide the distribution functionG.(5) UseG to form an estimatey ofY and the standard uncertaintyu(y) associated withy.(6) UseG to form an appropriate coverage interval forYfor a stipulated coverage probabilityp.


### 3.2.2. Uncertainty Estimation of End-Effector Position Based on AMCS

It can be seen that the numberM of MCS needs to be selected in advance and there is no direct control over the quality of the numerical results provided by MCS. It cannot guarantee that any specific preassigned numberM will suffice. As the Monte Carlo trials increases, the average quality of the system is improved by decreasing the statistical uncertainty of the quantities of interest. The statistical reliability of the measurement uncertainty evaluated by MCS is of utmost importance and AMCS is an effective method to ensure the calculation reliability.From (5), (6), and (7) we can see that the relationship of the end-effector position P
(
P
x
,
P
y
,
P
z
) and the geometric parametersa
i,d
i, θ
i, α
i, and β
i
(
i
=
1,2
,
…
,
H
) is highly nonlinear. Considering that a fixed numberMcannot guarantee the reliability of the calibration uncertainty for robot end-effector position, AMCS is proposed to estimate the uncertainty of robot end-effector position. A basic implementation of AMCS involves carrying out an increasing number of MCS trials until the important results such as expected valueP, standard uncertaintyu(P), and lower and higher coverage interval endpoints P
l
o
w and P
h
i
g
h have stabilized if twice their standard deviations are equal or smaller than the associated numerical toleranceδ, which relies on the number n
d
i
g of significant digits of the results [17, 20]. The AMCS flowchart for uncertainty estimation of robot end-effector position is shown in Figure 3.Figure 3
AMCS for uncertainty estimation of end-effector positions.


## 3.2.1. MCS

Typically, a measurandY is not measured directly and it is determined bym other quantities X
1
,
X
2
,
…
,
X
m through a functional relationshipf:(9)
Y
=
f
X
1
,
X
2
,
…
,
X
mIf the functional relationship betweenY and its input quantities X
i
(
i
=
1
,
…
,
m
) is highly nonlinear and a first-order Taylor expansion of the relationship is not an acceptable approximation, in such cases, MCS is required to evaluate the uncertainty [20]. The heart of MCS is repeatedly sampling from probability density function (PDF) for X
i
(
i
=
1
,
…
,
m
) and evaluating model (9) in each case.MCS can be summarized as follows.(1) Set the numberMof MCS to be conducted.(2) GenerateMvectors by sampling from the assigned PDFs, as realizations of the set ofm input quantitiesX
i.(3) Form the corresponding model value ofY and yieldM model values for each such vector.(4) Sort theseM model values into strictly increasing order and use the sorted model values to provide the distribution functionG.(5) UseG to form an estimatey ofY and the standard uncertaintyu(y) associated withy.(6) UseG to form an appropriate coverage interval forYfor a stipulated coverage probabilityp.


## 3.2.2. Uncertainty Estimation of End-Effector Position Based on AMCS

It can be seen that the numberM of MCS needs to be selected in advance and there is no direct control over the quality of the numerical results provided by MCS. It cannot guarantee that any specific preassigned numberM will suffice. As the Monte Carlo trials increases, the average quality of the system is improved by decreasing the statistical uncertainty of the quantities of interest. The statistical reliability of the measurement uncertainty evaluated by MCS is of utmost importance and AMCS is an effective method to ensure the calculation reliability.From (5), (6), and (7) we can see that the relationship of the end-effector position P
(
P
x
,
P
y
,
P
z
) and the geometric parametersa
i,d
i, θ
i, α
i, and β
i
(
i
=
1,2
,
…
,
H
) is highly nonlinear. Considering that a fixed numberMcannot guarantee the reliability of the calibration uncertainty for robot end-effector position, AMCS is proposed to estimate the uncertainty of robot end-effector position. A basic implementation of AMCS involves carrying out an increasing number of MCS trials until the important results such as expected valueP, standard uncertaintyu(P), and lower and higher coverage interval endpoints P
l
o
w and P
h
i
g
h have stabilized if twice their standard deviations are equal or smaller than the associated numerical toleranceδ, which relies on the number n
d
i
g of significant digits of the results [17, 20]. The AMCS flowchart for uncertainty estimation of robot end-effector position is shown in Figure 3.Figure 3
AMCS for uncertainty estimation of end-effector positions.


## 4. Results and Discussion

### 4.1. Simulation

In order to analyze the effect of robot geometric parameters errors on the uncertainty distributions of end-effector position, the geometric parameters nominal values of AIFUTE ER10L-C10 industrial robot are shown in Table1. When six-joint angles are 0°, -90°, 0°, 0°, 0°, 0°, respectively, the robot is at the initial zero position. 1000 random errors obeying the Gaussian distribution with mean 0 and standard deviation 0.1mm are generated and added to the nominal valuesd
i anda
i, respectively. 1000 random errors obeying Gaussian distribution with mean 0 and standard deviation 0.01rad and 0.001rad are generated and added to the nominal values θ
i, α
i, and β
i, respectively. End-effector positions are calculated according to (5) and (6). And the uncertainty distributions of three-dimensional space of robot end-effector are shown in Figures 4(a) and 4(b), respectively. Comparing Figures 4(a) and 4(b), it can be seen that geometric parameters errors directly influence the absolute positioning accuracy of robot end-effector. The smaller angle parameters errors result in larger uncertainty distribution of robot end-effector position.Table 1
The calibrated values, uncertainty, and the endpoints of 95% coverage intervals of geometric parameters.
Parameters
Nominal values
Calibrated values
Uncertainty
Expanded uncertainty
U-endpoint
L-endpoint
θ
1(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
θ
2(°)
-90.00
-90.0620
0.0112
0.0224
-90.0396
-90.0844
θ
3(°)
0.00
-0.1825
0.0188
0.0376
-0.1449
-0.2201
θ
4(°)
0.00
-0.3125
0.0935
0.1870
-0.1255
-0.4995
θ
5(°)
0.00
-0.5217
0.0146
0.0292
-0.4925
-0.5509
θ
6(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
d
1(mm)
504.00
504.0000
0.0000
0.0000
504.0000
504.0000
d
2(mm)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
d
3(mm)
0.00
0.7937
0.7233
1.4466
2.2403
-0.6529
d
4(mm)
1060.00
1059.7833
0.4637
0.9274
1060.7107
1058.8559
d
5(mm)
0.00
0.7176
0.0580
0.1160
0.8336
0.6016
d
6(mm)
125.00
125.0000
0.0000
0.0000
125.0000
125.0000
a
1(mm)
170.00
166.5975
0.1664
0.3328
166.9303
166.2647
a
2(mm)
780.00
780.2600
0.0637
0.1274
780.3874
780.1326
a
3(mm)
140.00
140.3917
0.2980
0.5960
140.9877
139.7957
a
4(mm)
0.00
0.1621
0.1490
0.2980
0.4601
-0.1359
a
5(mm)
0.00
0.1015
0.1778
0.3556
0.4571
-0.2541
a
6(mm)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
α
1(°)
-90.00
-90.0003
0.0108
0.0216
-89.9787
-90.0219
α
2(°)
0.00
-0.0261
0.0081
0.0162
-0.0099
-0.0423
α
3(°)
-90.00
-89.9810
0.0187
0.0374
-89.9436
-90.0184
α
4(°)
90.00
90.0605
0.3552
0.7104
90.7709
89.3501
α
5(°)
-90.00
-90.0150
0.3361
0.6722
-89.3428
-90.6872
α
6(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
β
2(°)
0.00
0.0234
0.0099
0.0198
0.0432
0.0036Figure 4
The uncertainty distributions of robot end-effector by simulation.
(a)
Standard deviations ofθ
i, α
i, and β
i are 0.01rad
(b)
Standard deviations ofθ
i, α
i, and β
i are 0.001rad
(c)
Standard deviations ofθ
i are 0.01rad
(d)
Standard deviations ofα
i are 0.01radFor further analyzing the effect of each parameter on the end-effector position, 1000 random errors obeying above Gaussian distribution are added to each nominal geometric parameter separately and the end-effector positions are calculated. We find that the angle parameters have much more effect on the uncertainty distribution of end-effector positions than the length parameters. Moreover, the joint angleθ
i has much more effect on end-effector positions than joint torsion angle α
i. The separate effects of joint angle θ
i and torsion angle α
i on the uncertainty distribution of the end-effector position are shown in Figures 4(c) and 4(d), respectively.


### 4.2. Experimental Setup

As case study, AIFUTE ER10L-C10 robot was used for calibration experiments. The experimental setup was shown in Figure5, where Leica 930 laser tracker was used and the target ball was installed at the robot end-effector. Assume that the method was able to locate at the exact position of the robot end-effector in the base frame coordinates. Therefore, the orientation of robot end-effector was not considered [22].Figure 5
Experimental setup.Point-by-point sampling was adopted and the stable point measurement mode was set in order to reduce the effect of motor vibration. During the process of the robot calibration, the experiments environment such as temperature, humid, vibration, and dust was well controlled [2]. Control each joint to rotate in the same environmental conditions. The target ball positions corresponding to calibrated points were measured by Leica 930 laser tracker.


### 4.3. The Calibration of Robot Geometric Parameters

Based on the spatial projection method the center and plane normal vector of the space circle are computed and the joint coordinate systems are established. The MDH model is established and the actual parameters are calibrated according to Section2.1. For comparison the calibrated geometric parameters are also listed in Table 1.


### 4.4. The Calibration Uncertainty of Geometric Parameters

When laser tracker is used to measure and calibrate the robot geometric parameters, the calibration uncertainty mainly results from measurement equipment (Leica AT930), sampling strategy, joint movement range, and measurement environment.The standard uncertainties of calibrated geometric parameters are expressed asu
a
i
,
u
d
i
,
u
α
i
,
u
θ
i
,
u
β
i
(
i
=
1,2
,
…
,
6
) and their uncertainty contributors significant include the following [29–31].(1
) Repeatability. 30 times repeated measurements are recommended [19] and have been conducted. Type-A uncertainty based upon repeated measurements from a controlled process is used to estimate the standard uncertainties u
r
θ
i
,
u
r
d
i
,
u
r
a
i
,
u
r
α
i
,
u
r
β
i
(
i
=
1,2
,
…
,
6
) from repeatability for the set and their values are shown in Table 2.Table 2
The standard uncertainties from repeatability.
joints
u
r
θ
i (°)
u
r
d
i (mm)
u
r
a
i (mm)
u
r
α
i (°)
u
r
β
i (°)
1
0.0000
0.0000
0.0301
0.0007
-
2
0.0043
0.0000
0.0066
0.0014
0.0007
3
0.0019
0.0115
0.0145
0.0010
-
4
0.0032
0.0174
0.0218
0.0060
-
5
0.0060
0.0054
0.0423
0.0059
-
6
0.0000
0.0000
0.0000
0.0000
-(
2) Sampling Strategy. Point-by-point sampling strategy is adopted. The robot geometric parameters can be calibrated with a theoretically necessary minimum number of measured points. But how many measured points are needed reliably to calibrate the geometric parameters is an unsolved problem. In our experiments, the number of measured points is set to 100, 80, 60, 50, 40, 25, 20, and 10, respectively. The standard uncertainties u
s
θ
i
,
u
s
d
i
,
u
s
a
i
,
u
s
α
i
,
u
s
β
i
(
i
=
1,2
,
…
,
6
) resulting from sampling strategy are estimated and their values are shown in Table 3.Table 3
The standard uncertainties from sampling strategy.
joints
u
s
θ
i (°)
u
s
d
i (mm)
u
s
a
i (mm)
u
s
α
i (°)
u
s
β
i (°)
1
0.0000
0.0000
0.0218
0.0004
-
2
0.0014
0.0000
0.0128
0.0005
0.0004
3
0.0012
0.0136
0.0053
0.0005
-
4
0.0004
0.0049
0.0016
0.0012
-
5
0.0050
0.0062
0.0012
0.0015
-
6
0.0000
0.0000
0.0000
0.0000
-(3
) Joint Movement Range. Considering the limitations of the robot space motion range and the laser tracker measurable range, six joints are controlled to rotate in the range of ±10°, ±20°, ±30°, ±40°, ±50°, ±60°, respectively. The standard uncertainties u
m
a
i
,
u
m
d
i
,
u
m
α
i
,
u
m
θ
i
,
u
m
β
i
(
i
=
1,2
,
…
,
6
) resulting from joint movement range are estimated and their values are shown in Table 4.Table 4
The standard uncertainties from joint movement range.
joints
u
m
θ
i (°)
u
m
d
i (mm)
u
m
a
i (mm)
u
m
α
i (°)
u
m
β
i (°)
1
0.0000
0.0000
0.1622
0.0108
-
2
0.0103
0.0000
0.0621
0.0080
0.0099
3
0.0187
0.7231
0.2976
0.0187
-
4
0.0934
0.4633
0.1474
0.3551
-
5
0.0123
0.0574
0.1727
0.3360
-
6
0.0000
0.0000
0.0000
0.0000
-(4
) Measurement Environment. During the robot calibration experiments, the environment such as temperature, humid, and dust was well controlled. In addition, Leica AT930 is equipped with temperature, humidity, and pressure sensors and small changes of temperature, humid, and dust can be compensated in real time. Therefore, the uncertainty from environment is neglected.Above contributors are considered as uncorrelated; the standard uncertaintiesu
a
i
,
u
d
i
,
u
α
i
,
u
θ
i
,
u
β
i of calibrated geometric parameters can be estimated and shown in Table 1. For example, the uncertainty u
a
i of link lengtha
i can be calculated by (10)
u
a
i
=
u
r
a
i
2
+
u
s
a
i
2
+
u
m
a
i
2The expanded uncertaintyUcan be obtained by multiplying the standard uncertaintyu by coverage factork:(11)
U
=
k
uFor calibration and test reports k
=
2 is commonly used and it provides 95% confidence. When Type-A components have been calculated with 30 repeated measurements or more this value is appropriate [19]. For comparison, the upper endpoint (U-endpoint) and lower endpoint (L-endpoint) of 95% coverage intervals for calibrated geometric parameters are also listed in Table 1.


### 4.5. The Uncertainty of the Robot End-Effector Position

In order to investigate the calibration uncertainty of geometric parameters effect on the end-effector position, AMCS is carried out to estimate the uncertainty of the robot end-effector position. Geometric parameter errors obeying the Gaussian distribution in which the mean is the calibrated value of each geometric parameter and standard deviation is its uncertainty are generated. The robot end-effector positions can be calculated by(6). In the application of AMCS, the confidence probabilityp and numerical toleranceδare set to 0.95 and 0.0001, respectively. The numberM = max (100/(1 −p), 104) = 104. AMCS procedure is run. The frequency distributions ofX,Y, andZcoordinates of the robot end-effector are illustrated in Figures 6, 7, and 8, respectively, whereX
i,Y
i, and Z
i
(
i
=
1,2
,
…
,
M
) are assembled into a histogram and it provides an approximation to the PDF ofX,Y, andZ. Left and right vertical solid lines express the 95% confidence intervals. From Figures 6, 7, and 8 it can be seen that AMCS involves carrying out an increasing number of MCS trials until the calculation results have stabilized and the frequency distributions ofX,Y, andZcoordinates are different. The uncertainty distributions of three-dimensional space,xy,xz, andyz section errors of robot end-effector are shown in Figures 9(a), 9(b), 9(c), and 9(d), respectively.Figure 6
Frequency distributions of X coordinates of the end-effector position.Figure 7
Frequency distributions of Y coordinates of the end-effector position.Figure 8
Frequency distributions of Z coordinates of the end-effector position.Figure 9
The uncertainty distributions of robot end-effector by experiments.
(a)
Three-dimensional space
(b)
xy section
(c)
xz section
(d)
yz sectionComparing Figures4 and 9, it can be seen that the uncertainty distributions of end-effector position are different due to different geometric parameter errors. And the robot absolute positioning accuracy can be improved by geometric parameters calibration.


## 4.1. Simulation

In order to analyze the effect of robot geometric parameters errors on the uncertainty distributions of end-effector position, the geometric parameters nominal values of AIFUTE ER10L-C10 industrial robot are shown in Table1. When six-joint angles are 0°, -90°, 0°, 0°, 0°, 0°, respectively, the robot is at the initial zero position. 1000 random errors obeying the Gaussian distribution with mean 0 and standard deviation 0.1mm are generated and added to the nominal valuesd
i anda
i, respectively. 1000 random errors obeying Gaussian distribution with mean 0 and standard deviation 0.01rad and 0.001rad are generated and added to the nominal values θ
i, α
i, and β
i, respectively. End-effector positions are calculated according to (5) and (6). And the uncertainty distributions of three-dimensional space of robot end-effector are shown in Figures 4(a) and 4(b), respectively. Comparing Figures 4(a) and 4(b), it can be seen that geometric parameters errors directly influence the absolute positioning accuracy of robot end-effector. The smaller angle parameters errors result in larger uncertainty distribution of robot end-effector position.Table 1
The calibrated values, uncertainty, and the endpoints of 95% coverage intervals of geometric parameters.
Parameters
Nominal values
Calibrated values
Uncertainty
Expanded uncertainty
U-endpoint
L-endpoint
θ
1(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
θ
2(°)
-90.00
-90.0620
0.0112
0.0224
-90.0396
-90.0844
θ
3(°)
0.00
-0.1825
0.0188
0.0376
-0.1449
-0.2201
θ
4(°)
0.00
-0.3125
0.0935
0.1870
-0.1255
-0.4995
θ
5(°)
0.00
-0.5217
0.0146
0.0292
-0.4925
-0.5509
θ
6(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
d
1(mm)
504.00
504.0000
0.0000
0.0000
504.0000
504.0000
d
2(mm)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
d
3(mm)
0.00
0.7937
0.7233
1.4466
2.2403
-0.6529
d
4(mm)
1060.00
1059.7833
0.4637
0.9274
1060.7107
1058.8559
d
5(mm)
0.00
0.7176
0.0580
0.1160
0.8336
0.6016
d
6(mm)
125.00
125.0000
0.0000
0.0000
125.0000
125.0000
a
1(mm)
170.00
166.5975
0.1664
0.3328
166.9303
166.2647
a
2(mm)
780.00
780.2600
0.0637
0.1274
780.3874
780.1326
a
3(mm)
140.00
140.3917
0.2980
0.5960
140.9877
139.7957
a
4(mm)
0.00
0.1621
0.1490
0.2980
0.4601
-0.1359
a
5(mm)
0.00
0.1015
0.1778
0.3556
0.4571
-0.2541
a
6(mm)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
α
1(°)
-90.00
-90.0003
0.0108
0.0216
-89.9787
-90.0219
α
2(°)
0.00
-0.0261
0.0081
0.0162
-0.0099
-0.0423
α
3(°)
-90.00
-89.9810
0.0187
0.0374
-89.9436
-90.0184
α
4(°)
90.00
90.0605
0.3552
0.7104
90.7709
89.3501
α
5(°)
-90.00
-90.0150
0.3361
0.6722
-89.3428
-90.6872
α
6(°)
0.00
0.0000
0.0000
0.0000
0.0000
0.0000
β
2(°)
0.00
0.0234
0.0099
0.0198
0.0432
0.0036Figure 4
The uncertainty distributions of robot end-effector by simulation.
(a)
Standard deviations ofθ
i, α
i, and β
i are 0.01rad
(b)
Standard deviations ofθ
i, α
i, and β
i are 0.001rad
(c)
Standard deviations ofθ
i are 0.01rad
(d)
Standard deviations ofα
i are 0.01radFor further analyzing the effect of each parameter on the end-effector position, 1000 random errors obeying above Gaussian distribution are added to each nominal geometric parameter separately and the end-effector positions are calculated. We find that the angle parameters have much more effect on the uncertainty distribution of end-effector positions than the length parameters. Moreover, the joint angleθ
i has much more effect on end-effector positions than joint torsion angle α
i. The separate effects of joint angle θ
i and torsion angle α
i on the uncertainty distribution of the end-effector position are shown in Figures 4(c) and 4(d), respectively.


## 4.2. Experimental Setup

As case study, AIFUTE ER10L-C10 robot was used for calibration experiments. The experimental setup was shown in Figure5, where Leica 930 laser tracker was used and the target ball was installed at the robot end-effector. Assume that the method was able to locate at the exact position of the robot end-effector in the base frame coordinates. Therefore, the orientation of robot end-effector was not considered [22].Figure 5
Experimental setup.Point-by-point sampling was adopted and the stable point measurement mode was set in order to reduce the effect of motor vibration. During the process of the robot calibration, the experiments environment such as temperature, humid, vibration, and dust was well controlled [2]. Control each joint to rotate in the same environmental conditions. The target ball positions corresponding to calibrated points were measured by Leica 930 laser tracker.


## 4.3. The Calibration of Robot Geometric Parameters

Based on the spatial projection method the center and plane normal vector of the space circle are computed and the joint coordinate systems are established. The MDH model is established and the actual parameters are calibrated according to Section2.1. For comparison the calibrated geometric parameters are also listed in Table 1.


## 4.4. The Calibration Uncertainty of Geometric Parameters

When laser tracker is used to measure and calibrate the robot geometric parameters, the calibration uncertainty mainly results from measurement equipment (Leica AT930), sampling strategy, joint movement range, and measurement environment.The standard uncertainties of calibrated geometric parameters are expressed asu
a
i
,
u
d
i
,
u
α
i
,
u
θ
i
,
u
β
i
(
i
=
1,2
,
…
,
6
) and their uncertainty contributors significant include the following [29–31].(1
) Repeatability. 30 times repeated measurements are recommended [19] and have been conducted. Type-A uncertainty based upon repeated measurements from a controlled process is used to estimate the standard uncertainties u
r
θ
i
,
u
r
d
i
,
u
r
a
i
,
u
r
α
i
,
u
r
β
i
(
i
=
1,2
,
…
,
6
) from repeatability for the set and their values are shown in Table 2.Table 2
The standard uncertainties from repeatability.
joints
u
r
θ
i (°)
u
r
d
i (mm)
u
r
a
i (mm)
u
r
α
i (°)
u
r
β
i (°)
1
0.0000
0.0000
0.0301
0.0007
-
2
0.0043
0.0000
0.0066
0.0014
0.0007
3
0.0019
0.0115
0.0145
0.0010
-
4
0.0032
0.0174
0.0218
0.0060
-
5
0.0060
0.0054
0.0423
0.0059
-
6
0.0000
0.0000
0.0000
0.0000
-(
2) Sampling Strategy. Point-by-point sampling strategy is adopted. The robot geometric parameters can be calibrated with a theoretically necessary minimum number of measured points. But how many measured points are needed reliably to calibrate the geometric parameters is an unsolved problem. In our experiments, the number of measured points is set to 100, 80, 60, 50, 40, 25, 20, and 10, respectively. The standard uncertainties u
s
θ
i
,
u
s
d
i
,
u
s
a
i
,
u
s
α
i
,
u
s
β
i
(
i
=
1,2
,
…
,
6
) resulting from sampling strategy are estimated and their values are shown in Table 3.Table 3
The standard uncertainties from sampling strategy.
joints
u
s
θ
i (°)
u
s
d
i (mm)
u
s
a
i (mm)
u
s
α
i (°)
u
s
β
i (°)
1
0.0000
0.0000
0.0218
0.0004
-
2
0.0014
0.0000
0.0128
0.0005
0.0004
3
0.0012
0.0136
0.0053
0.0005
-
4
0.0004
0.0049
0.0016
0.0012
-
5
0.0050
0.0062
0.0012
0.0015
-
6
0.0000
0.0000
0.0000
0.0000
-(3
) Joint Movement Range. Considering the limitations of the robot space motion range and the laser tracker measurable range, six joints are controlled to rotate in the range of ±10°, ±20°, ±30°, ±40°, ±50°, ±60°, respectively. The standard uncertainties u
m
a
i
,
u
m
d
i
,
u
m
α
i
,
u
m
θ
i
,
u
m
β
i
(
i
=
1,2
,
…
,
6
) resulting from joint movement range are estimated and their values are shown in Table 4.Table 4
The standard uncertainties from joint movement range.
joints
u
m
θ
i (°)
u
m
d
i (mm)
u
m
a
i (mm)
u
m
α
i (°)
u
m
β
i (°)
1
0.0000
0.0000
0.1622
0.0108
-
2
0.0103
0.0000
0.0621
0.0080
0.0099
3
0.0187
0.7231
0.2976
0.0187
-
4
0.0934
0.4633
0.1474
0.3551
-
5
0.0123
0.0574
0.1727
0.3360
-
6
0.0000
0.0000
0.0000
0.0000
-(4
) Measurement Environment. During the robot calibration experiments, the environment such as temperature, humid, and dust was well controlled. In addition, Leica AT930 is equipped with temperature, humidity, and pressure sensors and small changes of temperature, humid, and dust can be compensated in real time. Therefore, the uncertainty from environment is neglected.Above contributors are considered as uncorrelated; the standard uncertaintiesu
a
i
,
u
d
i
,
u
α
i
,
u
θ
i
,
u
β
i of calibrated geometric parameters can be estimated and shown in Table 1. For example, the uncertainty u
a
i of link lengtha
i can be calculated by (10)
u
a
i
=
u
r
a
i
2
+
u
s
a
i
2
+
u
m
a
i
2The expanded uncertaintyUcan be obtained by multiplying the standard uncertaintyu by coverage factork:(11)
U
=
k
uFor calibration and test reports k
=
2 is commonly used and it provides 95% confidence. When Type-A components have been calculated with 30 repeated measurements or more this value is appropriate [19]. For comparison, the upper endpoint (U-endpoint) and lower endpoint (L-endpoint) of 95% coverage intervals for calibrated geometric parameters are also listed in Table 1.


## 4.5. The Uncertainty of the Robot End-Effector Position

In order to investigate the calibration uncertainty of geometric parameters effect on the end-effector position, AMCS is carried out to estimate the uncertainty of the robot end-effector position. Geometric parameter errors obeying the Gaussian distribution in which the mean is the calibrated value of each geometric parameter and standard deviation is its uncertainty are generated. The robot end-effector positions can be calculated by(6). In the application of AMCS, the confidence probabilityp and numerical toleranceδare set to 0.95 and 0.0001, respectively. The numberM = max (100/(1 −p), 104) = 104. AMCS procedure is run. The frequency distributions ofX,Y, andZcoordinates of the robot end-effector are illustrated in Figures 6, 7, and 8, respectively, whereX
i,Y
i, and Z
i
(
i
=
1,2
,
…
,
M
) are assembled into a histogram and it provides an approximation to the PDF ofX,Y, andZ. Left and right vertical solid lines express the 95% confidence intervals. From Figures 6, 7, and 8 it can be seen that AMCS involves carrying out an increasing number of MCS trials until the calculation results have stabilized and the frequency distributions ofX,Y, andZcoordinates are different. The uncertainty distributions of three-dimensional space,xy,xz, andyz section errors of robot end-effector are shown in Figures 9(a), 9(b), 9(c), and 9(d), respectively.Figure 6
Frequency distributions of X coordinates of the end-effector position.Figure 7
Frequency distributions of Y coordinates of the end-effector position.Figure 8
Frequency distributions of Z coordinates of the end-effector position.Figure 9
The uncertainty distributions of robot end-effector by experiments.
(a)
Three-dimensional space
(b)
xy section
(c)
xz section
(d)
yz sectionComparing Figures4 and 9, it can be seen that the uncertainty distributions of end-effector position are different due to different geometric parameter errors. And the robot absolute positioning accuracy can be improved by geometric parameters calibration.


## 5. Conclusions

The MDH method is developed to establish the model of six-joint series robot and the joint movement trajectory method is applied to calibrate the robot geometric parameters. It not only can solve the singularity problem for adjacent joints, but also uses the minimum number of parameters to describe the robot completely. The robot positioning accuracy can be improved by geometric parameters calibration. A series of measurements are conducted and the standard deviation of tests is obtained. The uncertainty components significant for robot geometric parameters calibration by Leica 930 laser tracker are analyzed and the uncertainties are calculated based on GUM. Considering that the functional relationship of the robot end-effector position and geometric parameters is highly nonlinear and the computation time and stability are determined by the numberM of MCS trials, AMCS is developed to estimate the uncertainty distribution of end-effector position. It not only can make up for the limitations of GUM for highly nonlinear model, but also has the features that the number of MCS trials can be selected adaptively and the quality of the numerical results can be controlled directly. Simulation and practical example are illustrated and the experiments results verify the proposed method not only can evaluate the uncertainty of geometric parameters, but also estimate the uncertainty distribution of end-effector position in the whole robot workspace. And the uncertainty distribution of end-effector position is greatly affected by geometric parameter errors. It is an effective method to improve the robot absolute positioning accuracy by geometric parameters calibration, especially to reduce angle parameters errors. The proposed method conforms to new generation GPS requirements in which the calibration uncertainty characterizing the reliability of the results is given together when the calibration result is given. It can be popularized to the uncertainty estimation of other kinds of robots’ geometric parameters calibration.



---

*Source: 7830489-2019-06-09.xml*