1. Introduction
The motion accuracy of a manipulator is the basis for a robot to accurately perform its operations, and motion reliability is the key index to evaluate the motion performance of a manipulator [Reference Sun and Chen1,Reference Nazari and Notash2]. The end posture of a manipulator is affected by many influencing factors, and there is a certain coupling relationship among them. Therefore, establishing an accurate and effective multifactors influence model and determining the method of motion reliability calculation can effectively improve the work performance of a manipulator. In addition, fault diagnosis based on the operating data can not only compensate for the deficiencies of that based on an empirical model, but is also effective for identifying new influencing factors.
The error of end-effector pose accuracy of a manipulator mainly includes structural and joint motion parameter errors, such as temperature error, arm shaft deformation, joint flexibility, gear gap, and wear, which indirectly affect the kinematic parameters of a robot [Reference Sun, Liu and Liu3–Reference Jawale and Thorat5]. Jia et al. [Reference Jia, Li, Chen, Sun and Zhang6] proposed a multilayer mapping model to analyze whether motion reliability can quantitatively describe the integrated motion reliability of a space manipulator, and they examined the sensitivity of the reliability factor by conducting simulation experiments. Kim et al. [Reference Kim, Song and Kang7] developed a robot motion uncertainty model including component dimension tolerance and kinematic pair clearance and applied the improved first-order second-moment method and the Monte Carlo method to solve the kinematic accuracy reliability. Pandey and Zhang [Reference Pandey and Zhang8] obtained the extremum distribution and probability density function of a robot position error based on the maximum entropy principle and subsequently conducted a systematic reliability analysis of the robot movement. Du [Reference Du9] developed an envelope approach to analyze the time-dependent mechanism reliability and verified the significant accuracy improvement by a numerical example. Zhang and Du [Reference Zhang and Du10] applied random clearances to mechanism joints in the reliability method and subsequently used an envelope method to calculate the time-dependent reliability. The proposed method was verified by the reliability analysis of three 4-bar function generation mechanisms. Zhao et al. [Reference Zhao, Yan and Wu11] developed a new reliability analysis method to analyze the dynamic reliability and parameter sensitivity of a space manipulator with a harmonic drive. Yu and Wang [Reference Yu and Wang12] proposed a two-stage failure process decomposition method, using which they transformed the time-varying reliability problem into a time-invariant problem having a dynamic structure with uncertainty factors. The effectiveness of the proposed method was verified by examples. To provide an effective tool for the reliability design of complex products, Yu et al. [Reference Yu and Wang13] proposed a general decoupling method for time- and space-variant system reliability (TS-VSR) constraint problems. Moreover, they established a mapping between the design variables and the TS-VSR constraints in combination with exponential function approximation and moving least square method in each cycle. Zhang et al. [Reference Zhang, Xiao, Gao and Chu14] developed the Projection-outline-based active learning Kriging and adaptive importance sampling algorithm based on a combination of the active learning Kriging method and artificial intelligence to improve the computational efficiency of artificial intelligence in the case of a time-consuming performance function. Wu et al. [Reference Wu, Zhang, Liu and Han15] adopted a new calculation method to comprehensively evaluate the reliability of the position accuracies of single-coordinate, single-point, multipoint, and trajectory accuracies of industrial robots. Van et al. [Reference Van, Kang, Suh and Shin16] studied a robust fault diagnosis algorithm for uncertain robot systems based on a neural sliding mode observer, which can reduce the chattering of the sliding mode and ensure the convergence of neural networks in a very short time. In Ref. [Reference Zhang, Xiao, Gao and Fu17], the problem of mixed reliability analysis of random variables and interval variables was studied, and an hybrid reliability analysis under random and interval variables method combining a POAL metamodel and a Kriging metamodel (POAL–Kriging) that does not require establishing a high-quality metamodel was proposed. Zhang et al. [Reference Zhang, Wang and Du18] proposed a probabilistic mechanism synthesis method with interval (or time-varying) motion reliability, which can ensure not only the motion reliability requirements in the entire time interval but also the motion reliability requirements of each point in the interval. The authors verified the method by the synthesis of a four-bar function generator mechanism. Korayem et al. [Reference Korayem, Azimirad, Binabaji and Korayem19] described the configuration of a mobile manipulator with flexible joints (wheels and manipulator) and its full dynamic model, which was based on the coordinated motion of a vehicle and a manipulator. Subsequently, they used the Sobol sensitivity analysis method to determine the optimal value of the flexible joint constant, which was affected by the impact minimization. Siyu et al. [Reference Siyu, Jinyuan, Caiwang and Qibo20] derived the relationship between the center distance error and the tooth gap and subsequently established a multi-degree-of-freedom (DOF) nonlinear dynamic tooth drive system with a moving tooth gap, friction force, and time-varying stiffness. Finally, they used the Runge–Kutta method to solve the nonlinear dynamic system. Erkaya [Reference Erkaya21] and Pandey and Zhang [Reference Pandey and Zhang22] studied the kinematics and dynamics of a robot with six DOFs in a joint space compared to that without a joint space. Ping et al. [Reference Ping, Han, Jiang and Xiao23] proposed a time-varying extremum event evolution method to obtain the time-varying reliability of an arbitrary time interval and an arbitrary failure threshold. Jiang et al. [Reference Jiang, Qiu, Gao, Wang, Yang and Chen24] developed a real-time error estimation guided sampling method, which overcomes the problem that the substitution training strategy of the existing proxy model algorithm does not directly consider the error estimate of the failure probability. This issue leads to the waste or problem of some samples with a large number of computations.
In this study, a UR10 manipulator was taken as the research object, and the influences of the link flexibility factor, joint flexibility factor, and joint clearance factor on the kinematic transformation matrix were comprehensively considered. These factors were combined with the Denavit–Hartenberg (DH) parameters of a robot to establish the mapping model of the influencing factors of motion reliability for directly reflecting the coupling relationship among the multiple factors. In addition, the nonlinear response surface method was used to determine the sensitivity of the influencing factors of motion reliability [Reference Du9], and UR10 was utilized as an example for the simulation test. Finally, the motion accuracy and sensitivity of the DH parameters were verified by the data-driven kernel principal component analysis (KPCA) fault diagnosis method, and the joint rotation fault was used as an example to verify the accuracy of this diagnosis method.
2. Mapping model of influencing factors
2.1. Definition of motion accuracy reliability
${{\textit{P}}^{\rm{A}}}\left( \Theta \right)$ and 
${{\textit{P}}^{\rm{I}}}\left( \Theta \right)$ are defined as the actual and ideal positions at the end-effector of a manipulator, respectively. Positions at the end-effector of a manipulator include space position and space attitude. Assuming that there are m trajectory via points on the entire track, the mathematical description of the motion reliability of the continuous trajectory tracking of an industrial manipulator can be expressed as
\begin{equation}{R_p} = \Pr \left\{ {\bigcap_{i = 1}^{m \to \infty } {\left( {\left| {{P}_i^{\rm{A}}\left( \Theta \right) - {P}_i^{\rm{I}}\left( \Theta \right)} \right| < \lambda } \right)} } \right\}\end{equation}
where 
$\lambda $ is the position accuracy threshold of each trajectory point. x 1, x 2,…, x n are the influencing factors of motion reliability. Θ = (x 1, x 2,…, x n) is the set of influencing factors.
p(x 1, x 2,…x m) is the probability density function related to the position accuracy of the trajectory point [Reference Wu, Zhang, Liu and Han25]. Thus, Eq. (1) can be expressed as
\begin{equation}\mathrm{Pr}\Bigg\{\bigcap_{i=1}^{m\rightarrow\infty}(\big|P_{i}^{A}(\Theta)-P_{i}^{I}(\Theta)\big|<\lambda) \Bigg\}=\int_{0}^{\lambda}\left(\cdots\left(\int_{0}^{\lambda}p(x_1,x_2\cdots x_m)dx_1\right)\cdots\right)dx_m \end{equation}
2.2. The influencing factors of motion accuracy reliability
The end position accuracy of the manipulator is affected by many factors. To analyze its motion reliability, each influencing factor is required to be categorized. The classified influencing factors can effectively characterize the influences of multiple factors on the motion reliability. The connection factor between two joints can be considered as joint clearance, including installation clearance and transmission clearance. The effects of the joint clearance on the motion accuracy can represent the comprehensive influence of various transmission accuracy factors. Based on the above example, a set of factors affecting the motion reliability, represented by link flexibility factors, joint flexibility factors, joint clearance factors, and kinematics parameters of the manipulator, can be formed. Each factor is analyzed separately below.
2.2.1. Geometric parameters of the manipulator
For the 6-DOFs manipulator, the transformation matrix between the coordinate system of the 
${\left( {i - 1} \right)^{th}}$ and 
${i^{th}}$ connecting links can be constructed based on DH parameters 
${\alpha _i}$, 
${a_i}$, 
${d_i}$, and 
${\theta _i}$. The transformation matrix between the i th and i-1th links of robot 
${}_{\ \ \ \ i}^{\ i - 1}{\textit{T}}$ is expressed in Eq. (3).
\begin{equation}{{\textit{T}}_i}^N = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{\cos {\theta _i}} & { - \sin {\theta _i}\cos {\alpha _i}} & {\sin {\theta _i}\sin {\alpha _i}} & {{a_i}\cos {\theta _i}}\\[3pt]{\sin {\theta _i}} & {\cos {\theta _i}\cos {\alpha _i}} & { - \cos {\theta _i}\sin {\alpha _i}} & {{a_i}\sin {\theta _i}}\\[3pt]0 & {\sin {\alpha _i}} & {\cos {\alpha _i}} & {{d_i}}\\[3pt]0 & 0 & 0 & 1\end{array}} \right]\end{equation}
The actual and nominal transformation matrices between the coordinate system of the 
${\left( {i - 1} \right)^{th}}$ and 
${i^{th}}$ connecting links can be represented by 
$T_i^A$ and 
$T_i^N$, respectively, and the deviation can be expressed as follows:
\begin{equation}{\textit{T}}_i^A - {\textit{T}}_i^N = \frac{{\partial {{\textit{T}}_i}}}{{\partial {\alpha _i}}}\Delta {\alpha _i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {a_i}}}\Delta {a_i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {d_i}}}\Delta {d_i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {\theta _i}}}\Delta {\theta _i} = {\textit{T}}_i^N \cdot \delta {{\textit{T}}_i}\end{equation}
\begin{equation}\qquad\ \delta {{\textit{T}}_i} = {\left( {{\textit{T}}_i^N} \right)^{ - 1}} \cdot \left( {\frac{{\partial {{\textit{T}}_i}}}{{\partial {\alpha _i}}}\Delta {\alpha _i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {a_i}}}\Delta {a_i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {d_i}}}\Delta {d_i} + \frac{{\partial {{\textit{T}}_i}}}{{\partial {\theta _i}}}\Delta {\theta _i}} \right)\end{equation}
$\delta {\textit{T}_i}$ is calculated by Eq. (5). By substituting partial derivatives of 
${{\textit{T}}_i}$ with respect to each parameter and 
${\textit{T}}_i^N$ into the above equation, the specific analytic form of 
$\delta {T_i}$ can be obtained as follows [Reference Jia, Li, Chen, Sun and Zhang6]:
\begin{equation}\delta {T_i} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0 & { - \cos {\alpha _i}\Delta {\theta _i}} & {\sin {\alpha _i}\Delta {\theta _i}} & {\Delta {a_i}}\\[3pt]{\cos {\alpha _i}\Delta {\theta _i}} & 0 & { - \Delta {\alpha _i}} & {{a_i}\cos {\alpha _i}\Delta {\theta _i} + \sin {\alpha _i}\Delta {d_i}}\\[3pt]{ - \sin {\alpha _i}\Delta {\theta _i}} & {\Delta {\alpha _i}} & 0 & { - {a_i}\sin {\alpha _i}\Delta {\theta _i} + \cos {\alpha _i}\Delta {d_i}}\\[3pt]0 & 0 & 0 & 0\end{array}} \right]\end{equation}
As 
$\delta {T_i} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0 & { - {\delta _{iz}}} & { - {\delta _{iy}}} & {{\rm{P}}_{ix}^{\rm{d}}}\\[3pt]{{\delta _{iz}}} & 0 & { - {\delta _{ix}}} & {{\rm{P}}_{iy}^{\rm{d}}}\\[3pt]{ - {\delta _{iy}}} & {{\delta _{ix}}} & 0 & {{\rm{P}}_{iz}^{\rm{d}}}\\[3pt]0 & 0 & 0 & 0\end{array}} \right]$, the corresponding elements are equivalent, and the following Eq. (7) can be obtained:
\begin{equation}\left[ {\begin{array}{c}{{\rm{P}}_i^{\rm{d}}}\\[4pt]{{\delta _i}}\end{array}} \right] = \left[ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{0_{3 \times 1}} & k_i^1 & k_i^2 & k_i^3\\[4pt]k_i^4 & {0_{3 \times 1}} & {0_{3 \times 1}} & k_i^5\end{array}\right]\left[ \begin{array}{c}{\Delta {\alpha _i}}\\[4pt]{\Delta {a_i}}\\[4pt]{\Delta {d_i}}\\[4pt]{\Delta {\theta _i}}\end{array} \right]\end{equation}
where 
${\bf{P}}_i^{\rm{d}}{\rm{ = [P}}_{ix}^{\rm{d}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{P}}_{iy}^{\rm{d}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{P}}_{iz}^{\rm{d}}{{\rm{]}}^{\rm{T}}}$ and 
${\delta _i}{\rm{ = [}}{\delta _{ix}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\delta _{iy}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\delta _{iz}}{{\rm{]}}^{\rm{T}}}$, respectively, represent the position and attitude deviations of the 
$\left\{ i \right\}$ coordinate system relative to the 
$\left\{ {i - 1} \right\}$ coordinate system and k il(l = 1,…,5) 
$\in$ R 3×1.
Where 
$k_i^1 = {\left[ {\begin{array}{c@{\quad}c@{\quad}c}1 & 0 & 0\end{array}} \right]^{\rm T}}$, 
$k_i^2 = {\left[ {\begin{array}{c@{\quad}c@{\quad}c}0 & {\sin {\alpha _i}} & {\cos {\alpha _i}}\end{array}} \right]^{\rm T}}$, 
$k_i^3 = {\left[ {\begin{array}{c@{\quad}c@{\quad}c}0 & {{a_i}\cos {\alpha _i}} & { - {a_i}\sin {\alpha _i}}\end{array}} \right]^{\rm T}}$, 
$k_i^4 = {\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0\end{array}} \right]^{\rm T}}$, and 
$k_i^5 = {\left[ {\begin{array}{c@{\quad}c@{\quad}c} 0 & {\sin {\alpha _i}} & {\cos {\alpha _i}}\end{array}} \right]^{\rm T}}$.
2.2.2. Link flexibility factors
The ideal link of the manipulator is considered as a rigid body, and the coordinate transformation matrix from joint i to joint i + 1 can be represented by 
${\textit{T}}_i^N$. However, the position of joint i + 1 will be changed owing to the flexibility of the link. Oi and Oi* are the ideal and actual positions of joint i + 1, respectively. The flexible factors of the link mainly include its flexible torsion and flexible bending deformations. A schematic representation of the influence of the link flexibility factors is shown in Fig. 1.
Schematic representation of influence of link flexibility factors.
In the figure, 
${{\textit{E}}_i}$ is defined as the transformation matrix from the ideal position to the actual position of the i + 1 joint, and it is expressed in Eq. (8). In Eq. (8), 
${\Phi _i}$ is the antisymmetric matrix of 
${\varphi _i}$. 
${{\boldsymbol{\varepsilon }}_i}$ is the position vector in homogeneous transformation matrix 
${{\textit{E}}_i}$, and 
${\phi _i}{\rm{ = [}}{\varphi _{ix}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _{iy}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _{iz}}{{\rm{]}}^{\rm{T}}}$ and 
${\varepsilon _i}{\rm{ = [}}{\varepsilon _{ix}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _{iy}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _{iz}}{{\rm{]}}^{\rm{T}}}$ are the flexible torsion and flexible bending deformation vectors, respectively.
\begin{equation}{{\textit{E}}_i} = \left[ \begin{array}{c@{\quad}c}{{\Phi _i}} & {{\varepsilon _i}}\\[4pt]0 & 1\end{array} \right]= \left[ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1 & { - {\varphi _{iz}}} & {{\varphi _{iy}}} & {{\varepsilon _i}_x}\\[4pt]{{\varphi _{iz}}} & 1 & { - {\varphi _{ix}}} & {{\varepsilon _i}_y}\\[4pt]{ - {\varphi _{iy}}} & {{\varphi _{ix}}} & 1 & {{\varepsilon _i}_z}\\[4pt]0 & 0 & 0 & 1\end{array} \right]\end{equation}
When there is a link flexibility factor, the homogeneous transformation matrix between adjacent joints is 
${\textit{T}}_{\!\!i}^{A*} = {\textit{T}}_{\!\!i}^N \cdot \delta {\textit{T}}_{\!\!i}^{E*} = {\textit{T}}_{\!\!i}^N \cdot {{\textit{E}}_{\!\!i}} - {\textit{T}}_{\!\!i}^N$, and the specific analytic form of 
$\delta {\textit{T}}_{\!\!i}^{E*}$ can be expressed as
\begin{equation}\delta {\textit{T}}_{\!\!i}^{E*} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - {\varphi _{iz}}} & {{\varphi _{iy}}} & {{\varepsilon _i}_x}\\[3pt]{{\varphi _{iz}}} & 0 & { - {\varphi _{ix}}} & {{\varepsilon _i}_y}\\[3pt]{ - {\varphi _{iy}}} & {{\varphi _{ix}}} & 0 & {{\varepsilon _i}_z}\\[3pt]0 & 0 & 0 & 0\end{array}} \right]\end{equation}
When there is no deviation in geometric parameters, the link flexibility factor can be transformed into the error of kinematic parameters, and there is
\begin{align}\delta {\textit{T}}_{\!\!i}^{E*} & = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - {\varphi _{iz}}} & {{\varphi _{iy}}} & {{\varepsilon _i}_x}\\[4pt]{{\varphi _{iz}}} & 0 & { - {\varphi _{ix}}} & {{\varepsilon _i}_y}\\[4pt]{ - {\varphi _{iy}}} & {{\varphi _{ix}}} & 0 & {{\varepsilon _i}_z}\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\nonumber\\[15pt]& = \delta {T_{\!i}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - \cos {\alpha _i}\Delta {\theta _i}} & {\sin {\alpha _i}\Delta {\theta _i}} & {\Delta {a_i}}\\[4pt]{\cos {\alpha _i}\Delta {\theta _i}} & 0 & { - \Delta {\alpha _i}} & {{a_i}\cos {\alpha _i}\Delta {\theta _i} + \sin {\alpha _i}\Delta {d_i}}\\[4pt]{ - \sin {\alpha _i}\Delta {\theta _i}} & {\Delta {\alpha _i}} & 0 & { - {a_i}\sin {\alpha _i}\Delta {\theta _i} + \cos {\alpha _i}\Delta {d_i}}\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\end{align}
The corresponding elements of the matrix are equal, and Eq. (11) can be obtained.
\begin{equation}\left[ {\begin{array}{c}{\Delta {\alpha _i}}\\[4pt]{\Delta {a_i}}\\[4pt]{\Delta {d_i}}\\[4pt]{\Delta {\theta _i}}\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 1 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & {\sin {\alpha _i}} & {\cos {\alpha _i}}\\[4pt]0 & {\sin {\alpha _i}} & {\cos {\alpha _i}} & 0 & 0 & 0\end{array}} \right]\left[ {\begin{array}{c}{{\varphi _{ix}}}\\[4pt]{{\varphi _{iy}}}\\[4pt]{{\varphi _{iz}}}\\[4pt]{{\varepsilon _i}_x}\\[4pt]{{\varepsilon _i}_y}\\[4pt]{{\varepsilon _i}_z}\end{array}} \right]\end{equation}
2.2.3. Joint flexibility factor
The structure of each joint of the manipulator is shown in Fig. 2. In the case of a small deformation, a reducer gear is equivalent to a torsion spring with a certain stiffness coefficient. 
${J_{mi}}$ is the moment of inertia of a motor rotor at joint 
$i$, 
${N_i}$ is the reduction ratio of a joint reducer, 
${K_i}$ is the stiffness coefficient of a linear torsional spring, 
${\theta _i}$ is the position of a joint output shaft, and 
${q_i}$ is the rotation angle of the mechanical arm connecting link.
Elastic equivalent model of joint.
Owing to the existence of joint elasticity, 
${\theta _i}$ and 
${q_i}$ are inconsistent. 
${\xi _i}$ is defined to characterize the angle influence caused by the joint flexibility. Thus, the transformation matrix 
${{{\textit{B}}}_i}$ caused by the joint flexibility can be expressed as
\begin{equation}{{{\textit{B}}}_i} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1 & { - {\xi _i}} & 0 & 0\\[4pt]{{\xi _i}} & 1 & 0 & 0\\[4pt]0 & 0 & 1 & 0\\[4pt]0 & 0 & 0 & 1\end{array}} \right]\end{equation}
Similarly, the specific analytic form of 
$\delta {\textit{T}}_{\!\!i}^{B*}$ can be expressed as
\begin{equation}\delta {\textit{T}}_{\!\!i}^{B*} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - {\xi _i}} & 0 & 0\\[4pt]{{\xi _i}} & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\end{equation}
When there is no deviation in geometric parameters, the joint flexibility factor can be transformed into the error of kinematic parameters, and there is
\begin{align}\delta {\textit{T}}_{\!\!i}^{B*} & = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - {\xi _i}} & 0 & 0\\[4pt]{{\xi _i}} & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\nonumber\\[10pt]& = \delta {T_{\!i}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - \cos {\alpha _i}\Delta {\theta _i}} & {\sin {\alpha _i}\Delta {\theta _i}} & {\Delta {a_i}}\\[4pt] {\cos {\alpha _i}\Delta {\theta _i}} & 0 & { - \Delta {\alpha _i}} & {{a_i}\cos {\alpha _i}\Delta {\theta _i} + \sin {\alpha _i}\Delta {d_i}}\\[4pt] { - \sin {\alpha _i}\Delta {\theta _i}} & {\Delta {\alpha _i}} & 0 & { - {a_i}\sin {\alpha _i}\Delta {\theta _i} + \cos {\alpha _i}\Delta {d_i}}\\[4pt] 0 & 0 & 0 & 0\end{array}} \right]\end{align}
The corresponding elements of the matrix are equal and can be obtained.
\begin{equation}\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\Delta {\alpha _i}}\\[4pt]{\Delta {a_i}}\\[4pt]{\Delta {d_i}}\\[4pt]{\Delta {\theta _i}}\end{array}} \right] = \left[ {\begin{array}{c}0\\[4pt]0\\[4pt]0\\[4pt]{1/\cos {\alpha _i}}\end{array}} \right] \cdot {\xi _i}\end{equation}
2.2.4. Joint clearance factor
A manipulator typically consists of a series of kinematic pairs, which include moving and rotating pairs. The motion pairs formed by connecting the links are the joints of the manipulator. The existence of a joint clearance affects the motion accuracy of the manipulator. The influences of a joint clearance on the coordinate system are shown in Fig. 3. Due to the existence of a clearance, the fit between the connecting link and the axle hole will inevitably have deviations, which include rotation error 
${{\textit{T}}_{{r_i}}}$ and movement error 
${{\textit{T}}_{{t_i}}}$. Because of these deviations, the coordinate system of the link 
$\left\{ i \right\}$ will inevitably skew and shift from original coordinate system 
$\left\{ i \right\}$ to 
$\left\{ {{i^{\prime}}} \right\}$.
Influences of joint clearance on coordinate system.
Assuming that the components of movement error in the directions of 
${x_i}$, 
${y_i}$, and 
${z_i}$ are 
$\Delta {x_i}$, 
$\Delta {y_i}$, and 
$\Delta {z_i}$, respectively, the rotation error of the rotation about multiple axes is the product of rotation errors 
${T_{{r_x}}}\left( {{x_i},\Delta {\theta _{{x_i}}}} \right)$, 
${T_{{r_y}}}\left( {{y_i},\Delta {\theta _{{y_i}}}} \right)$, and 
${T_{{r_z}}}\left( {{z_i},\Delta {\theta _{{z_i}}}} \right)$ of the rotation about a single axis. The combined consideration of the movement and rotation errors can be expressed as follows:
\begin{equation}{T_r}_{_i} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1 & { - \Delta {\theta _{{z_i}}}} & {\Delta {\theta _{{y_i}}}} & {\Delta {x_i}}\\[4pt]{\Delta {\theta _{{z_i}}}} & 1 & { - \Delta {\theta _{{x_i}}}} & {\Delta {y_i}}\\[4pt]{ - \Delta {\theta _{{y_i}}}} & {\Delta {\theta _{{x_i}}}} & 1 & {\Delta {z_i}}\\[4pt]0 & 0 & 0 & 1\end{array}} \right]\end{equation}
where 
${T_{{r_x}}}\left( {{x_i},\Delta {\theta _{{x_i}}}} \right)$, 
${T_{{r_y}}}\left( {{y_i},\Delta {\theta _{{y_i}}}} \right)$, and 
${T_{{r_z}}}\left( {{z_i},\Delta {\theta _{{z_i}}}} \right)$, respectively, represent the rotation errors of angle 
$\Delta {\theta _{{x_i}}}$ around axis 
${x_i}$, angle 
$\Delta {\theta _{{y_i}}}$ around axis 
${y_i}$, and angle 
$\Delta {\theta _{{z_i}}}$ around axis 
${z_i}$.Where, 
$\Delta {\theta _{{x_i}}}$, 
$\Delta {\theta _{{y_i}}}$, and 
$\Delta {\theta _{{z_i}}}$, respectively, represent the error angles of rotation about the 
${x_i}$, 
${y_i}$, and 
${z_i}$ axis.
Similarly, the specific analytic form of 
$\delta {\textit{T}}_{\!\!i}^{T{\rm{r}}*}$ can be expressed as
\begin{equation}\delta {\textit{T}}_{\!\!i}^{T{\rm{r}}*} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - \Delta {\theta _{{z_i}}}} & {\Delta {\theta _{{y_i}}}} & {\Delta {x_i}}\\[4pt]{\Delta {\theta _{{z_i}}}} & 0 & { - \Delta {\theta _{{x_i}}}} & {\Delta {y_i}}\\[4pt]{ - \Delta {\theta _{{y_i}}}} & {\Delta {\theta _{{x_i}}}} & 0 & {\Delta {z_i}}\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\end{equation}
When there is no deviation in geometric parameters, the joint clearance factor can be transformed into the error of kinematic parameters, and there is
\begin{align}\delta {\textit{T}}_{\!\!i}^{T{\rm{r}}*} & = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - \Delta {\theta _{{z_i}}}} & {\Delta {\theta _{{y_i}}}} & {\Delta {x_i}}\\[4pt]{\Delta {\theta _{{z_i}}}} & 0 & { - \Delta {\theta _{{x_i}}}} & {\Delta {y_i}}\\[4pt]{ - \Delta {\theta _{{y_i}}}} & {\Delta {\theta _{{x_i}}}} & 0 & {\Delta {z_i}}\\[4pt]0 & 0 & 0 & 0\end{array}} \right]\nonumber\\[4pt] & = \delta {T_i} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & { - \cos {\alpha _i}\Delta {\theta _i}} & {\sin {\alpha _i}\Delta {\theta _i}} & {\Delta {a_i}}\\[4pt] {\cos {\alpha _i}\Delta {\theta _i}} & 0 & { - \Delta {\alpha _i}} & {{a_i}\cos {\alpha _i}\Delta {\theta _i} + \sin {\alpha _i}\Delta {d_i}}\\[4pt] { - \sin {\alpha _i}\Delta {\theta _i}} & {\Delta {\alpha _i}} & 0 & { - {a_i}\sin {\alpha _i}\Delta {\theta _i} + \cos {\alpha _i}\Delta {d_i}}\\[4pt] 0 & 0 & 0 & 0\end{array}} \right]\end{align}
The corresponding elements of the matrix are equal and can be obtained.
\begin{equation}\left[ {\begin{array}{c}{\Delta {\alpha _i}}\\[4pt]{\Delta {a_i}}\\[4pt]{\Delta {d_i}}\\[4pt]{\Delta {\theta _i}}\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 1 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & {\sin {\alpha _i}} & {\cos {\alpha _i}}\\[4pt]0 & {\sin {\alpha _i}} & {\cos {\alpha _i}} & 0 & 0 & 0\end{array}} \right]\left[ {\begin{array}{c}{\Delta {\theta _{{x_i}}}}\\[4pt]{\Delta {\theta _{{y_i}}}}\\[4pt]{\Delta {\theta _{{z_i}}}}\\[4pt]{\Delta {x_i}}\\[4pt]{\Delta {y_i}}\\[4pt]{\Delta {z_i}}\end{array}} \right]\end{equation}
2.3. Comprehensive mapping model of reliability influencing factors
Taking the mapping relationship between the kinematic parameter deviation and the end-effector pose accuracy as the intermediate transfer layer, that between each influencing factor and the kinematic accuracy reliability of the manipulator can be determined, as shown in Fig. 4.
Multilayer mapping model.
Based on the analysis of each influencing factor, the influence model of a single factor on the motion accuracy and reliability of the mechanical arm was established. The influence matrix of a multifactor synthesis on the motion accuracy of the manipulator is expressed in Eq. (20).
\begin{equation}T_{\!i}^A = T_{\!i}^N \cdot {E_i} \cdot {B_i} \cdot {T_{{r_i}}}\end{equation}
The deviation can be calculated as expressed in Eq. (21).
\begin{equation}d{T_{\!i}} = T_{\!i}^A - T_{\!i}^N = T_{\!i}^N \cdot \delta {T_{\!i}} = T_{\!i}^N\left( {{E_i} \cdot {B_i} \cdot {T_{{r_i}}} - {I_4}} \right)\end{equation}
3. Sensitivity analysis based on the nonlinear response surface function
Position accuracy of the end-effector is defined as the motion reliability evaluation index. The end position accuracy is obtained from the absolute error, 
$\Delta P$ between the ideal position, 
${P^{\rm{I}}}$ and the actual position, 
${P^{\rm{A}}}$. The limit state function is 
$g(\Theta ) = \left| {{P^{\rm{A}}} - {P^{\rm{I}}}} \right| = \Delta P$. Θ = [x 1 … x n] 
$\in$ R n×1 represents the influencing factors of the motion reliability, and 
$n$ is the number of influencing factors.
By solving the unknown coefficients, the response surface function is developed to estimate the limit state function. Thereby, it can be used to characterize the relationship among the various influencing factors and the accuracy of the end position. For the sake of simplicity, the quadratic polynomials without cross terms were selected as response surface function [Reference Jia, Li, Chen, Sun and Zhang6]. The response surface function [Reference Bucher and Bourgund26] to characterize the relationship among the influencing factors and the pose deviation can be established, as expressed in Eq. (22).
\begin{equation}\hat g(\Theta ) = a + \sum\limits_{i = 1}^n {{b_i}{x_i}} + \sum\limits_{i = 1}^n {{c_i}} x_i^2\end{equation}
where 
$a$,
${b_i}$, and 
${c_i}$ are the coefficients and 
${x_i}$ represents the influencing factor. 
$\hat g(\Theta )$ is the response surface function.
The unknown coefficients 
$a,{b_i},{c_i}$ are defined in vector form as 
${\boldsymbol{\kappa }}$. v is defined as an impact factor matrix. Eq. (22) is expressed in matrix form, as shown in Eq. (23).
\begin{equation}\left\{\begin{array}{l}\boldsymbol{\kappa}=\left[a,b_1,\cdots,b_n,c_1,\cdots,c_n\right]\\[4pt]{\textit{v}}=\left[1 x_1\cdots x_n\ x_1^2\cdots x_n^2\right]^T\\[4pt]\hat{g}(\Theta)=\boldsymbol{\kappa}\cdot{\textit{v}} \end{array}\right.\end{equation}
To solve unknown coefficients vector 
${\boldsymbol{\kappa }}$, test points 
$(\overline x {}_i,g(\overline x {}_i))$ are constructed. 
$\overline x {}_i$ is the variable value of the constructed test point, and 
$g(\overline x {}_i)$ is the limit state function value of the test point. When constructing test points, rule Bucher’s law [Reference Chen, Tang, Luo and Wang27] should be obeyed. Assuming each influencing factor in 
$\Theta $ is normally distributed and there is 
${x_i}\sim({u_i},{\sigma _i})$, the test point is constructed according to the following rules:
\begin{equation}\left\{\begin{array}{c}\overline{x}_1=(u_1,u_2,\cdots,u_n)\\[4pt]\overline{x}_2=(u_1+f\alpha_1,u_2,\cdots,u_n)\\[4pt]\vdots\\[4pt]\overline{x}_{n+1}=(u_1,u_2,\cdots,u_n+f\alpha_n)\\[4pt]\overline{x}_{n+2}=(u_1-f\alpha_n,u_2,\cdots,u_n)\\[4pt]\vdots\\[4pt]\overline{x}_{2n+1}=(u_1,u_2,\cdots,u_n-f\alpha_n)\end{array}\right.\end{equation}
where 
${\overline x _1}$ is the first generation central test point, therefore, 
${\overline v _1} = {[{\begin{array}{c@{\quad}c@{\quad}c}1 & {{{\overline x }_1}} & {{{\overline x }_1^2}}\end{array}}]^{\rm T}}$ and 
${\overline v _m} = {[{\begin{array}{c@{\quad}c@{\quad}c}1 & {\overline {{x_m}} } & {\overline {{x_m^2}} }\end{array}}]^{\rm T}}$.
The response surface function 
$\hat g(\Theta )$ should be as close as possible to the function at the limit state 
$g\left( \Theta \right)$. Namely, the deviation between the two functions is minimal. Thereby, the deviation of the two functions should be minimized as follows:
\begin{equation}{\min } \quad S(\kappa ) = {\boldsymbol{\kappa }} \cdot {{\textit{v}}} - {{\textit{g}}} = {\sum\limits_{j = 1}^{{\rm{2}}n + 1} {\left( {\sum\limits_{i = 1}^{2n + 1} {{\kappa _i}\overline {{v_{ij}}} - g\left(\overline {{x_j}} \right)} } \right)} ^2}\end{equation}
where the extremum condition is 
$\frac{{\partial S}}{{\partial {\kappa _i}}} = 0$, i = 1,2,…,2n + 1. The extremum of Eq. (25) can be expressed as
\begin{equation}\left\{\begin{array}{c}(V\cdot\kappa-g)^T V=0\\[8pt]V=\left[\overline{v}_ij\right]^T = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 & \overline{x}_{11}& \cdots& \overline{x}_{1n} &\overline{x}_{11}^2&\cdots&\overline{x}_{1n}^2\\[4pt]1 & \overline{x}_{21}& \cdots& \overline{x}_{2n} &\overline{x}_{21}^2&\cdots&\overline{x}_{2n}^2\\[4pt]\vdots & \vdots & & \vdots &\vdots&&\vdots\\[4pt]1 & \overline{x}_{m1}& \cdots& \overline{x}_{mn} &\overline{x}_{m1}^2&\cdots&\overline{x}_{mn}^2\end{array}\right]\end{array}\right.\end{equation}
where 
$V \in {R^{(2n + 1) \times (2n + 1)}}$ is the regression coefficient matrix.
Based on the least square principle, taking Eq. (25) as the optimization objective, Eq. (27) can be solved, and the response surface function coefficient 
$\kappa $ is
\begin{equation}\kappa = {\left({V^{\rm T}}WV\right)^{ - 1}}{V^{\rm T}}Wg\end{equation}
where g = (
$g({\overline x _1})$,…,
$g({\overline x _m})$)T, 
$W \in {R^{m \times m}}$ is the weight matrix, W = diag(w 1…w m), and 
${w_i} = \dfrac{{\mathop {\min }\limits_{i = 1}^m \left| {g({{\overline x }_i})} \right|}}{{\left| {g({{\overline x }_i})} \right|}}$.
Because of the limitation on the test points, the response surface function cannot approach the limit state function completely in a single time. Therefore, the corresponding convergence criteria are necessary to be established in the process of solving the response surface function. The above issue can be solved by the multiple iteration method.
Flowchart of solution of analytic expression of response surface function.
The design point of the s generation is assumed as 
$\left( {\overline x _D^{(s)},\hat g(\overline x _D^{(s)})} \right)$. The construction method of 
$\overline x _D^{(s)}$ is shown in Eq. (24). 
$\overline x _D^{(s)}$ conforms to the normal distribution (
${\mu _i},{\sigma _i}$).
\begin{equation}\left\{ {\begin{array}{c}\overline x _D^{(s{\rm{ + 1}})} = {\mu _i}{\rm{ + }}{\sigma _i}{\lambda _i}^{(s)}{\beta ^{(s)}}\\[4pt]{\beta ^{(s)}} = {\mu _{{{\hat g}^{(s)}}}}/{\sigma _{{{\hat g}^{(s)}}}}\\[4pt]{\lambda _i^{(s)} = - {{\left( {\dfrac{{\partial {{\hat g}^{(s)}}}}{{\partial {x_i}}}} \right)}^2} \cdot {\sigma _i}/{{\left( {\sum\limits_{i = 1}^n {{{\left( {\dfrac{{\partial {{\hat g}^{(s)}}}}{{\partial {x_i}}}} \right)}^2} \cdot {\sigma _i}^2} } \right)}^{1/2}},i = 1 \sim n}\end{array}} \right.\end{equation}
where 
${\beta ^{(s)}}$ is the reliability index of the s th generation, and it is the ratio of expectation and standard deviation of s th generation limit state function.
Because the design point satisfies the limit state equation, 
$g\left(\overline x _D^{(s)}\right) = 0$, 
${\beta ^{(s)}}$ can be solved by Eq. (29).
\begin{equation}g\left({u_i} + {\sigma _i}\lambda _i^{(s)}{\beta ^{(s)}}\right) = 0\end{equation}
${\beta ^{(s)}}$ can be substituted into Eq. (30) to obtain the 
$s + 1$ generation test center point, as expressed in Eq. (30) [Reference Jia, Li, Chen, Sun and Zhang6].
\begin{equation}{\overline x ^{(s + 1)}} = \mu + \left( {\overline x _D^{(s + 1)} - \mu } \right)\frac{{\hat g(\mu )}}{{\hat g(\mu ) - \hat g\left(\overline x _D^{(s + 1)}\right)}}\end{equation}
If the initial design point is defined as 
$\left( {\overline x _D^{(1)},\hat g(\overline x _D^{(1)}} \right)$, 
$\overline x _D^{(1)} = \mu $ = (
${\mu _1}$,
${\mu _2}$,…,
${\mu _n}$) can be obtained. The mathematical expression and reliability index of the response surface function can be obtained using the initial design points and test points. On the basis of the above, the next generation of the design and test points can also be obtained. During the iterative loop, when 
$\left| {{\beta ^{(s + 1)}} - {\beta ^{(s)}}} \right| < \xi $ is satisfied, convergence is considered to be achieved, and the loop ends. So the difference of reliability index 
${\beta ^{(s)}}$ between adjacent moments is less than 
$\xi $, which is the convergence criterion for solving the coefficient process. At this point, 
$g\left( \cdot \right) = \hat g\left( \cdot \right)$, and the calculation process is shown in Fig. 5. Due to the high computational efficiency of the nonresponse surface method in dealing with the limit state of complex implicit functions, the nonresponse surface method is adopted to solve the sensitivity quickly and avoid the complex multidimensional integral continuous partial derivative problem.
DH parameters of UR10 manipulator.
Diagram of UR robot kinematic coordinate systems.
4. Sensitivity simulation
4.1. Kinematic analysis of manipulator
Kinematics analysis is the basis of trajectory planning. When a robot is working, it is necessary to obtain the displacement of each joint based on the position and posture of the tool at the end of the robot to achieve the specified position. The kinematic coordinate system of the UR10 manipulator was established according to the rules of the DH method, as shown in Fig. 6. The DH parameters are listed in Table I.
According to the given DH parameters, the transformation matrix,
${}_1^0T$, of the first link can be obtained as follows:
\begin{equation}{}_1^0T = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\cos {\theta _i}} & { - \sin {\theta _i}} & 0 & 0\\[4pt]{\sin {\theta _i}} & {\cos {\theta _i}} & 0 & 0\\[4pt]0 & 0 & 0 & {{d_i}}\\[4pt]0 & 0 & 0 & 1\end{array}} \right]\end{equation}
${}_2^1T$, 
${}_3^2T$, 
${}_4^3T$, 
${}_5^4T$, 
${}_6^5T$ can be obtained similarly to 
${}_1^0T$. The transformation matrix of each link can be multiplied to obtain the transformation matrix between the end and base of the manipulator. The base coordinates 
${}_6^0T$ can be obtained by multiplying the transformation matrix of each link in turn.
\begin{equation}{}_6^0T= {}_1^0T\ {}_2^1T\ {}_3^2T\ {}_4^3T\ {}_5^4T\ {}_6^5T = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{n_x}} & {{o_x}} & {{a_x}} & {{p_x}}\\[4pt]{{n_y}} & {{o_y}} & {{a_y}} & {{p_y}}\\[4pt]{{n_z}} & {{o_z}} & {{a_z}} & {{p_z}}\\[4pt]0 & 0 & 0 & 1\end{array}} \right]\end{equation}
4.2. Sensitivity simulation of influencing factors
Based on the multilayer mapping model, the sensitivity of each influencing factor on the motion reliability of UR10 was analyzed. Based on the sensitivity analysis, the factors affecting the sensitivity of the robot end position are determined, which is convenient to increase the priority of the suppression of sensitive factors in the control and provides the basis for the trajectory planning and precision control of the robot.
4.2.1. DH kinematic parameters
First, the sensitivity of the joint angles to the motion reliability was investigated. The ideal joint angle sequence is assumed as 
$\overline \theta = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{{15}^\circ }} & {{{15}^\circ }} & {{{15}^\circ }} & {{{15}^\circ }} & {{{15}^\circ }} & {{{15}^\circ }}\end{array}} \right]$. Based on the response surface method, the constructed limit state function is expressed in Eq. (33), where 
$\Theta = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{\theta _1}} & {{\theta _2}} & {{\theta _3}} & {{\theta _4}} & {{\theta _5}} & {{\theta _6}}\end{array}} \right]$,
${\theta _i} \,:\, \left( {{\mu _i},{\sigma _i}} \right)$, and 
$\mu = \overline \theta $. Assuming that the standard deviations of all the joint angles are the same, that is, 
${\sigma _i} = 0.01$, first difference coefficient 
$f = 3$ and subsequent difference coefficient 
$f = 1$ are defined. Fifteen groups of test points are included, and the motion reliability response surface function under the influence of the joint angles can be obtained as follows:
\begin{align}\widehat g\left( \theta \right) = 1.4859 - 0.0034{\theta _1} + 0.0075{\theta _2} + 0.0074{\theta _3} + 0.0076{\theta _4} + 0.0190{\theta _5} + 0.0062{\theta _6}\nonumber\\[4pt] + {10^{ - 5}}\left( {6.1610\theta _1^2 + 36.108\theta _2^2 + 34.462\theta _3^2 + 35.911\theta _4^2 - 7.3552\theta _5^2 + 29.076\theta _6^2} \right) \end{align}
In order to facilitate the comparison of the sensitivity of influencing factors to reliability, the sensitivity results are normalized. The normalization method is shown in Eq. (34). All influence factor analysis results are normalized in this study. In each set of parameter sensitivity analysis, there are always two sets of sensitivity 0 and 1.
\begin{equation}{\varsigma ^ * } = \frac{{\left( {\varsigma - {\varsigma _{\min }}} \right)}}{{\left( {{\varsigma _{\max }} - {\varsigma _{\min }}} \right)}}\end{equation}
where 
${\varsigma ^ * }$ is the sensitivity after normalization, 
${\varsigma _{\min }}$ is the minimum sensitivity, and 
${\varsigma _{\max }}$ is the maximum sensitivity.
The sensitivity analysis of the geometric parameters of the mechanical arm to the motion reliability was conducted similar to that of joint angles, in which the distribution standard deviation of 
$\alpha $ was 0.01 rad and that of 
$d$ and 
$a$ was 0.01 mm. The analysis results of the sensitivity of each joint angle to the motion reliability can be obtained, as shown in Fig. 7. The influence of different joint parameters on motion reliability is different under the same DH parameter error. As shown in Fig. 7(a), the joint angle error of joint 4 has the highest sensitivity to motion reliability, and joint 1 has the lowest sensitivity. When controlling the robot, the angle error of joint 4 should be controlled first to improve the motion reliability of the robotic arm. After normalization, the sensitivity of joint 1 angle error is 0, but this parameter still has certain influence on motion reliability. As shown in Fig. 7(b), the rod twist angle error of joint 1 has the highest sensitivity to motion reliability, and joint 4 has the lowest sensitivity. As shown in Fig. 7(c), the rod length error of joint 6 has the highest sensitivity to motion reliability, and joint 5 has the lowest sensitivity. As shown in Fig. 7(c), the rod length error of joint 6 has the highest sensitivity to motion reliability, and joint 5 has the lowest sensitivity. As shown in Fig. 7(d), the rod offset distance error of joint 6 has the highest sensitivity to motion reliability, and joint 1 has the lowest sensitivity.
Sensitivity of DH parameter errors to motion reliability.
The results present that the sensitivity of the geometrical parameters of the manipulator and the parameters of each link is different. Even though the geometric parameters are not controllable in the movement of the manipulator, the parameter errors with high sensitivity can be suppressed preferentially by parameter calibration and other methods, which can, in turn, improve the motion reliability.
4.2.2. Link flexibility factor
When the sensitivity of the link flexibility factor of UR10 to the motion reliability is analyzed, the influencing factor vector is 
$\Theta = \left( {{\varphi _{ix}},{\varphi _{iy}},{\varphi _{iz}},{\varepsilon _{ix}},{\varepsilon _y},{\varepsilon _z}} \right)$. Ideally, the links are rigid, that is, 
$\Theta = \left( {0,0,0,0,0,0} \right)$. During the actual movement, the links will generate flexible deformation, and all the deformations obey the normal distribution of (0,0.1). The sensitivity of the link bending variable and link torsion deformation to the motion reliability are displayed in Fig. 8, respectively.
Sensitivity of arm flexibility factor to motion reliability.
The results show that the motion reliability sensitivity of all the flexible factors varies significantly with the same error of each deformation factor. Joint 6 has the greatest motion reliability sensitivity on all three bending factors. For three torsion factors, joint 1, joint 4, and joint 3 have the highest motion reliability sensitivity.
4.2.3. Joint flexibility factor
As the joint flexibility variable is subject to normal distribution 
${\xi _i} \in (0,0.01)$, its sensitivity to the motion reliability can be calculated, as shown in Fig. 9. The sensitivity of the joint flexibility to the motion reliability in descending order is joint 4, joint 3, joint 1, joint 2, joint 5, and joint 6.
Sensitivity of joint flexibility to motion reliability.
4.2.4. Joint clearance
For UR10, the influencing factors of the six joint clearances are defined as 
$\Theta = \left( \Delta {x_i},\Delta {y_i},\right.$ 
$\left.\Delta {z_i},\Delta {\theta _{xi}},\Delta {\theta _{yi}},\Delta {\theta _{zi}} \right)$. Ideally, the joint clearance of each joint is 
$\Theta = \left( {0,0,0,0,0,0} \right)$. In the actual movement process, each joint will generate a joint clearance. It is assumed that the influencing factors of each joint clearance will be generated, and the movement and rotation errors are subject to the normal distributions of 
$\left( {0,0.01} \right){\rm{mm}}$ and 
$\left( {0,0.01} \right){\rm{rad}}$, respectively.
The sensitivity of the moving and rotating variables to the motion reliability of each joint is displayed in Fig. 10, respectively. The results show that the sensitivity of joint clearance factors to motion reliability is basically consistent with link flexibility factors. The reason is that these two factors have the same effect on robot kinematics.
Sensitivity of joint clearance factors to motion reliability.
5. Fault diagnosis based on KPCA method
5.1. Introduction of KPCA method
In the traditional principal component analysis (PCA) method, the monitoring process space can be divided into the principal component subspace and the residual subspace, and the fault detection and isolation are realized in the two subspaces separately.
In the preceding sections, the sensitivity relationships are expressed linearly depending on the deviations of the involved entities from their nominal values. However, the purpose of sensitivity analysis is to analyze the influence between error of end-effector and influencing factors. The relationship between error of end-effector and influencing factors is nonlinear. There are numerous nonlinear processes in the operation of a manipulator, particularly with a large range of input values. In this case, many process information and such nonlinear relations cannot be described by the PCA method. However, the KPCA can be used to deal with such nonlinear problems. The KPCA method is an extension of the nonlinear form of the PCA method. In this study, the Gaussian kernel function, 
$k(x,y) = \exp \left( {\frac{{ - {{\left\| {x - y} \right\|}^2}}}{{2{\sigma ^2}}}} \right)$, is used to replace the original data for the calculation, and the process is shown in Fig. 11.
Fault diagnosis flowchart of KPCA method.
A mixed indicator 
$\varphi $ is introduced to comprehensively reflect the changes in the 
${T^2}$ and 
$SPE$ statistics [Reference Yue and Qin28], and its calculation formula is as follows:
\begin{equation}\varphi = c\frac{{SPE}}{{{\delta ^2}}} + (1 - c)\frac{{{T^2}}}{{\chi _1^2}}\end{equation}
In the formula, 
${\delta ^2}$ is the control limit of the 
$SPE$ and 
$\chi _1^2$ is the control limit of 
${T^2}$. T 2 term, Hotelling’s T2, is used to detect whether the principal subspace of data is abnormal. For different purposes, 
$c$ has different selection criteria. Here, the complementarity of and is to be utilized. Therefore, their weights of the two can be the same, that is, c = 0.5.
5.2. Simulation verification of fault diagnosis
In this section, KPCA fault diagnosis algorithm that was applied to verify the DH parameter sensitivity of the directly affected layer is described. The sensitivity verification process for joint angular error 
$\theta $ is as follows:
Step 1: The training and test data were constructed using MATLAB software. In the training data, 600 groups of joint angle values (
$\theta _{1 - 6}^i$,i = 1,2,…,600) were randomly generated, which were substituted into the positive kinematic model of UR10 to obtain the corresponding end positions (${P_i}$,i = 1,2,…,600). The data sets were divided into six groups, namely ${N_1}{\rm{ \sim }}{N_6}$, with 100 samples in each group.Step 2: In the test data, the i th joints in
${N_j}\left( {j = 1 - 6} \right)$ were given the same error E from the angle values of the 51–100 samples, and subsequently, they were substituted into the positive kinematic model of the manipulator to obtain the comprehensive end position (${P^{\prime}_i},i = 1{\rm{ \sim }}600$). The data sets were also divided into six groups, namely ${N^{\prime}_1} {\rm{ \sim }}{N^{\prime}_6}$.Step 3: The training and test data were introduced into the KPCA fault diagnosis program to obtain the comprehensive changes of
${T^2}$ and the $SPE$ statistics of the joint rotation errors.
The comprehensive changes of 
${T^2}$ and 
$SPE$ statistics of the joint rotation errors are shown in Fig. 12. Both 
${T^2}$ and 
$SPE$ present great fluctuations and changes from the 50th sample data. The PCA and KPCA methods are used to analyze the contribution rate of faults. The mixed contribution rate of the joint 3 angle error is shown in Fig. 13. The results show that the two methods can achieve consistent results.
Change of fault detection index.
Joint 3 angle fault diagnosis result. (a) Changes of 
${T^2}$. (b) Changes of 
$SPE$ statistics.
The comprehensive contribution of the DH parameter is shown in Fig. 13. Both 
${T^2}$ and the 
$SPE$ have significant fluctuations and variations due to the 50th sample data. By comparing the contribution rates of all the joint errors, it can be found that joint 3 has the largest comprehensive contribution rate, which corresponds to joint 5 in Fig. 7(a) with the largest sensitivity.
5.3. Experimental verification of joint angular fault diagnosis
The validity of the KPCA fault diagnosis method was verified in this section by experiments. A camera (Optitrack V120:Trio) was used to measure and collect the pose data of the UR10 manipulator end. The experimental scheme is shown in Fig. 14. The intelligent camera and the UR10 manipulator are communicated with the computer through transmission control protocol/internet protocol, and the joint rotation angles of UR10 manipulator were controlled by the computer. Five marked points were installed at the end of the manipulator. When the manipulator reached a new position, the camera collects the end-effector pose data automatically.
Experimental scheme of joint angle fault diagnosis.
5.3.1. Training data (no fault) are built
Within the visual range of the camera, UR10 was controlled by the computer to move randomly to 100 different positions, and the six joint angles were recorded. Simultaneously, the end positions of manipulator were collected by the camera. The above steps were repeated six times, and 600 sets of training data corresponding to the six joint angles under no fault condition were obtained.
5.3.2. Test data (with fault) are generated
For the six joint angles, the test data were generated, respectively. Joint 1 was taken as an example, the first 50 joint angle data and the end-effector pose data were constant, and the data are the same as the first 50 groups in the training data. The last 50 joint angle data corresponding to joint 1 were, respectively, applied by 1°. The new joint angles were controlled by the computer, and the new positions were collected, respectively. From this, 50 sets of test data corresponding to joint 1 were obtained with faults. The test data of the other five joints were constructed similarly. Thus, 100 sets of fault test data corresponding to the six joint angles were constructed successfully.
5.3.3. Fault diagnosis
The training and test data were input into the KPCA fault diagnosis program, and the fault diagnosis results are displayed in Fig. 15. When there was a fault on a joint, KPCA could produce the correct diagnosis, and the comprehensive contribution rate of the corresponding joint angle was the highest.
Joint angle fault diagnosis result.
6. Conclusion
In this study, a six-DOF industrial robot was applied as the research object, and the influencing factor mapping model of its motion accuracy reliability was established comprehensively considering various factors. Moreover, the reliability of the end-effector pose was solved based on the response surface method. In addition, the sensitivity was verified from the data angle, and the joint rotation error was diagnosed in terms of the KPCA fault diagnosis principle. The obtained conclusions are obtained as follows:
(1) An extensible mapping model for the influencing factors of the motion reliability of the manipulator was developed. When there are new influencing factors, it can be directly extended on the basis of the existing model.
-
(2) Based on the sensitivity analysis of the motion reliability of this model, the factors with a large controllable influence should be considered as the priority control target or assigned a large weight. The uncontrollable factors should be suppressed by an indirect control strategy.
-
(3) The fault diagnosis based on data can not only compensate for the deficiency of empirical model diagnosis but also play a role in identifying new influencing factors.
In this study, only single-joint fault diagnosis is studied. In the future, multijoint fault diagnosis based on deep learning method will be studied. Effective diagnosis can be made when the manipulator has multiple faults.
Funding Support
This work was supported by the National Key R&D Program of China(2018YFB1308100), Zhejiang Provincial Natural Science Foundation under Grant LQ21F020026, Science foundation of Zhejiang Sci-Tech University(ZSTU) under Grant No. 19022104-Y, General Scientific Research Project of Zhejiang Provincial Department of Education (Grant Nos.19020038-F, 19020033-F, Y202045464), National Natural Science Foundation (Grant 51805488), and Beijing Satellite Environmental Engineering Research Institute (CAST-BISEE2019-001).
Conflict of interest
The authors have no conflicts of interest to declare.
