A hybrid least-squares genetic algorithm–based algorithm for simultaneous identification of geometric and compliance errors in industrial robots

Abstract

Due to the flexibility of robot joints and links, industrial robots can hardly achieve the accuracy required to perform tasks when a payload is attached at their end-effectors. This article presents a new technique for identifying and compensating compliance errors in industrial robots. Within this technique, a comprehensive error model consisting of both geometric and compliance errors is established, where joint compliance is modeled as a piecewise linear function of joint torque to approximate the nonlinear relation between joint torque and torsional angle. A hybrid least-squares genetic algorithm–based algorithm is then developed to simultaneously identify the geometric parameters, joint compliance values, and the transition joint torques. These identified geometric and non-geometric parameters are then used to compensate geometric and joint compliance errors. Finally, the developed technique is applied to a 6 degree-of-freedom industrial serial robot (Hyundai HA006). Experimental results are presented that demonstrate the effectiveness of the identification and compensation techniques.

Introduction

Robot manipulators play an important role in the industrial field. Industrial applications such as assembly, welding, and machining operations require highly accurate robot manipulators. However, present industrial robots have high repeatability but low accuracy. Therefore, robot manipulators need to be calibrated to meet the accuracy demands of industrial applications.

Error sources that affect the robot positioning accuracy can be divided into geometric errors and non-geometric errors. ¹ Geometric errors may result from manufacturing imperfections, misalignments, and encoder offsets. Non-geometric errors may come from joint and link compliance, temperature variation, gear transmission error, and backlash in gear transmission. To the best of our knowledge, practical techniques to identify and compensate for all the geometric and non-geometric errors are not yet developed. For example, some research^2–6 focused on identifying and compensating geometric errors related to the deviation of geometrical parameters with respect to their nominal values. However, for industrial robot manipulators, compliance errors related to the deflections of joints and links caused by gravity (payload and robot links) or applied forces contribute significantly to robot inaccuracy, in addition to geometric errors. Therefore, compensating compliance errors to further enhance robot positioning accuracy has been an active area of research.

The problem of modeling and identifying compliance errors has been studied in many publications. Some works have considered both joint and link flexibilities.^7–9 However, as mentioned in Elatta et al., ¹ Becquet, ¹⁰ and Gong et al., ¹¹ compliance errors due to joint deflections are much more significant than those due to link deflections. Other research has focused on modeling and identifying joint compliance to compensate compliance errors due to joint deflections^11–21 under the assumption that robot links are much stiffer than robot joints. For simplicity, robot joints in these works are mainly modeled as a linear torsional spring; joint deflections are the result of multiplication of constant joint compliance (inverse of joint stiffness) and joint torques. However, the experimental study in Kircanski et al. ²² reveals that the joint torque–torsion relation is nonlinear; joint torque can be expressed by a third-order polynomial function of joint deflections. Not much attention has been paid to identification and compensation of compliance errors related to nonlinear joint compliance. Jang et al. ²³ divided the robot workspace into several local regions. In each local region, joint angular errors were modeled as the sum of geometric error and joint deflections, and then calibrated. A radial basis function network was then used to obtain the continuous error function of joint angular errors at any specified joint coordinate space. Similarly, Meggiolaro et al. ²⁴ defined compliance errors (elastic errors) and geometric errors in a unified manner; they are approximated by a polynomial function of robot configurations and applied wrench. As these two methods used nonlinear error models, compliance errors related to nonlinear joint stiffness may be considered. However, these two methods require a large number of measurements in the calibration process, which is time consuming. In addition, the effect of geometric errors and compliance errors on the robot positioning accuracy cannot be individually observed in these two methods.

This article is devoted to identifying and compensating both geometric errors and compliance errors in industrial robots. To achieve this, a comprehensive error model is derived for combining geometric errors and compliance errors. In contrast to previous works, this error model takes into account the nonlinear relation between joint torque and joint deflection presented in Kircanski et al. ²² by modeling joint stiffness as a piecewise linear function of joint torque. For each joint of the robot manipulator, three non-geometric parameters must be identified in addition to geometric parameters: joint stiffness (denoted by $k_l$ and $k_h$ ) and transition joint torque (denoted by ${\tau }_t$ ). Because there is no direct relation between the transition joint torque and the position errors of the end-effector of a robot manipulator, the proposed comprehensive error model is established under the assumption that transition joint torques are given. Based on the comprehensive error model with the given transition joint torques, an iterative least-squares algorithm is used to identify the geometric parameters and joint stiffness values. For identification of the transition joint torque, the genetic algorithm (GA) is adopted for its advantages such as the ability to explore many regions of search space simultaneously, lower sensitivity to local minima, and not requiring the calculation of derivatives. By combing GA and the least-squares algorithm (i.e. the least-squares algorithm for identifying geometric parameters and joint stiffness values is considered as the fitness function of GA), the defined geometric and non-geometric parameters can be identified simultaneously. To demonstrate the advantages of the developed techniques, experimental calibration for a 6-degree-of-freedom (dof) serial robot (Hyundai HA006) was carried out. The results show that the position accuracy of the robot is significantly improved at a rate of 98% after compensating both geometric and compliance errors based on the identified parameters.

The remainder of this article is organized as follows. Sections “Modeling of kinematic error” and “Modeling of compliance error” describe the modeling of kinematic errors and compliance errors, respectively. Section “The proposed hybrid algorithm for robot calibration” presents the proposed hybrid algorithm for identifying these errors. Section “Experimental setup and measurements” describes the experiment design. Section “Experimental calibration results and validation” presents experimental calibration results for the Hyundai robot HA006. In addition, the calibration results are validated. The final section outlines some conclusions.

Modeling of kinematic error

The Denavit–Hartenberg (D-H) method ²⁵ is used to model the kinematics of the Hyundai robot HA006, which allows one to model the robot joint with four parameters: joint angle ( ${\theta }_i$ ), offset distance ( $d_i$ ), link length ( $a_{i-1}$ ), and twist angle ( ${\alpha }_{i-1}$ ). Based on D-H method, the coordinate frames for each link of the Hyundai robot HA006 are assigned as shown in Figure 1; its nominal D-H parameters are derived and listed in Table 1. The relative transformation (transition and rotation) between link frame $\{i-1\}$ and link frame $\{i\}$ can be described by a homogenous transformation matrix as follows

$${}_i^{i-1}\mathit{T}=\mathrm{Rot}(x_{i-1},{\alpha }_{i-1})\mathrm{Tr}(x_{i-1},a_{i-1})\mathrm{Tr}(z_i,d_i)\mathrm{Rot}(z_i,{\theta }_i)$$

(1)

Figure 1.

A sketch of Hyundai robot HA006 and its attached link frames.

Table 1.

Nominal D-H parameters of Hyundai robot HA006.

i	${\alpha }_{i-1}$ (°)	$a_{i-1}$ (m)	$d_i$ (m)	${\theta }_i$ (°)
1	0	0	0.36	${\theta }_1$
2	90	0.2	0	${\theta }_2$
3	0	0.56	0	${\theta }_3$
4	90	0.13	0.62	${\theta }_4$
5	−90	0	0	${\theta }_5$
6	90	0	0.1	${\theta }_6$

In the case of two consecutive parallel or near parallel joints (see joints 2 and 3 in Figure 1), small errors in the end-effector position cannot be modeled by small errors in the D-H link parameters because equation (1) is singular, so a small rotation of $\beta$ about the y-axis is added. ²⁶ Specifically, the homogenous transformation matrix ${}_3^2\mathit{T}$ can be defined as follows

$${}_3^2\mathit{T}=\mathrm{Rot}(x_2,{\alpha }_2)\mathrm{Tr}(x_2,a_2)\mathrm{Rot}(y_2,{\beta }_2)\mathrm{Tr}(z_3,d_3)\mathrm{Rot}(z_3,{\theta }_3)$$

(2)

Generally, the world fixed frame is defined arbitrarily by the user. Thus, the robot base frame $\{0\}$ can be defined by six parameters $({\alpha }_S,a_S,{\beta }_S,b_S,{\theta }_0,d_0)$ with respect to the laser sensor frame $\{S\}$ . The transformation ${}_0^S\mathit{T}$ is described as follows

$$\begin{array}{cc}\multicolumn{1}{c}{{}_0^S\mathit{T}=}& \mathrm{Rot}(x_S,{\alpha }_S)\mathrm{Tr}(x_S,a_S)\mathrm{Rot}(y_S,{\beta }_S)\hfill \\ \multicolumn{1}{c}{}& \mathrm{Tr}(y_S,b_S)\mathrm{Tr}(z_0,d_0)\mathrm{Rot}(z_0,{\theta }_0)\hfill \end{array}$$

(3)

Since ${\alpha }_0=0$ , $a_0=0$ , the transformation ${}_1^S\mathit{T}$ can be given as follows

$$\begin{array}{cc}\multicolumn{1}{c}{{}_1^S\mathit{T}=}& \mathrm{Rot}(x_S,{\alpha }_S)\mathrm{Tr}(x_S,a_S)\mathrm{Rot}(y_S,{\beta }_S)\hfill \\ \multicolumn{1}{c}{}& \mathrm{Tr}(y_S,b_S)\mathrm{Tr}(z_1,d{\prime }_1)\mathrm{Rot}(z_1,\theta {\prime }_1)\hfill \end{array}$$

(4)

where $d{\prime }_1=d_0+d_1$ and $\theta {\prime }_1={\theta }_0+{\theta }_1$ . As the laser tracker (see section “Experimental setup and measurements”) only measures the position data, the translation between frame $\{6\}$ and that of the end-effector (reflector) $\{E\}$ can be described as follows

$${}_E^6\mathit{T}=\mathrm{Tr}(x_6,a_E)\mathrm{Tr}(y_6,b_E)\mathrm{Tr}(z_6,d_E)$$

(5)

where the vector $[a_E,b_E,d_E]$ is the location of reflector with respect to the frame $\{6\}$ .

By successive multiplications of the transformation matrices, the position ( ${}_E^S\mathit{P}$ ) and orientation ( ${}_E^S\mathit{R}$ ) of the reflector with respect to the sensor base frame $\{S\}$ can be determined as follows

$${}_E^S\mathit{T}={}_1^S\mathit{T}{}_2^1\mathit{T}{}_3^2\mathit{T}{}_4^3\mathit{T}{}_5^4\mathit{T}{}_6^5\mathit{T}{}_E^6\mathit{T}=\left[\begin{array}{cc}\hfill {}_E^S\mathit{R}\hfill & \hfill {}_E^S\mathit{P}\hfill \\ \multicolumn{1}{c}{\mathit{O}}& \hfill 1\hfill \end{array}\right]$$

(6)

where $\mathit{O}$ is a $(3\times 1)$ zero vector.

The resultant position errors of the end-effector ( $\Delta {\mathit{X}}_{\mathit{kin}}$ ) due to small errors in kinematic parameters (denoted by $\Delta \varnothing =[\Delta \mathit{\alpha },\Delta \mathit{a},\Delta \mathbf{\beta },\Delta \mathit{b},\Delta \mathit{d},\Delta {\mathbf{\theta }}_{\mathit{o}}](n\times 1)$ ) can be obtained by differentiating the homogenous transformation ${}_E^S\mathit{T}$ of equation (6) and is described as follows

$$\Delta {\mathit{X}}_{\mathit{kin}}={\mathit{J}}_{\mathit{kin}}\Delta \varnothing$$

(7)

where ${\mathit{J}}_{\mathrm{kin}}=[{\mathit{J}}_{\alpha }:{\mathit{J}}_a:{\mathit{J}}_{\beta }:{\mathit{J}}_b:{\mathit{J}}_d:{\mathit{J}}_{\theta }](3\times n)$ is the extended kinematic Jacobian matrix and $n$ is the number of kinematic parameters that will be calibrated. Each column in the matrix ${\mathit{J}}_{\mathit{kin}}$ corresponds to each parameter error in $\Delta \varnothing$ and can be calculated directly as follows ¹⁴

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mathit{J}}_{{\alpha }_{i-1}}=[{\mathit{x}}_{i-1}\times {\mathit{L}}_{i-1}],{\mathit{J}}_{a_{i-1}}=\left[{\mathit{x}}_{i-1}\right],\hfill \\ \multicolumn{1}{c}{}& {\mathit{J}}_{{\beta }_{i-1}}=[{\mathit{y}}_{i-1}\times {\mathit{L}}_{i-1}],{\mathit{J}}_{b_{i-1}}=\left[{\mathit{y}}_{i-1}\right],\hfill \\ \multicolumn{1}{c}{}& {\mathit{J}}_{{\theta }_i}=[{\mathit{z}}_i\times {\mathit{L}}_i],{\mathit{J}}_{d_i}=\left[{\mathit{z}}_i\right]\hfill \end{array}$$

(8)

where ${\mathit{x}}_i$ , ${\mathit{y}}_i$ , and ${\mathit{z}}_i$ are $(3\times 1)$ directional vectors of the link frame $\{i\}$ with respect to frame $\{S\}$ , and $L_i$ is a $(3\times 1)$ vector between the origin of frame $\{i\}$ and $\{E\}$ with respect to frame $\{S\}$ .

There are 32 kinematic parameters in equation (6); however, some kinematic parameters are dependent on others: $\{{\theta }_1,{\theta }_0\}$ , $\{d_1,d_0\}$ , $\{d_2,d_3\}$ , $\{d_E,d_6\}$ , and $\{(a_E,b_E),{\theta }_6\}$ . These dependent parameters cannot be calibrated at the same time. Therefore, the first parameter in each pair is selected and calibrated. In addition, the error of the not-selected parameter is forced to zero. Thus $n$ is equal to 27.

Modeling of compliance error

This section first presents a new nonlinear joint stiffness model that allows us to take into account the nonlinear relation between joint torque and torsion. Based on this joint stiffness model, a mathematical formulation is then derived to describe the related compliance error.

Nonlinear joint stiffness modeling

The experimental study in Kircanski et al. ²² reveals that the behavior of the robot joint is nonlinear; the joint stiffness coefficient can be modeled as a second-order polynomial in joint torsion. Based on this nonlinear joint stiffness model, the torque–torsion relation is defined as ²²

$$\tau \left(q\right)=a_{f1}{q}^3+a_{f2}q+K_{\mathrm{sw}}\left(q\right)q$$

(9)

where $a_{f1}$ and $a_{f2}$ are coefficients; $K_{\mathrm{sw}}(q)$ is the correction factor of the soft windup that causes lower joint stiffness at lower torques. For more detailed information, refer to Kircanski et al. ²² A linearization of this nonlinear torque–torsion function has been used in motion control design. ²⁷

Simpler joint stiffness models of the robot correlate to easier identification of its stiffness parameters. We model the robot joint as a piecewise linear system to approximate its nonlinear behavior, which allows us to obtain a linear form of the identification equations that can be easily solved using a standard least-squares technique (see section “The proposed hybrid algorithm for robot calibration”). This linearization is conceptually explained in Figure 2. As the torque applied to the joint increases, the joint torsion increases at a diminishing rate. Accordingly, the joint stiffness as a function of applied torque can be approximated as

$$k=\{\begin{array}{c}\hfill k_l\hfill \\ \multicolumn{1}{c}{k_h}\end{array}\begin{array}{c}\hfill \left|\tau \right|<{\tau }_t\hfill \\ \multicolumn{1}{c}{\left|\tau \right|\ge {\tau }_t}\end{array}$$

(10)

Figure 2.

Conceptual drawing for linearization of the nonlinear torque–torsion relation.

where ${\tau }_{\mathrm{t}}$ is the transition torque; $k_1$ and $k_h$ are joint stiffness values ( $k_1<k_h$ ). These three parameters and kinematic parameters must be calibrated for each joint.

Compliance error model

Assuming the robot’s links are much stiffer than its joints, compliance errors are due to the flexibility of robot joints under the link self-gravity and external payload. Effective torque in joint i is due not only to the gravity force related to link $i$ but also to the gravity force related to the subsequent links and the payload. Thus, the total effective torque in joint $i$ is given as

$${\tau }_i={\displaystyle \sum _{j=1}^{N+1}}{\tau }_{i,j}={\displaystyle \sum _{j=i}^{N+1}}{\mathit{J}}_{{\theta }_{i,j}}^T{\mathit{F}}_{\mathit{j}}$$

(11)

where ${\mathit{F}}_{\mathit{j}}=[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -M_{j}g\hfill \end{array}{]}^{\mathrm{T}}$ is a gravity force vector, $M_j$ is the mass of link $j$ , $g$ is the gravity coefficient, $N$ is the number of dof of the robot, and ${\mathit{F}}_{N+1}$ is the gravity force vector related to the payload. It is worth mentioning that ${\tau }_i$ should be computed at the deflected equilibrium position, which is explained in section “The proposed hybrid algorithm for robot calibration.” The transpose of the Jacobian matrix ${\mathit{J}}_{{\theta }_{i,j}}$ is used as a force transformation to find the effective joint torque ${\tau }_{i,j}$ in joint $i$ due to the gravity force ${\mathit{F}}_{\mathit{j}}$ . Based on the definition of the Jacobian matrix, ${\mathit{J}}_{{\mathbf{\theta }}_{\mathit{i},j}}$ is given as

$${\mathit{J}}_{{\mathbf{\theta }}_{\mathit{i},j}}={\mathit{z}}_{\mathit{i}}\times {\mathit{l}}_{\mathit{i},j}$$

(12)

where ${\mathit{l}}_{\mathit{i},j}$ is the $3\times 1$ vector between the origin of the frame $\{i\}$ and the mass center of link $j$ .

Based on the nonlinear joint stiffness model presented in section “Nonlinear Joint stiffness modeling,” the small joint deflection in joint $i$ due to the applied torque ${\tau }_i$ can be calculated as follows

$$\Delta {\theta }_c^i=\{\begin{array}{c}\hfill {\tau }_i/k_l^i\hfill \\ \multicolumn{1}{c}{{\tau }_i/k_h^i}\end{array}\begin{array}{c}\hfill \left|{\tau }_i\right|<{\tau }_t^i\hfill \\ \multicolumn{1}{c}{\left|{\tau }_i\right|\ge {\tau }_t^i}\end{array}$$

(13)

The Cartesian position errors of the robot end-effector due to the small joint deflections $\Delta {\theta }_c$ are given as

$$\Delta {\mathit{X}}_{\mathit{c}}={\mathit{J}}_{\mathbf{\theta }}\Delta {\mathbf{\theta }}_{\mathit{c}}={\displaystyle \sum _i^N}\Delta {\mathit{X}}_c^i={\displaystyle \sum _i^N}{\mathit{J}}_{{\theta }_i}\Delta {\theta }_c^i$$

(14)

where $\Delta {\mathit{X}}_{\mathit{c}}^{\mathit{i}}$ is the compliance errors due to the deflection in joint $i$ , ${\mathit{J}}_{\mathbf{\theta }}=[{\mathit{J}}_{{\mathbf{\theta }}_1},\dots ,{\mathit{J}}_{{\mathbf{\theta }}_N}]$ is the standard Jacobian matrix, and $\Delta {\mathbf{\theta }}_{\mathit{c}}=[\Delta {\theta }_c^1,\dots ,\Delta {\theta }_c^N{]}^{\mathrm{T}}$ .

It is the fact that z-axes (the direction of joint rotation) of joints 2 and 3 lie within the gravity force flow and are coplanar with the gravity forces acting on the succeeding joint 4, 5, and 6 which leads to joints 2 and 3 suffering larger torques. Therefore, the compliance errors related to joints 2 and 3 are the most significant. From equations (13) to (14), the compliance errors due to the deflections in joints 2 and 3 can be modeled as

$$\begin{array}{cc}\multicolumn{1}{c}{}& \Delta {\mathit{X}}_{\mathit{c}}\cong \Delta {\mathit{X}}_{\mathit{c}}^2+\Delta {\mathit{X}}_{\mathit{c}}^3\hfill \\ \multicolumn{1}{c}{}& =\{\begin{array}{c}\hfill \left[\begin{array}{cc}\hfill {\mathit{J}}_{{\theta }_2}^l{\tau }_2^l\hfill & \hfill {\mathit{J}}_{{\theta }_3}^l{\tau }_3^l\hfill \end{array}\right]{\left[\begin{array}{cc}\hfill c_2^l\hfill & \hfill c_3^l\hfill \end{array}\right]}^{\mathrm{T}}\left|{\tau }_2\right|<{\tau }_t^2,\left|{\tau }_3\right|<{\tau }_t^3\hfill \\ \multicolumn{1}{c}{\left[\begin{array}{cc}\hfill {\mathit{J}}_{{\theta }_2}^h{\tau }_2^h\hfill & \hfill {\mathit{J}}_{{\theta }_3}^h{\tau }_3^h\hfill \end{array}\right]{\left[\begin{array}{cc}\hfill c_2^h\hfill & \hfill c_3^h\hfill \end{array}\right]}^{\mathrm{T}}\left|{\tau }_2\right|\ge {\tau }_t^2,\left|{\tau }_3\right|\ge {\tau }_t^3}\\ \multicolumn{1}{c}{\left[\begin{array}{cc}\hfill {\mathit{J}}_{{\theta }_2}^h{\tau }_2^h\hfill & \hfill {\mathit{J}}_{{\theta }_3}^l{\tau }_3^l\hfill \end{array}\right]{\left[\begin{array}{cc}\hfill c_2^h\hfill & \hfill c_3^l\hfill \end{array}\right]}^{\mathrm{T}}\left|{\tau }_2\right|\ge {\tau }_t^2,\left|{\tau }_3\right|<{\tau }_t^3}\\ \multicolumn{1}{c}{\left[\begin{array}{cc}\hfill {\mathit{J}}_{{\theta }_2}^l{\tau }_2^l\hfill & \hfill {\mathit{J}}_{{\theta }_3}^h{\tau }_3^h\hfill \end{array}\right]{\left[\begin{array}{cc}\hfill c_2^l\hfill & \hfill c_3^h\hfill \end{array}\right]}^{\mathrm{T}}\left|{\tau }_2\right|<{\tau }_t^2,\left|{\tau }_3\right|\ge {\tau }_t^3}\end{array}\hfill \end{array}$$

(15)

where $c_2^l$ , $c_2^h$ , $c_3^l$ , and $c_3^h$ are inversions of $k_2^l$ , $k_2^h$ , $k_3^l$ , and $k_3^h$ , respectively. Usually, parameter estimation in robot calibration requires many actual measurements over different robot configurations. For a given number of measurement configurations ( $m$ ), the effective torques ${\tau }_2$ and ${\tau }_3$ can be computed based on equations (11) and (12). These $m$ measurements are respectively distributed to each case in equation (15) under the assumption that the transition torques ${\tau }_t^2$ and ${\tau }_t^3$ are the known parameters. Accordingly, a compact compliance error model for $m$ measurements can be derived as

$$\Delta {\mathit{X}}_{\mathit{c}}\cong \left[\begin{array}{c}\hfill \begin{array}{cccc}\hfill {\mathit{J}}_{{\theta }_2}^l{\tau }_2^l\hfill & \hfill \mathit{O}\hfill & \hfill {\mathit{J}}_{{\theta }_3}^l{\tau }_3^l\hfill & \hfill \mathit{O}\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{O}& \hfill {\mathit{J}}_{{\theta }_2}^h{\tau }_2^h\hfill & \hfill O\hfill & \hfill {\mathit{J}}_{{\theta }_3}^h{\tau }_3^h\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cccc}\hfill O\hfill & \hfill {\mathit{J}}_{{\theta }_2}^h{\tau }_2^h\hfill & \hfill {\mathit{J}}_{{\theta }_3}^l{\tau }_3^l\hfill & \hfill O\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{J}}_{{\theta }_2}^l{\tau }_2^l}& \hfill O\hfill & \hfill O\hfill & \hfill {\mathit{J}}_{{\theta }_3}^h{\tau }_3^h\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \end{array}}\end{array}\right]\left[\begin{array}{c}\hfill c_2^l\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\hfill c_2^h\hfill \\ \multicolumn{1}{c}{c_3^l}\end{array}}\\ \multicolumn{1}{c}{c_3^h}\end{array}\right]={\mathit{J}}_{\mathit{c}}\mathit{C}$$

(16)

where $\mathit{O}$ is a $3\times 1$ zeros vector, ${\mathit{J}}_c$ is a $3m\times 4$ transformation matrix relating the joint compliance parameters and the deflection of the robot end-effector, and $\mathit{C}=[c_2^l,c_2^h,c_3^l,c_3^h{]}^{\mathrm{T}}$ is the joint compliance vector.

Similarly, equation (16) can be extended by including the compliance errors related to other joints. For example, if the compliance errors related to $N$ joints are considered, the dimensions of $J_c$ and $\mathbf{C}$ are $3m\times 2N$ and $2N\times 1$ , respectively.

The proposed hybrid algorithm for robot calibration

The real position vector ( ${\mathit{P}}_{\mathit{real}}$ ) of the robot end-effector can be given as

$${\mathit{P}}_{\mathit{real}}={\mathit{P}}_{\mathit{kin}}+\Delta {\mathit{X}}_{\mathit{kin}}+\Delta {\mathit{X}}_{\mathit{c}}+\Delta {\mathit{X}}_{\mathit{extra}}$$

(17)

where ${\mathit{P}}_{\mathit{kin}}$ is the results of forward kinematics with the nominal kinematic parameters; $\Delta {\mathit{X}}_{\mathit{kin}}$ , $\Delta {\mathit{X}}_{\mathit{c}}$ , and $\Delta {\mathit{X}}_{\mathit{extra}}$ are position errors due to errors in kinematic parameters, joint deflections, and errors from other sources (unmodeled errors), respectively. Under the assumption that all robot position errors come from the effect of kinematic errors and compliance errors (i.e. $\Delta {\mathit{X}}_{\mathit{extra}}=0$ ), by substituting equations (7) and (16) into equation (17), the position error model can be given as

$$\begin{array}{cc}\multicolumn{1}{c}{\Delta \mathit{X}}& ={\mathit{P}}_{\mathit{real}}-{\mathit{P}}_{\mathit{kin}}=\Delta {\mathit{X}}_{\mathit{kin}}+\Delta {\mathit{X}}_{\mathit{c}}\hfill \\ \multicolumn{1}{c}{}& =\left[\begin{array}{cc}\hfill {\mathit{J}}_{\mathit{kin}}\hfill & \hfill {\mathit{J}}_{\mathit{c}}\hfill \end{array}\right]\left[\begin{array}{c}\hfill \Delta \varnothing \hfill \\ \multicolumn{1}{c}{\mathit{C}}\end{array}\right]=\mathit{J}\Delta {\varnothing }_{\mathit{ext}}\hfill \end{array}$$

(18)

where $\mathit{J}=\left[\begin{array}{cc}\hfill {\mathit{J}}_{\mathit{kin}}\hfill & \hfill {\mathit{J}}_{\mathit{c}}\hfill \end{array}\right]$ is a $3m\times (4+n)$ matrix and $\Delta {\varnothing }_{\mathrm{ext}}={\left[\begin{array}{cc}\hfill \Delta \varnothing \hfill & \hfill C\hfill \end{array}\right]}^T$ is a $(4+n)\times 1$ vector consisting of $n$ kinematic parameter errors and four joint compliance parameters ( $c_2^l$ , $c_2^h$ , $c_3^l$ , and $c_3^h$ ). Equation (18) is established under the assumption that the transition joint torques ${\tau }_t^2$ and ${\tau }_t^3$ are given (i.e. the transition joint torques must be estimated). Once the transition joint torques are estimated, equation (18) can be solved by the least-squares approach, with the pseudo-inverse given by

$$\Delta {\varnothing }_{\mathit{ext}}=[{({\mathit{J}}^{\mathrm{T}}\mathit{J})}^{-1}{\mathit{J}}^{\mathrm{T}}]\Delta \mathit{X}$$

(19)

where $\Delta \mathit{X}={\mathit{P}}_{\mathit{M}}-{\mathit{P}}_{\mathit{kin}}$ and ${\mathit{P}}_{\mathit{M}}$ is the measured position vector of the end-effector attached with a payload. Equation (19) is applied iteratively until convergence is reached. This iterative process is shown in Figure 3(b). Initially, kinematic parameters (denoted by ${\varnothing }_{\mathrm{kin}}^0$ ) and joint deflections (denoted by $\Delta {\mathbf{\theta }}_{\mathit{c}}^0$ ) are set to nominal values and zero, respectively. ${\mathit{P}}_{\mathrm{kin}}^{i-1}$ is the results of forward kinematics with the current kinematic parameters ( ${\varnothing }_{\mathrm{kin}}^{i-1}$ ), where the kinematic values of joint angles describing the preloading robot pose are given as follows

$${\mathbf{\theta }}_{\mathrm{kin}}^{i-1}={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o^{i-1}$$

(20)

Figure 3.

Flowchart of the hybrid least-squares-GA-based algorithm: (a) GA and (b) iterative least-squares algorithm.

The values of joint angles describing the deflected equilibrium robot pose are given as follows

$${\mathbf{\theta }}_e^{i-1}={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o^{i-1}+\Delta {\mathbf{\theta }}_c^{i-1}$$

(21)

Thus, the matrices ${(J_{\mathrm{kin}})}_{\mathrm{kin}}$ and ${(J_c)}_e$ should be computed with joint angles given by equations (20) and (21), respectively. After applying equation (19), $\Delta {\varnothing }^i$ is obtained. The kinematic parameters are then updated as follows

$${\varnothing }_{\mathrm{kin}}^i={\varnothing }_{\mathrm{kin}}^{i-1}+\Delta {\varnothing }^i$$

(22)

Simultaneously, joint compliance parameters $C^i$ are obtained from equation (19). Based on these identified joint compliance parameters, another iterative process to find the equilibrium position ( ${\mathbf{\theta }}_e^i$ ), torque, and the resulting joint deflections ( $\Delta {\mathbf{\theta }}_c^i$ ) is applied as follows

1. Compute $\tau$ based on equations (11) and (12) at equilibrium joint position ${\mathbf{\theta }}_e^i={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o^i+{(\Delta {\mathbf{\theta }}_c^i)}^{j-1}$ .

2. Compute joint deflections ${(\Delta {\mathbf{\theta }}_c^i)}^j$ based on equation (13).

3. If $\mathrm{norm}({(\Delta {\mathbf{\theta }}_c^i)}^j-{(\Delta {\mathbf{\theta }}_c^i)}^{j-1})\le {\omega }_0$ , iteration terminates and set $\Delta {\mathbf{\theta }}_c^i={(\Delta {\mathbf{\theta }}_c^i)}^j$ , else go to step 1, where ${\omega }_0$ is a very small predefined value.

The resulting joint equilibrium position is updated as

$${\mathbf{\theta }}_e^i={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o^i+\Delta {\mathbf{\theta }}_c^i$$

(23)

After compensating for both kinematic errors and compliance errors, the residual position errors of the robot end-effector are given as follows

$$\Delta {\mathit{X}}_{\mathrm{res}}^i={\mathit{P}}_C^i-{\mathit{P}}_M$$

(24)

where ${\mathit{P}}_C^i$ is computed based on the forward kinematics with ${\varnothing }_{\mathrm{kin}}^i$ and ${\mathbf{\theta }}_e^i$ . The convergence condition can be set as follows

$$\epsilon =\mathrm{norm}(\Delta {\mathit{X}}_{\mathrm{res}}^i-\Delta {\mathit{X}}_{\mathrm{res}}^{i-1})$$

(25)

The overall iterative process (Figure 3(b)) terminates when $\epsilon \le {\epsilon }_0$ is satisfied, where ${\epsilon }_0$ is a very small predefined value.

As mentioned before, the transition joint torques should be estimated. More precise estimation of the transition joint torques leads to better approximation of the nonlinear joint behavior, which results in a good modeling of the compliance errors. Consequently, the robot position accuracy is improved by compensating the compliance errors. Therefore, the performance of the transition torques estimation can be evaluated by the absolute position errors of the calibrated robot. Accordingly, GA is adopted to estimate the transition torques due to its advantages, such as exploring many regions of search space simultaneously, lowering sensitivity to local minima, and not requiring calculation of the derivatives. The flowchart of the basic GA is shown in Figure 3(a).

By combining GA and the least-squares algorithm, the defined kinematic parameters, joint compliance, and transition torques can be identified simultaneously. The principle of the proposed hybrid algorithm is shown in Figure 3 and can be described as follows:

1. Load data set, which in our case includes the nominal kinematic parameters, the masses of the load and links, the mass centers of links, and the measurement data consisting of the position vectors ${\mathit{P}}_M$ of the robot end-effector and the joint encoder readings ${\mathbf{\theta }}_E$ related to $m$ measurement poses.

2. Generate initial population pool; the individual is the transition joint torque vector $[\begin{array}{cc}\hfill {\tau }_2\hfill & \hfill {\tau }_3\hfill \end{array}]$ .

3. Compute the fitness value related to each individual based on the iterative least-squares algorithm, as shown in Figure 3(b). The fitness value is the absolute position error of the calibrated robot, that is, $\Vert {\mathit{P}}_C^i-{\mathit{P}}_M\Vert$ .

4. Evaluate the fitness of each individual.

5. If the stopping criterion is satisfied, then stop. The outputs are the transition joint torque vector yielding the maximum fitness value and the simultaneously identified kinematic parameters and joint compliance values.

6. Generate new population by applying genetic operators, selection, crossover, and mutation based on the fitness values from step 4.

7. Repeat steps 3–6 until the stopping criterion is satisfied.

Experimental setup and measurements

As shown in Figure 4, the experimental setup composes of a robot, a laser tracker, and a reflector. The robot manipulator is a 6-dof serial robot (Hyundai HA006 with repeatability ±0.05 mm), and a payload of 6 kg is attached on its end-effector, which is the maximum loading allowed. The laser tracker (Laser Tracker LTD800 with accuracy of ±5 μm/m) and reflector are arbitrarily located near the robot base and attached on the payload, respectively.

Figure 4.

Experimental setup.

The positions of the robot end-effector were measured over 90 robot configurations evenly distributed within the reachable workspace. Due to the limitation of the measurement volume of the laser tracker, the reachable workspace where data acquisition was carried out is a sub-workspace of the entire workspace. As mentioned before, more precise estimation of the transition joint torques leads to better approximation of the nonlinear joint behavior. Thus, the effective torques in each joint computed at the measurement configurations used for calibration should be a diverse group with a large range. The computed torques over 90 measurement configurations in joint 2 are more diverse and have a wider range than those in joint 3 (see Figure 5). Accordingly, the 90 measurement configurations were arranged based on sorting their relevant torques in joint 2 in ascending order. Every other pose within the arranged 90 measurement poses is then selected. The resultant 45 poses (denoted by ${\mathit{Q}}_1$ ) are used for robot calibration, while the other 45 poses (denoted by ${\mathit{Q}}_2$ ) are used to validate the calibration results. The torques shown in Figure 5 were computed based on the nominal kinematic parameters under the assumption that those nominal kinematic parameters are close to the real kinematic parameters.

Figure 5.

Computed torques in joints 2 and 3.

Experimental calibration results and validation

This section presents the calibration results obtained by the proposed hybrid calibration algorithm. In addition, to validate the calibration results, the positions of the end-effector predicted by the calibrated robot over the arbitrary robot configurations (not used in calibration) are compared with the measured positions. Moreover, to show the benefits of compensating the joint compliance errors, the results of both calibration and validation are compared with those obtained by the conventional kinematic calibration.

Experimental calibration results

The measurements related to the pose set ${\mathit{Q}}_1$ were used to calibrate the robot HA006. As mentioned before, the laser tracker and reflector are arbitrarily located around the robot base and on the payload, respectively. First, the base and tool alignment was carried out by calibrating the kinematic parameters in equations (4) and (5). The mean of the residual position errors after base and tool alignment is considered the position accuracy of the robot before calibration. As shown in Figure 6, the robot position accuracy over the pose set ${\mathit{Q}}_1$ before calibration is 3.994 mm, and the maximum position error is 8.936 mm.

Figure 6.

Residual position errors over pose set ${\mathit{Q}}_1$ before and after conventional kinematic calibration.

Second, in order to show the benefits of the calibration dealing with both kinematic and compliance errors, the conventional kinematic calibration was carried out to identify 27 kinematic parameters, that is, equation (7) is solved by the least-squares approach, with the pseudo-inverse that is similar to equation (19). This operation is applied iteratively until convergence. Figure 6 shows the residual position errors over the pose set ${\mathit{Q}}_1$ before and after the conventional kinematic calibration. The mean of residual position errors is reduced from 3.994 to 0.171 mm; the standard deviation is 0.085 mm, and the maximum position error is reduced from 8.936 to 0.379 mm.

Third, the proposed calibration algorithm was used to calibrate the Hyundai robot HA006. We employed Optimization Tool Box (MATLAB R2008a) as the implementation of GA. As shown in Figure 5, the range of $\left|{\tau }_2\right|$ and $\left|{\tau }_3\right|$ is $[\begin{array}{cc}\hfill 0.6\hfill & \hfill 308.014\hfill \end{array}]$ (N m) and $[\begin{array}{cc}\hfill 15\hfill & \hfill 106.115\hfill \end{array}]$ (N m), respectively. To achieve a rapid convergence, we set the bounds for the solutions as $30\le \left|{\tau }_2\right|\le 300(\mathrm{N}\mathrm{m})$ and $30\le \left|{\tau }_3\right|\le 100(\mathrm{N}\mathrm{m})$ . Consequently, the mutation function is set as the “Adaptive feasible” option. The initial range for the entries of the vectors in the initial population (population size: 50) is set as $[30;300]$ . The crossover rate is tuned experientially and set to 0.7. The stopping criteria are set to allow 50 generations and 20 stall generations. Default settings are used for the other options. We observed that GA stopped within few minutes. The simultaneously identified kinematic parameters and joint compliance parameters are shown in Tables 2 and 3, respectively. Figure 7 compares the residual position errors over the pose set ${\mathit{Q}}_1$ after applying the proposed calibration algorithm with those after applying the conventional kinematic calibration. The mean and maximum of the residual position errors are further reduced from 0.171 to 0.080 mm and from 0.379 to 0.133 mm, respectively; the standard deviation is 0.026 mm.

Table 2.

Identified D-H parameters of Hyundai robot HA006 (–: undefined; x: unselected).

i	${\alpha }_{i-1}$ (°)	$a_{i-1}$ (mm)	${\beta }_{i-1}$ (°)	$b_{i-1}$ (mm)	$d_i$ (mm)	$\Delta {\theta }_i$ (°)
S-1	0.1	0.02	0.2	0.01	0.50	0.1
2	90.2	0.322	–	–	0.001	0.2
3	0.2	0.868	0.1	–	0 (x)	0.1
4	89.9	0.203	–	–	1.032	0.2
5	−90.2	0.002	–	–	0.001	0.1
6	89.8	0.001	–	–	0.185 (x)	0 (x)
6-E	–	−0.03	–	0.02	0.05	–

Table 3.

Identified transition joint torques and joint stiffness.

	${\tau }_t$ (N m)	$k_l$ (N m/rad)	$k_h$ (N m/rad)
Joint 2	167.228	225,674	243,369
Joint 3	80.425	89,238	101,248

Figure 7.

Residual position errors over pose set ${\mathit{Q}}_1$ after calibration using the proposed algorithm.

Experimental validation

In the robot calibration experiment, the mean of the end-effector position errors over the arbitrary configurations (not used in robot calibration) indicates how accurately the calibrated robot manipulator will perform. Thus, to validate the calibration accuracy obtained by the proposed calibration algorithm, the calibrated robot was used to predict the positions of the end-effector over the pose set ${\mathit{Q}}_2$ . First, the joint deflections in joints 2 and 3 were interactively computed by using equations (11)–(13) with the identified parameters (Tables 2 and 3) and were shown in Figure 8. The positions of the end-effector predicted by the calibrated robot were then computed based on the forward kinematics with the identified kinematic parameters (Table 2) and the joint angles ${\mathbf{\theta }}_e$ ( ${\mathbf{\theta }}_e={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o+\Delta {\mathbf{\theta }}_c$ ). The mean of the residual position errors over the pose set ${\mathit{Q}}_2$ after compensating for both kinematic and compliance errors is considered the position accuracy of the calibrated robot. As shown in Figure 9, the position accuracy of the calibrated robot is 0.091 mm with standard deviation of 0.033 mm; the maximum position error is 0.153 mm. The robot position accuracy is significantly improved with a rate of 98% ((4.699 − 0.091)/4.699).

Figure 8.

Computed joint deflections over pose set $Q_2$ .

Figure 9.

Residual position errors over pose set $Q_2$ after calibration using the proposed algorithm (validation).

The proposed calibration algorithm allows one to observe the effects of kinematic errors and compliance errors on the robot position accuracy individually. The positions of the end-effector after only compensating kinematic parameters can be computed based on the forward kinematics with the identified kinematic parameters (Table 2) and the joint angles ${\mathbf{\theta }}_{\mathrm{kin}}$ ( ${\mathbf{\theta }}_{\mathrm{kin}}={\mathbf{\theta }}_E+\Delta {\mathbf{\theta }}_o$ ). Figure 10 shows the gaps between the residual position errors after compensating for both kinematic and compliance errors, and those after only compensating kinematic errors. These gaps are the compliance errors due to the joint deflections in Figure 8. In other words, the robot end-effector attached with a 6 kg payload was deflected 1.3 mm on average.

Figure 10.

Residual position errors over pose set ${\mathit{Q}}_2$ before and after compensating compliance errors (validation).

Figure 11 shows robot position accuracy obtained by the conventional kinematic calibration, and the comparison with that obtained by the proposed calibration algorithm. The position accuracy of the robot calibrated by the conventional kinematic calibration is 0.156 mm with standard deviation of 0.078 mm; the maximum position error is 0.421 mm. The effect of compliance errors on robot position accuracy is not as clear in Figure 11 as in Figure 10 because the unmodeled errors (e.g. compliance errors) are propagated to the calibration results, that is, the kinematic parameters identified by the conventional kinematic calibration absorb these errors. Position accuracy of the robot calibrated using the proposed algorithm was improved with a rate of 41.67% compared to that obtained using the conventional kinematic calibration.

Figure 11.

Comparison of position accuracy using conventional kinematic calibration and using the proposed calibration algorithm.

Conclusion

This article presents a practical technique for enhancing robot position accuracy by compensating both kinematic errors and compliance errors in industrial serial robots. In contrast to previous works, the robot joint is modeled as a piecewise linear system to approximate its nonlinear behavior. Based on this new joint stiffness model, a compliance error model is established. A hybrid least-squares-GA-based algorithm is then developed that allows one to simultaneously identify the defined geometrical and non-geometrical parameters through a comprehensive error model consisting of kinematic and compliance errors. In addition, to increase identification accuracy, measurement data used in robot calibration are selected properly. Finally, the advantages of the developed techniques are illustrated by experimental study of the Hyundai industrial 6-dof serial robot. The results show that the position accuracy of the robot after compensating for both kinematic and compliance errors is significantly improved at a rate of 98%. This represents an improvement of more than 41% in robot position accuracy obtained by the proposed calibration algorithm compared with the conventional kinematic calibration.