Kinematic calibration and feedforward control of a heavy-load manipulator using parameters optimization by an ant colony algorithm

Abstract

Most of the currently available three-degree-of-freedom manipulators are light load and cannot achieve full continuous rotation; given this, we designed a heavy-load manipulator that achieves unrestricted and continuous rotation. Due to manufacturing and assembly errors, parameter deviations between the real manipulator and its underlying theoretical model were unavoidable. Because of the lack of high-precision, high-frequency, and real-time closed-loop detection methods, we proposed a type of kinematics calibration of parameterized ant colony optimization and feedforward control methods. This was done to achieve high-precision motion control. First, an error model combining structural parameters and joint output angles was established, and the global sensitivity of each error source was analyzed to distinguish both primary and secondary sources. Based on the measured data of a laser tracker, the ant colony optimization was then used to identify six error sources. This resulted in both link length and joint driving errors of the designed manipulator. As it is a type of systematic error, the rounding error of the theoretical trajectory was carefully analyzed, and feedforward control methods with different coefficients were designed to further improve positioning accuracy based on the kinematic calibration. Experimental results showed that the proposed kinematic calibration and feedforward control methods achieved relatively precise motion control for the designed manipulator.

1. Introduction

Compared to humans, designed manipulators have many benefits including high adaptability, high accuracy, low-production cost, and low-error rate. However, in practice, manipulators have varying parameter inaccuracies compared to their underlying theoretical models [1–6]. These discrepancies are due to a host of reasons, including tolerance design, machining errors, assembly errors, and mechanical wear [7–9]. In traditional control systems, the kinematic analysis is always based on theoretical models [10]. Both geometric and nongeometric deviations between the theoretical model and actual manipulator will severely affect motion accuracy as well as system performance. Therefore, it is necessary to quantify and compensate actual errors to improve motion accuracy. The kinematic calibration and precise control of manipulators has become one of the key problems that need to be solved in the industrial application of manipulators [11, 12]. Currently, the dominant solution has been to optimally estimate the specific structural parameters of manipulators [13]. The basic process for this approach includes establishing an error model as well as data measurement, error identification, and error compensation.

According to different research objectives, both domestic and foreign scholars have adopted different calibration methods. For instance, Liu et al. established an error model that considered geometry parameters and attitude angles, which reduced the motion error of the UR10 robot from a maximum error of 5.60 mm to 0.5 mm [14]. Wang et al. proposed a relative pose-error model that integrated relative distance error and rotation error. This approach achieved a better calibration effect compared to the distance error model [15]. Chen et al. developed an error model including the positioning error of base coordinate system, the kinematic parameter error, and the end positioning error of the IRB120 robot [16]. Nguyen et al. built an error model considering the identified error sources and calibrated an industrial HH800 robot [17]. Ye et al. proposed a geometric error modeling method for super-constrained hybrid robots based on the product of exponential formulas, which reduced approximately 90% of all pose errors [18]. Xiao et al. proposed a linear error model based on geometric error parameters. This method greatly reduced the error residuals of the delta robot in two different working domains [19]. Taken together, the above studies focused on commercial robot calibration, highlighting the lack of investigation into calibration of self-designed manipulators.

Data measurement also varies given these measurements are taken using different approaches. For instance, Kamali et al. identified the exact kinematic model and joint flexibility by measuring robot position at optimal design points. This method significantly improved overall positioning accuracy [20]. Filion et al. used a commercial portable photogrammetry system along with a laser tracker to calibrate industrial robots and compared their calibration performance [21]. Separately, the accuracy of different measurement tools is an essential factor impacting the effectiveness of kinematic calibration. Therefore, high-precision and high-resolution measurement tools such as laser trackers [22, 23] and monocular cameras to the systematic error caused by these approaches to data measurement, which can weaken the overall [24, 25], are typically selected for data measurement. In fact, few scholars have paid much attention to the systematic error caused by these approaches to data measurement, which can weaken the overall precision control of robot motion.

For error identification, domestic and foreign scholars have also proposed different methods. For example, Li et al. proposed a calibration method for industrial robots based on the principle of circumferential closure. This approach used the LM algorithm to identify robot errors [26]. Sun et al. proposed a general error modeling method based on finite and instantaneous spiral theory and designed two error identification algorithms to effectively improve robot accuracy. focused on the calibration of structural parameters, including length of links and mounting angle errors [27]. In short, most previous work has focused on the calibration of structural parameters, including length of links and mounting angle errors.

After identifying the error sources, Nguyen et al. tested the effect of planar local calibration and six-dimensional pose calibration through extensive experiments on three different structured robots [28]. Wu et al. developed an accurate model of a 5-degree-of-freedom hybrid painting robot based on Lagrange’s equation and the principle of friction discrimination. They then designed a multichannel feedforward controller to compensate for the trajectory tracking errors from different sources [29]. Other scholars have proposed other solutions [30–32]. However, few studies have paid attention to the rounding error of the robot servo control process.

From the above analysis, it is clear that many studies have focused on the calibration of structural parameters, including length of links and mounting angle errors. However, the calibration of nongeometric errors has been less well documented. Based on the analysis of the designed manipulator, nongeometric errors mainly include the flexibility of the robot joints, rounding error, joint clearance, and thermal effects. In practice, nongeometric parameter errors may have a smaller impact in the early stages of application; however, their impacts become more pronounced as the operational time increases. Therefore, a comprehensive study of geometric and nongeometric parameter errors is of research and engineering value for the high-precision control of a designed manipulator.

Numerous studies have shown that there is a specific relationship between joint clearance, working time, and rotational speed. The thermal effect is affected by temperature, contact surface, and other factors. During calibration, it is difficult to find an accurate identification value for rounding error, joint clearance, and thermal effect; given this, only an average value is identified. To improve control operability, we equate here the effects arising from rounding error and joint clearance as joint driving errors.

In this paper, we designed a heavy-load manipulator with the ability of unrestricted continuous rotation. First and as shown in Section 3, an error model considering link length and joint driving errors was developed for a designed three-degree-of-freedom manipulator. Global sensitivity analysis was then performed for each error source to distinguish both primary and secondary sources. Second, we identified error sources by applying an ant colony optimization algorithm, and the calibration effect was experimentally verified as shown in Section 4. Third and based on the correlation analysis of actual trajectory errors, we designed a feedforward controller to further improve the motion accuracy as shown in Section 5. Finally, experimental results showed that these proposed methods achieved relatively precise control of the three-degree-of-freedom manipulator. This approach also provided a solution for the kinematics calibration and feedforward control of similar manipulators.

2. Model of designed manipulator

A three-dimensional view of the designed manipulator is shown in Fig. 1.

Figure 1. The three-dimensional view of designed manipulator.

The manipulator was mainly composed of a base, a first link, a second link, a third link, a first motor, a second motor, a third motor, a first reducer, a second reducer, a third reducer, and a support platform. The electric slip ring is not included in Fig. 1, which was used to realize the capability of unrestricted continuous rotation of designed manipulator. Of these, the role of the third link, motor, and reducer was to ensure that the laser tracker could follow the full continuous rotation of the designed manipulator without dead-end tracking.

3. Error model of designed manipulator

3.1. Establishment of the error model

The calibration error is actually the differences between the theoretical length and the mean projected length in the x-y plane of the manipulator. Errors caused by bending due to overload, machining, and assembly are considered and reflected in the projected length, and the vertical deformation can be compensated with an additional fourth joint for spatial vertical motion. By neglecting the effect of the structural bending of the manipulator links and joints due to the heavy load, the error model is established for the designed manipulator.

As shown in Fig. 2(a), let the theoretical coordinate vector of the manipulator endpoint be O ₁ A and the actual coordinate vector be O ₁ A′. Due to machining and assembly errors, the actual parameters of the manipulator deviated from the underlying theoretical parameters. Let the theoretical length and error of the first link be $l_{1}$ and $\delta l_{1}$ , respectively, the theoretical length and error of the second link be $l_{2}$ and $\delta l_{2}$ , respectively, and the theoretical length and error of the third link be $l_{3}$ and $\delta l_{3}$ , respectively. It was easy to know that in the x-y plane projection, there existed $l_{3}=0$ .

Figure 2. The error sources of designed manipulator. (a) Schematic diagram of plane projection of six error sources. (b) Three-dimensional schematic diagram of six error sources.

Let the theoretical angles and errors of the first link with respect to the y-axis be $\theta _{1}$ and $\delta \theta _{1}$ , respectively, the theoretical angles and errors of the second link with respect to the y-axis be $\theta _{2}$ and $\delta \theta _{2}$ , respectively, and the theoretical angles and errors of the third link with respect to the y-axis be $\theta _{3}$ and $\delta \theta _{3}$ , respectively. The three-dimensional schematic diagram of the six error sources of designed manipulator is shown in Fig. 2(b).

There are output rounding errors of the servo motor and gear backlash of the planetary reducer. These will lead to joint driving errors $(\delta \theta _{1},\delta \theta _{2},\delta \theta _{3})$ of the three manipulator joints [9, 13]. Due to the presence of both machining and assembly errors, the link length of the manipulator had varying parameter inaccuracies $(\delta l_{1},\delta l_{2},\delta l_{3})$ compared to the theoretical models. The generation process of both joint driving and link length errors is shown in Fig. 3.

Figure 3. The generation process of joint driving errors and link length errors.

The section diagram of the designed first link is shown in Fig. 4, which demonstrates the assembly details of the first link. One can see that the locking assembly connecting planetary reducer and flange can transmit the torque by preloading the bolt. When the bolt preload is weakened, it is easy to generate joint clearance. In addition, although interference fit is carefully designed, the bearing connecting the flange and bearing block can also generate joint clearance. This is due to machining errors and mechanical wear. To improve the control operability, this paper equates the effects arising from joint clearance and rounding errors as joint driving errors. The analysis of rounding error will be shown in Section 5.1.

Figure 4. The section diagram of designed first link.

Under a desirable condition, the end position of the manipulator can be expressed as

(1) \begin{equation} \left\{\begin{array}{l} \left(O_{1}A\right)_{x}=l_{1}\sin \theta _{1}+l_{2}\sin \theta _{2}\\[5pt] \left(O_{1}A\right)_{y}=l_{1}\cos \theta _{1}+l_{2}\cos \theta _{2} \end{array}\right. \end{equation}

When the six error sources set above exist in the manipulator, it can be obtained from equation (1) as

(2) \begin{equation} \left\{\begin{array}{l} \left(O_{1}A'\right)_{x}=\left(l_{1}+\delta l_{1}\right)\sin \left(\theta _{1}+\delta \theta _{1}\right)+\left(l_{2}+\delta l_{2}\right)\sin \left(\theta _{2}+\delta \theta _{2}\right)+\delta l_{3}\sin \left(\theta _{3}+\delta \theta _{3}\right)\\[5pt] \left(O_{1}A'\right)_{y}=\left(l_{1}+\delta l_{1}\right)\cos \left(\theta _{1}+\delta \theta _{1}\right)+\left(l_{2}+\delta l_{2}\right)\cos \left(\theta _{2}+\delta \theta _{2}\right)+\delta l_{3}\cos \left(\theta _{3}+\delta \theta _{3}\right) \end{array}\right. \end{equation}

Then, the global error vector can be expressed as

(3) \begin{equation} \left\{\begin{array}{l} \delta x=\left(l_{1}+\delta l_{1}\right)\sin \left(\theta _{1}+\delta \theta _{1}\right)+\left(l_{2}+\delta l_{2}\right)\sin \left(\theta _{2}+\delta \theta _{2}\right)\\[5pt]\qquad +\,\delta l_{3}\sin \left(\theta _{3}+\delta \theta _{3}\right)-l_{1}\sin \theta _{1}-l_{2}\sin \theta _{2}\\[5pt] \delta y=\left(l_{1}+\delta l_{1}\right)\cos \left(\theta _{1}+\delta \theta _{1}\right)+\left(l_{2}+\delta l_{2}\right)\cos \left(\theta _{2}+\delta \theta _{2}\right)\\[5pt]\qquad +\,\delta l_{3}\cos \left(\theta _{3}+\delta \theta _{3}\right)-l_{1}\cos \theta _{1}-l_{2}\cos \theta _{2} \end{array}\right. \end{equation}

After simplification and reorganization by the vector representations, an error model proposed in this paper can be obtained as

(4) \begin{equation} \left\{\begin{array}{l} \delta x=\boldsymbol{\delta }\textbf{L}\left(\textbf{P}_{\textbf{1}}\left(\boldsymbol{\delta \theta }\right)\textbf{-}\textbf{P}_{\textbf{2}}\right)\textbf{E}_{{\boldsymbol{\theta }_{1}}}\\[5pt] \delta y=\boldsymbol{\delta }\textbf{L}\left(\textbf{P}_{\textbf{3}}\left(\boldsymbol{\delta \theta }\right)\textbf{-}\textbf{P}_{2}\right)\textbf{E}_{{\boldsymbol{\theta }_{2}}} \end{array}\right. \end{equation}

Where $\boldsymbol{\delta }\textbf{L}=[l_{1},l_{2},l_{3},\delta l_{1},\delta l_{2},\delta l_{3}], \boldsymbol{\delta \theta }=\left[\begin{array}{c@{\quad}c} \delta \textbf{E}_{\textbf{c}\boldsymbol{\theta }} & \textbf{0}_{3\times 3}\\[5pt] \textbf{0}_{3\times 3} & \delta \textbf{E}_{\textbf{s}\boldsymbol{\theta }} \end{array}\right]$ ,

$\textbf{P}_{1}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1} & \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1}\\[5pt] \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1} & \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{E}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{E}_{1\times 1} \end{array}\right], \textbf{P}_{2}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1} & \textbf{0}_{2\times 2} & \textbf{0}_{2\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} \end{array}\right]$ ,

$ \textbf{P}_{3}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1} & -\textbf{E}_{2\times 2} & \textbf{0}_{2\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1} & \textbf{0}_{1\times 2} & \textbf{0}_{1\times 1}\\[5pt] \textbf{E}_{2\times 2} & \textbf{0}_{2\times 1} & -\textbf{E}_{2\times 2} & \textbf{0}_{2\times 1}\\[5pt] \textbf{0}_{1\times 2} & \textbf{E}_{1\times 1} & \textbf{0}_{1\times 2} & -\textbf{E}_{1\times 1} \end{array}\right]$ .

$\textbf{E}_{{\boldsymbol{\theta }_{1}}}=[\sin \theta _{1},\sin \theta _{2},\sin \theta _{3},\cos \theta _{1},\cos \theta _{2},\cos \theta _{3}]^{T}$ ,

$\textbf{E}_{{\boldsymbol{\theta }_{2}}}=[\cos \theta _{1},\cos \theta _{2},\cos \theta _{3},\sin \theta _{1},\sin \theta _{2},\sin \theta _{3}]^{T}$ . Where $\textbf{E}_{i\times i}$ denotes the unit matrix of order i, $\delta \textbf{E}_{s\theta }= Diag (\sin \delta \theta _{1},\sin \delta \theta _{2},\sin \delta \theta _{3}), \delta \textbf{E}_{c\theta }= Diag (\cos \delta \theta _{1},\cos \delta \theta _{2},\cos \delta \theta _{3})$ .

In conclusion, the effects of the six error sources were considered on the absolute positioning accuracy of the manipulator. These were the length errors of the first, second, and third links and the joint driving errors of the first, second, and third joints. The transfer function between each error source and the absolute positioning error is shown in equation (4).

3.2. Sensitivity analysis of error source

As shown in equation (4), both the length and joint driving errors had an important influence on the absolute positioning accuracy of the manipulator. However, the degree of influence of these six error sources on the positioning accuracy could not be quantified from equation (4) alone. Given this, it became necessary to conduct sensitivity analysis on six error sources to observe the degree of their influence on the positioning accuracy. It is assumed that a system $S(\delta )$ with six error sources $\delta l_{1},\delta l_{2},\delta l_{3},\delta \theta _{1},\delta \theta _{2},\delta \theta _{3}$ is existed corresponding to equation (4). If the $\delta _{i}$ is set to be 0, the system $S(\delta )$ will degenerate to a system $S(\delta _{0}),(\delta _{i}\notin \delta _{0})$ . Then, the system $S(\delta )$ can be seen as the disturbed system from $S(\delta _{0})$ . The following equation (5) can be derived by the perturbation method from equation (4).

As shown in equation (4), the absolute positioning error transmitted to the end by each error source can be expressed as follows:

(5) \begin{equation} \left\{\begin{array}{l} \delta x\left(l_{1}\right)=\delta l_{1}\sin \left(\theta _{1}+\delta \theta _{1}\right)\\[5pt] \delta y\left(l_{1}\right)=\delta l_{1}\cos \left(\theta _{1}+\delta \theta _{1}\right)\\[5pt] \delta x\left(l_{2}\right)=\delta l_{2}\sin \left(\theta _{2}+\delta \theta _{2}\right)\\[5pt] \delta y\left(l_{2}\right)=\delta l_{2}\cos \left(\theta _{2}+\delta \theta _{2}\right)\\[5pt] \delta x\left(l_{3}\right)=\delta l_{3}\sin \left(\theta _{3}+\delta \theta _{3}\right)\\[5pt] \delta y\left(l_{3}\right)=\delta l_{3}\cos \left(\theta _{3}+\delta \theta _{3}\right)\\[5pt] \delta x\left(\theta _{1}\right)=\left(l_{1}+\delta l_{1}\right)\left(\sin \left(\theta _{1}+\delta \theta _{1}\right)-\sin \theta _{1}\right)\\[5pt] \delta y\left(\theta _{1}\right)=\left(l_{1}+\delta l_{1}\right)\left(\cos \left(\theta _{1}+\delta \theta _{1}\right)-\cos \theta _{1}\right)\\[5pt] \delta x\left(\theta _{2}\right)=\left(l_{2}+\delta l_{2}\right)\left(\sin \left(\theta _{2}+\delta \theta _{2}\right)-\sin \theta _{2}\right)\\[5pt] \delta y\left(\theta _{2}\right)=\left(l_{2}+\delta l_{2}\right)\left(\cos \left(\theta _{2}+\delta \theta _{2}\right)-\cos \theta _{2}\right)\\[5pt] \delta x\left(\theta _{3}\right)=\delta l_{3}\left(\sin \left(\theta _{3}+\delta \theta _{3}\right)-\sin \theta _{3}\right)\\[5pt] \delta y\left(\theta _{3}\right)=\delta l_{3}\left(\cos \left(\theta _{3}+\delta \theta _{3}\right)-\cos \theta _{3}\right) \end{array}\right. \end{equation}

Define the functional relationship as

(6) \begin{equation} \boldsymbol{\Delta }\textbf{=}\textbf{M}\left(\boldsymbol{\Delta }\textbf{L}\textbf{,}\boldsymbol{\Delta }\boldsymbol{\Theta }\textbf{,L,}\boldsymbol{\Theta }\right) \end{equation}

where $\boldsymbol{\Delta }=[\delta x(L,\theta ),\delta y(L,\theta )], \delta x(L,\theta )=[\delta x(l_{1}),\delta x(l_{2}),\delta x(l_{3}),\delta x(\theta _{1}),\delta x(\theta _{2}),\delta x(\theta _{3})]$ , $\delta y(L,\theta )=[\delta y(l_{1}),$ $\delta y(l_{2}),\delta y(l_{3}),\delta y(\theta _{1}),\delta y(\theta _{2}),\delta y(\theta _{3})]$ , $\textbf{L}=(l_{1},l_{2},l_{3})$ , $\boldsymbol{\Theta }=(\theta _{1},\theta _{2},\theta _{3})$ , $\boldsymbol{\Delta }\textbf{L}=(\delta l_{1},\delta l_{2},\delta l_{3})$ , $\boldsymbol{\Delta \Theta }=(\delta \theta _{1}$ , $\delta \theta _{2},\delta \theta _{3})$ .

The coefficient of local sensitivity can be defined as $k(\delta, x, y)$ , for example:

(7) \begin{equation} k\left(\delta x\left(l_{1}\right)\!,x,y\right)=\left| \delta x\left(l_{1}\right)/\delta x\right| \end{equation}

Other coefficients of sensitivity are similar, where $(\theta _{1},\theta _{2})=f_{1}^{-1}(x,y), \theta _{3}=-(\theta _{1}+\theta _{2})+\theta _{\theta }$ , and $\theta _{\theta }$ is the angle of the connecting line between the laser tracker and the end of third link with respect to the y-axis. Here, $f_{1}$ is defined as inverse kinematics equation.

As indicated, the coefficient of local sensitivity was closely related to the position coordinates (x,y). Consequently, the coefficient of local sensitivity for the same error source was different at various points. To fully reflect the degree of influence on the positioning error in the defined workspace, the coefficient of global sensitivity can be defined as [33, 34]:

(8) \begin{equation} K\left(\delta \right)=\left[\int _{y_{\text{lower}}}^{y_{\text{upper}}}\int _{x_{\text{lower}}}^{x_{\text{upper}}}k\left(\delta,x,y\right)dxdy\right]\bigg/\left[\int _{y_{\text{lower}}}^{y_{\text{upper}}}\int _{x_{\text{lower}}}^{x_{\text{upper}}}dxdy\right] \end{equation}

where $K(\delta )$ is the coefficient of global sensitivity corresponding to the error term $\delta$ . After numerical calculation, the global degree of influence of each error source on the positioning error was then obtained.

3.3. Analysis of global error sensitivity

We set the lengths of the first, second, and third links to 0.3 m, 0.3 m, and 0 m, respectively. The working space was set as follows: $x_{\text{lower}}=y_{\text{lower}}=-0.4\,\textrm{m}, x_{\text{upper}}=y_{\text{upper}}=0.4\,\textrm{m}$ . The link length errors were set to 0.001 m, and the joint driving errors were set to 0.001 rad. Then, the histograms of coefficient of global error sensitivity were plotted in Fig. 5 using equation (7). As control groups and as shown in Fig. 6, the link length error was set to 0.001 m and the joint driving error was set to 0.01 rad. As shown in Fig. 7, the link length error was set to 0.0001 m and the joint driving error was set to 0.001 rad.

Figure 5. Histogram of global sensitivity coefficient. (Link length error of 0.001 m, joint driving error of 0.001 rad). (a) Sensitivity coefficient of each error source to $\delta$ x. (b) Sensitivity coefficient of each error source to $\delta$ y.

Figure 6. Histogram of global sensitivity coefficient (Link length error of 0.001 m, joint driving error of 0.01 rad). (a) Sensitivity coefficient of each error source to $\delta x$ . (b) Sensitivity coefficient of each error source to $\delta y$ .

As shown in Figs. 5–7, the error sources that had the greatest impact on the absolute positioning accuracy were mainly from the first and second link length errors $(\delta l_{1},\delta l_{2})$ as well as the first and second joint driving errors $(\delta \theta _{1},\delta \theta _{2})$ . When the link length errors were 0.001 m and the angle errors were 0.001 rad (Fig. 5), the link length errors $(\delta l_{1},\delta l_{2})$ were the main influencing factors. When the length errors were 0.001 m and the angle errors were 0.01 rad (Fig. 6), the joint driving errors $(\delta \theta _{1},\delta \theta _{2})$ were the main influencing factors. When the length errors were 0.0001 m and the angle errors were 0.001 rad (Fig. 7), the angle errors $(\delta \theta _{1},\delta \theta _{2})$ were the main influencing factors. It is worth noting that according to the following identification experiment of error sources, the range of link length error was from 0.0005 mm to 0.002 mm. Moreover, the range of joint driving error was from 0.005 rad to 0.02 rad. Given this, the above six error sources needed to be carefully identified.

To clarify, the coefficients $\mu (x)$ and $\mu (y)$ are defined to describe the ratio of the influence degree of

the link length error to the joint driving error. This is calculated as follows equation (9).

(9) \begin{equation} \left\{\begin{array}{l} \mu \left(x\right)=\dfrac{\sum _{i=1}^{3}K\left(\delta x\left(\theta _{i}\right)\right)}{\sum _{i=1}^{3}K\left(\delta x\left(l_{i}\right)\right)}\\[15pt] \mu \left(y\right)=\dfrac{\sum _{i=1}^{3}K\left(\delta y\left(\theta _{i}\right)\right)}{\sum _{i=1}^{3}K\left(\delta y\left(l_{i}\right)\right)} \end{array}\right. \end{equation}

The link length error was selected from 0.0001 m to 0.002 m with 500 discrete points; similarly, the joint driving error was selected from 0.001 rad to 0.02 rad with 500 discrete points. The coefficients $\mu (x)$ and $\mu (y)$ varied with both the link length and joint driving errors (Fig. 8), where $\mu = 1$ is the reference plane.

Figure 7. Histogram of global sensitivity coefficient (Link length error of 0.0001 m, joint driving error of 0.01 rad). (a) Sensitivity coefficient of each error source to $\delta x$ . (b) Sensitivity coefficient of each error source to $\delta y$ .

Figure 8. Ratio of sensitivity coefficient of link length error to joint driving error. (a) In x-direction. (b) In y-direction.

As shown in Fig. 8(a), the coefficient $\mu (x)$ increased with increasing joint driving error and decreasing link length error in the x-direction. In the y-direction and as shown in Fig. 8(b), $\mu (y)$ also generally conformed to the trend observed in Fig. 8(a). However, there were some peaks and valleys where coefficients $\mu (x)$ and $\mu (y)$ changed abruptly in Fig. 8(b). This showed some randomness and did not necessarily follow the trend shown in Fig. 8(a). Taken together, this shows that the influence of the above six error sources needs to be fully considered when performing kinematic calibration.

4. Kinematic calibration based on the error model

The applications of kinematic calibration aim to improve the positioning accuracies of the manipulator. Global calibration and local calibration are introduced and the experiments are conducted, respectively. All points in the calibration process are measured continuously under the state of motion for the manipulator, where the laser tracker is used with the sampling frequency of 100 Hz. Before the experiments, the path is preset according to the expected trajectory, and a waiting time of 1 s is set at each point. The laser tracker collects the coordinates of the end continuously on the trajectory in the process of experiment. For the stage of data processing, the measured coordinate of each discrete node is used as the mean value, which is calculated by the 100 collected coordinates of this point within an interval of 0.5 s. Furthermore, the feedforward control will be added to reduce errors in the calibration. The influence of the controller on the manipulator is still reflected in the kinematic calibration, although it is not explicitly mentioned and analyzed as discrete error or rounding error. In fact, the six error sources in the kinematic calibration represent the calibrated mean values considering all influencing factors, including the influence of the controller, driver, backlash of the reducer, and other factors.

4.1. Global calibration

In the global calibration, the selection of calibration points aims to cover the entire workspace, where the center of the workspace is chosen as the center of the circle. It is the most effective way to avoid interference from the other factors. Choosing other positions as the center of the workspace will result in the loss of universality for global calibration and have certain impact on its positioning accuracy. Furthermore, the selection of calibration points should be closely determined by the workspace to ensure the universality of global calibration. A chessboard pattern is commonly used for a square workspace, where an equilateral triangle pattern is usually employed for an equilateral triangle workspace. Similarly, the selection of calibration points would form concentric circles for a circular workspace. The calibration path is set in a way between two adjacent points for the process of calibrations, and the manipulator returns to the origin point to eliminate accumulated errors for each calibration point.

Figure 9. Schematic diagram of the discrete points of global calibration.

According to equations (2) and (3), a certain number of discrete points were selected to identify the above six error sources. Considering the whole working domain of three-degree-of-freedom manipulator and as shown in Fig. 9, the discrete points of global calibration should be carefully selected. Specifically, the circular radius was set with equal distance $\Delta r$ around the circle center O, and the discrete points were set with equal circular angle $\Delta \theta$ . Considering eliminating cumulative errors, the calibration path is set as $O \rightarrow P_{1,1} \rightarrow O \rightarrow P_{1,2} \rightarrow \cdots O \rightarrow P_{i,j} \cdots O \rightarrow P_{I,J} $ , where $P_{i,j}$ is the discrete point. The polar diameter $i\Delta r$ and polar angle $j\Delta \theta$ can be expressed in rectangular coordinates as

(10) \begin{equation} \left\{\begin{array}{l} x_{i,j}=\mathit{i}\Delta \mathit{r}\cos j\Delta \theta \\[5pt] y_{i,j}=\mathit{i}\Delta \mathit{r}\sin j\Delta \theta \end{array}\right. \end{equation}

According to the parameters of the designed manipulator, the theoretical working domain was a circle with a radius of 0.6 m. Then, the equal circular angle $\Delta\theta$ was set to 20°, and the equal distance $\Delta r$ was set to 0.05 m. This resulted in 216 discrete points. The position coordinates of the 216 discrete points for global calibration were obtained from equation (10). Adding the position coordinates of every origin point, we obtained the coordinate data for 432 discrete points. As shown in Fig. 10, the theoretical position coordinates were sent to the control program and were processed by inverse kinematics to control the manipulator. Then, a laser tracker was used to measure the actual coordinate data of the 432 discrete points.

Figure 10. Schematic diagram of error identification.

Furthermore, we considered that the six error sources were the identification objective, and the optimization objective was determined as the sum of positioning errors between the fitted and calibration points. As an optimization algorithm with the advantages of distributed computing and high operational efficiency, the ant colony algorithm has been widely used in the field of robot control [35–37]. This paper adopted the ant colony optimization algorithm for the identification of six error sources. The specific data processing and error identification processes are shown in Fig. 10.

In the specific process of ant colony optimization, the objective function is set as

(11) \begin{equation} \text{fit}=\sum _{i=1}^{N}\sqrt{\left(\delta x_{i}\right)^{2}+\left(\delta y_{i}\right)^{2}} \end{equation}

The equation calculating the probability of state transition is demonstrated as

(12) \begin{equation} p_{j}=\frac{\tau _{j}\left(t\right)-\min \left(\tau _{j}\left(t\right)\right)}{\tau _{j}\left(t\right)} \end{equation}

When the probability of state transfer $p_j$ of ant j is less than the constant $p_{0}$ , the ant j belongs to the better one in ith iteration, and it performs local search. When the probability of state transfer $p_{j}$ is larger than the constant $p_{0}$ , the ant j belongs to the worse one in ith iteration, and it performs global search. The pheromone $\tau _{j}$ can be updated as

(13) \begin{equation} \tau _{j}\left(t+1\right)=\left(1-\rho \right)\tau _{j}\left(t\right)+\text{fit}\left(\delta x,\delta y,l_{1},l_{2},\theta _{1},\theta _{2}\right) \end{equation}

where $\rho$ is the pheromone volatility factor. During each iteration, the pheromone $\tau _{j}$ appears to be increased by the value of objective function fit. However, to avoid premature maturation and falling into local optima between generations, the pheromone $\tau _{j}$ will not be fully inherited. The ant colony optimization is taken from the previous works [35–37], and its specific progress is shown in Fig. 11, which outputs the optimal solution $\left(k_{\text{1best}},k_{\text{2best}},k_{\text{3best}},k_{\text{4best}},k_{\text{5best}},k_{\text{6best}}\right)$ as six error sources $(\delta l_{1},\delta l_{2},\delta l_{3},\delta l_{4},\delta l_{5},\delta l_{6})$ .

The ant colony size is set as $J = 10,000$ , with the maximum number of iterations $T = 1000$ , the pheromone volatility factor $\rho = 0.9$ , and the transfer probability constant $p 0 = 0.15$ . All adjustable parameters mentioned in the paper were selected as optimal parameters based on multiple simulation comparisons of their effects.

Table I. Identification results of global calibration.

Figure 11. Flow chart of ant colony algorithm.

The initial values of six error source $(\delta l_{1},\delta l_{2},\delta l_{3},\delta l_{4},\delta l_{5},\delta l_{6})$ are set to (0, 0, 0, 0, 0, 0). After 1000 iterations for optimization, the final identification results were obtained and are shown in Table I.

The identification process of ant colony optimization is shown in Fig. 12. The optimal value of fitness function is 64.4996 mm as shown in Fig. 12(a), where the value of 216 discrete calibration points was 60.0567 mm and the value of 216 origin points was 4.4428 mm. Overall, the average positioning error of each point after identification was 0.1641 mm, where the average positioning errors of the calibration point and origin point were 0.3049 mm and 0.0226 mm, respectively.

Figure 12. The identification process of ant colony optimization. (a) fit. (b) $\delta$ l ₁ and $\delta$ l ₂. (c) $\delta$ l ₃. (d) $\delta$ θ ₁ and $\delta$ θ ₂. (e) $\delta$ θ ₃.

Compared with the initial fitness function value of 425.1475 mm, the global calibration greatly improved positioning accuracy. From Fig. 12(b)–(e), it is clear that after approximately 100 iterations for optimization, these six error sources $(\delta l_{1},\delta l_{2},\delta l_{3},\delta \theta _{1},\delta \theta _{2},\delta \theta _{3})$ converged to constant values. These identification results showed that although the identified model was more consistent with the actual model than the theoretical model, the existence of systematic errors weakened the positioning accuracy. In view of the systematic errors caused by the convergence limit of the objective function of ant colony optimization (Fig. 12(a)), there should be some specific methods used to improve the positioning accuracy in a specific workspace.

4.2. Global calibration experiment

The experimental calibration setup is shown in Fig. 13. Briefly, the laser tracking ball was mounted on top of the third link, and the laser tracker quickly measured the continuous position coordinates of the manipulator endpoint. The electric slip ring was used to achieve continuous rotation, which means that all three joint motors are capable of completing a full circular motion. This continuous rotation capability allows joint 3 adjustment to ensure that the tracking ball remains aligned with the laser tracker.

Figure 13. The experimental setup.

First, the identification results of the ant colony optimization were substituted into the inverse kinematic solution of the designed manipulator. Then, the experimental trajectories were selected as the following four groups: (1) Circle center of (0.3 m, 0 m), radius of 0.2 m; (2) circle center of (0 m, 0.3 m), radius of 0.2 m; (3) circle center of (−0.3 m, 0 m), radius of 0.2 m; (4) circle center of (0 m, −0.3 m), radius of 0.2 m. After every experimental trajectory was sent into the modified control system, the joint input data obtained by the inverse kinematic solution was processed by the control card and servo motor. Finally, the actual trajectory data were collected by a laser tracker, and the global calibration error was compared with the original error by data processing (Fig. 14). The average error and standard deviation of the above four experimental trajectory groups are shown in Table II.

As shown in Fig. 14 and Table II, the global calibration achieved superior positioning accuracy at most points. The average error and the standard deviation were substantially reduced, indicating that a large improvement in accuracy was achieved by global calibration. However, the effect of global calibration varied greatly in different trajectories. For example, the reduction percentages of average errors of the four experimental trajectories were 32.5%, 16.3%, 76.7%, and 51.0%, respectively (Table II). Moreover, the average reduction percentage was 44.1%.

Table II. Average error and standard deviation of four groups of experimental trajectories.

Figure 14. Comparison of trajectory errors between original and global calibration. (a) Circle center of (0.3 m, 0 m), radius of 0.2 m. (b) Circle center of (0 m, 0.3 m), radius of 0.2 m. (c) Circle center of (−0.3 m, 0 m), radius of 0.2 m. (d) Circle center of (0 m, −0.3 m), radius of 0.2 m.

Due to the presence of systematic errors caused by the convergence limit of ant colony optimization and measurement errors, a few points have a worse positioning accuracy. This is shown in Fig. 14(b) and (d). To further improve positioning accuracy, a method of local calibration was next considered and adopted on the basis of our global calibration.

4.3. Local calibration

Local calibration is determined by the particular trajectories with the better positioning accuracy. The selection of calibration points for local calibration is related to the desired trajectory. However, if the desired trajectory is changed, the effectiveness of local calibration will be greatly restricted correspondingly. It is not necessary to preserve the complete circular sector shape for local calibrations. In practice, the discretized points of the expected trajectory can be used as calibration points for local calibration. The calibration path should include the equal number of origin points to eliminate cumulative errors, which means that the manipulator should returns to the origin point repeatedly in the process of calibrations.

The points on the edge of the workspace have one set of solutions for planar manipulators, where all other positions have two sets of solutions. The following rules are used in the experiment to determine the choice of “pose”. (1) Choose the “pose” with the smallest sum of the absolute increments of joint motor output angles (comparing with the previous point, which is considered as the origin by default for the starting point of the trajectory). (2)When multiple “poses” have the same sum of absolute increments of joint motor output angles, the “pose” with the smallest absolute increment of the first motor’s output angle is chosen. If the absolute increment of the first motor’s output angle is still equal, that of the second motor’s output angle is compared and the smallest one is chosen, and so on.

Given the deteriorating accuracy of the global calibration across a few points, we considered selecting calibration points from specific domains rather than the whole working domain (Fig. 15). The local calibration process was similar to that shown in Fig. 10, and the optimization method was also the ant colony optimization as shown in Fig. 11. The initial values of the six error sources $(\delta l_{1},\delta l_{2},\delta l_{3},\delta \theta _{1},\delta \theta _{2},\delta \theta _{3})$ were set to (0,0,0,0,0,0).

Table III. Identification results of local calibration.

Figure 15. Schematic diagram of the discrete points of local calibration.

Figure 16. The identification process of ant colony optimization. (a) fit. (b) $\delta$ l₁ and $\delta$ l₂. (c) $\delta$ l₃. (d) $\delta$ θ₁ and $\delta$ θ₂. (e) $\delta$ θ₃.

We conducted 1000 iterations for optimization and the identification results are shown in Table III. The identification process of ant colony optimization is shown in Fig. 16. The optimal value of fitness function was 26.1605 mm (Fig. 16(a)), and the average positioning error of each point after identification was 0.2565 mm. This was better than that achieved for global calibration (Fig. 12). Compared with the initial fitness function value of 127.9186 mm, the local calibration greatly improved the positioning accuracy. From Fig. 16(b)–(e) and after about 600 iterations of optimization, these six error sources $(\delta l_{1},\delta l_{2},\delta l_{3},\delta \theta _{1},\delta \theta _{2},\delta \theta _{3})$ converged to constant values. As discussed in Section 4.1, the systematic errors caused by the convergence limit of the objective function of the ant colony optimization resulted in a residual error of 26.1605 mm (Fig. 16(a)).

4.4. Local calibration experiment

The identification parameters of local calibration were used to modify the inverse kinematic solution of the designed manipulator. To test the effectiveness of the local calibration method, a circular trajectory with circle center of (0 m, 0.3 m) and radius of 0.3 m were selected. Separately, the identification parameters for both global and local calibration were used to control the manipulator. As shown in Fig. 13, the actual trajectory data were collected by a laser tracker. Comparison of trajectory errors of different calibration methods is shown in Fig. 17.

As shown in Fig. 17, both the global and local calibration methods reduced positioning errors at most points. The positioning accuracy of local calibration was also better than that of global calibration. As shown in Table IV and compared with the original errors, the reduction percentages of average errors between global and local calibration were 45.6% and 65.8%, respectively. However and different from global calibration, local calibration only improved positioning accuracy at a specific domain. This limits its applicability across the whole workspace.

Table IV. Average error and standard deviation of different calibration methods.

Figure 17. Comparison of trajectory errors of different calibration methods.

Further analysis of Fig. 17 showed that after the global and local calibration methods, the manipulator still had large positioning errors in the red box. One of the main reasons for this was that there were systematic errors caused by the convergence limit of the objective function of ant colony optimization. The other is that there were random errors introduced by data measurement. Therefore, to further improve the positioning accuracy, a feedforward control method should be analyzed and adopted on the basis of calibration.

5. Feedforward control method

5.1. Analysis of rounding error

Before the introduction of feedforward control, it is necessary to analyze the rounding error, which is a typical systematic error. First, the theoretical trajectory is set to $(x(t),y(t))$ , where x and y are both functions of time t. Let $n$ represents the number of discrete points of a continuous trajectory, which does not include the repeated paths. The starting point position is $(x(t_{q}),y(t_{q}))$ , with $t_{q}=0$ . The ending point position is $(x(t_{z}),y(t_{z}))$ , with $t_{z}=T$ . The impact of the start and stop positions of a circular trajectory will vary in different cases. For a circular trajectory with a center at (0.3, 0) and a radius of 0.3 m, the global calibration method can apply to any trajectory in the entire workspace to control the manipulator. In this situation, it has no impact on the experimental results. On the other hand, it will have an impact on the experimental results for the local calibration or feedforward control methods.

Then and according to the theoretical manipulator model, $(x(t),y(t))$ can be obtained as

(14) \begin{equation} \left\{\begin{array}{l} x\left(t\right)=l_{1}\sin \theta _{1}\left(t\right)+l_{2}\sin \theta _{2}\left(t\right)\\[5pt] y\left(t\right)=l_{1}\cos \theta _{1}\left(t\right)+l_{2}\cos \theta _{2}\left(t\right) \end{array}\right. \end{equation}

Defining the function $(x,y)=f(\theta _{1},\theta _{2})$ , where $(x,y)$ and $(\theta _{1},\theta _{2})$ are corresponding to each other. If it can be guaranteed that $\forall t\in [0,T]$ , the output angles satisfy $(\theta _{1}(t),\theta _{2}(t))=f^{-1}(x(t),y(t))$ , then the actual position $(x_{s}(t),y_{s}(t))$ will coincide with the theoretical position $(x(t),y(t))$ . In practice, it is not possible to send a continuous signal to drive the servo motor; rather, it is a number of discrete pulses with varying time intervals. Therefore, the introduction of rounding error for discrete trajectory is inevitable.

Theoretically, when the time interval approaches an infinitely small value, the trajectory error will also tend to be infinitely small. Assuming that the time interval of pulses sent to the servo motor is $\Delta T$ , it is observed that at the node of $t=n\Delta T$ , the actual output angles $(\theta _{1s}(t),\theta _{2s}(t))$ coincide with the theoretical output angles $(\theta _{1}(t),\theta _{2}(t))$ (Fig. 18).

Figure 18. The theoretical and actual output angles.

Between $n\Delta T$ and $(n+1)\Delta T$ , servo motors were under synchronous control mode with the time of motor start and stop being small enough. This drove the first and second joints to follow the specified angles in equal time. Between $n\Delta T$ and $(n+1)\Delta T$ , the actual output angles $(\theta _{1s}(t),\theta _{2s}(t))$ and theoretical output angles $(\theta _{1}(t),\theta _{2}(t))$ are shown in Fig. 18. After solving the equation of line segments in Fig. 18, the actual output angles $(\theta _{1s}(t),\theta _{2s}(t))$ and actual position $(x_{s}(t),y_{s}(t))$ can be obtained as follows:

(15) \begin{equation} \left\{\begin{array}{l} \theta _{1s}\left(t\right)=\theta _{1}\left(n\Delta t\right)+\dfrac{t-n\Delta t}{\Delta t}\left(\theta _{1}\left(\left(n+1\right)\Delta t\right)-\theta _{1}\left(n\Delta t\right)\right)\\[9pt] \theta _{2s}\left(t\right)=\theta _{2}\left(n\Delta t\right)+\dfrac{t-n\Delta t}{\Delta t}\left(\theta _{2}\left(\left(n+1\right)\Delta t\right)-\theta _{2}\left(n\Delta t\right)\right)\\[9pt] x_{s}\left(t\right)=l_{1}\sin \left(\theta _{1s}\left(t\right)\right)+l_{2}\sin \left(\theta _{2s}\left(t\right)\right)\\[5pt] y_{s}\left(t\right)=l_{1}\cos \left(\theta _{1s}\left(t\right)\right)+l_{2}\cos \left(\theta _{2s}\left(t\right)\right) \end{array}\right.\,\,\,t\in \left[n\Delta t,\left(n+1\right)\Delta t\right] \end{equation}

Between $n\Delta T$ and $(n+1)\Delta T$ , the trajectory error $e(t)$ can be expressed as

(16) \begin{equation} e\left(t\right)=\left\| f\left(\theta _{1s}\left(t\right),\theta _{2s}\left(t\right)\right)-f\left(\theta _{1}\left(t\right),\theta _{2}\left(t\right)\right)\right\| \end{equation}

When $\Delta T$ is small enough, and $\theta (n\Delta T)$ and $\theta ((n+1)\Delta T)$ are close enough, the trajectory error $e(t)$ will infinitely converge to 0. Here, $f$ is defined as inverse kinematics equation evaluated using the calibrated parameters. One can see that if and only if equation (17) is true, then the trajectory error $e(t)$ is 0. This is illustrated by red points at the intersection of the dashed and solid lines in Fig. 18. However, the theoretical output angles $(\theta _{1}(t),\theta _{2}(t))$ solved by $(\theta _{1}(t),\theta _{2}(t))=f^{-1}(x(t),y(t))$ do not always satisfy equation (17).

(17) \begin{equation} \left\{\begin{array}{l} \theta _{1}\left(t\right)=\theta _{1s}\left(t\right)\\[5pt] \theta _{2}\left(t\right)=\theta _{2s}\left(t\right) \end{array}\right. \end{equation}

Assuming that the maximum allowed value of the trajectory error is $e_{\text{Max}}$ , then the selection of $\Delta T$ should satisfy $\forall t\in [0,T], \exists e(t)\leq e_{\text{Max}}$ to meet the engineering requirement. In practice, the number of pulses sent to drive the servo motor must be an integer, but the values of $(\theta _{1}(t),\theta _{2}(t))$ solved by $(\theta _{1}(t),\theta _{2}(t))=f^{-1}(x(t),y(t))$ are generally not integers. Then, the introduction of a rounding error is inevitable, which can affect the control accuracy of the manipulator [38, 39].

For the ith servo motor, the theoretical input pulse $CD_{L,i}(n)$ between $(n-1)\Delta T$ and $n\Delta T$ can be expressed as:

(18) \begin{equation} CD_{L,i}\left(n\right)=\Delta \theta _{i}/\text{Resi} \end{equation}

Where Resi indicates the resolution of the servo motor, which is related to the incremental encoder, electronic tooth ratio, and reduction ratio.

The theoretical input pulse $CD_{L,i}(n)$ is not always an integer, and the rounding operation is required before being sent to the servo motor. Then, the rounded input pulse $CD_{S,i}(n)$ can be expressed as:

(19) \begin{equation} CD_{S,i}\left(n\right)=\text{round}\left(CD_{L,i}\left(n\right)\right) \end{equation}

When the value of n becomes larger and $\Delta T$ becomes smaller, it is necessary to perform more rounding operations. Ultimately, this will result in a larger rounding error that is introduced into the overall trajectory. In conclusion, the selection of $\Delta T$ should be neither too small nor too large, and should be suitable for different trajectories.

5.2. Influence of rounding error

To verify the above analysis of the rounding error, we next tested a real example. The test trajectory was set as a circle center of (0.3 m, 0 m) and radius of 0.2 m, and the rotational speed was 6 rpm. The electronic tooth ratio of servo motor was set to 20,000, with a 24-bit encoder and reduction ratio of 30. Assuming $\Delta T=0.1\,\text{s}$ and according to the given rotational speed and manipulator parameters, the position coordinates of discrete points $p_{1},p_{2},\cdots,p_{100}$ were obtained. The trajectory errors with and without rounding operation were obtained from equations (14) to (19), which are shown in Fig 19 and Table V.

As shown in Fig. 19(b), the theoretical trajectory error at the red node of $t=\Delta T$ was 0 without rounding operation. Due to the discretization, the theoretical trajectory error around the red node tended to increase and then decrease. This illustrated that it was closely related to the node coordinates and showed a certain regularity. However, after considering the rounding error, the average error and standard deviation increased substantially (Table V). As shown in Fig. 19(a), there was no symmetrical characteristic.

To further analyze the influence of time interval $\Delta T$ on the theoretical trajectory error, the value of $\Delta T$ was set from 0.01 s to 1 s. The effect of $\Delta T$ on the theoretical trajectory error is shown in Fig. 20. In Fig. 20(a) and without rounding operation, the average error decreased with increasing of the number of discrete points. Moreover, it tended to 0. After considering the rounding error, first the average error decreased with increasing of the number of discrete points (Fig. 20(b)). However, when the number of discrete points is higher than 100, the average error fluctuated reciprocally with increasing of the number of discrete points. As shown in Fig. 20(c), this overall tendency also gradually increased. Through the above test example, it is clear that the value of $\Delta T$ should be carefully selected to prevent the value from being too large or too small. This depends on the different manipulator parameters and test trajectories. However, as shown in Table V, the rounding operation does have a greater impact on the trajectory accuracy. This can increase the average error from 0.0631 mm to 0.2023 mm. In fact, the rounding error is a type of systematic error, and one can consider using the feedforward control to eliminate some systematic errors.

Table V. Data statistics of trajectory error.

Figure 19. Trajectory errors with and without rounding operation. (a) Complete figure. (b) Zoom in on the red box of (a).

Figure 20. Effect of $\Delta$ T on theoretical trajectory error. (a) Complete figure. (b) Zoom in on the 1/ $\Delta$ T axis from 0 to 100. (c) Zoom in on the 1/ $\Delta$ T axis from 100 to 1000.

5.3. Design of feedforward controller

Due to the systematic errors (e.g., rounding errors), the errors of repeated trajectory almost remain unchanged across time. In fact, random errors also accounted for some of the trajectory errors. However and in general, systematic errors were the major contributors to trajectory errors. To illustrate the above point, a circular trajectory with a circle center of (0.3 m, 0 m) and a radius of 0.2 m was used as an example for repeated experiments. Two sets of experimental data were selected and the corresponding trajectory errors were obtained (Fig. 21).

Figure 21. Trajectory errors of repeated experiments.

The correlation formula can be expressed as:

(20) \begin{equation} r_{XY}=\frac{\sum \left(X_{i}-\overline{X}_{i}\right)\left(Y_{i}-\overline{Y}_{i}\right)}{\sqrt{\sum \left(X_{i}-\overline{X}_{i}\right)^{2}}\sqrt{\sum \left(Y_{i}-\overline{Y}_{i}\right)^{2}}} \end{equation}

The correlation coefficient $r_{XY}=0.9558$ was obtained from the two sets of experimental data using equation (20), which showed that the two sets of trajectory errors $X_{i}$ and $Y_{i}$ were highly correlated. This high correlation indicates that systematic errors were the major contributors to the trajectory errors.

The trajectory errors can be largely eliminated by feedforward control, and the average errors obtained from the previous repeated experiments were used as the current trajectory error. The specific feedforward control based on kinematic calibration is shown in Fig. 22. As indicated, this is similar to Fig. 10 except for the feedforward control part shown in the red box.

Figure 22. Flowchart of feedforward control based on kinematic calibration.

It is assumed that the theoretical trajectory was $(x(t),y(t))$ , which is discretized into n points by time interval $\Delta T$ . The position coordinates are noted as $\{x_{0}(1),y_{0}(1),\cdots,x_{0}(n),y_{0}(n)\}$ . The actual trajectory data acquired by a laser tracker from the previous repeated experiments are also discretized into n points by time interval $\Delta T$ , and noted as $\{x_{1}(1),y_{1}(1),\cdots,x_{1}(n),y_{1}(n)\}$ . Then, the position coordinates of the actual discrete points $\{x_{11}(1),y_{11}(1),\cdots,x_{11}(n),y_{11}(n)\}$ were obtained after coordinate transformation. The kinematic model after calibration can be expressed as $(x,y)=f_{1}(\theta _{11},\theta _{21})$ , and the initial kinematic model is $(x,y)=f(\theta _{1},\theta _{2})$ . The corresponding inverse functions can be expressed as $(\theta _{11},\theta _{21})=f_{1}^{-1}(x,y)$ and $(\theta _{1},\theta _{2})=f^{-1}(x,y)$ , respectively. Finally, we define the joint input errors as $e_{1}(i)=\theta _{1}(i)-\theta _{11}(i)$ and $e_{2}(i)=\theta _{2}(i)-\theta _{21}(i)$ . Considering the existence of system cumulative errors, the previous cumulative errors needed to be eliminated before proceeding to the feedforward control of the next point, which can be expressed as:

(21) \begin{equation} \left\{\begin{array}{l} e_{11}\left(i\right)=e_{1}\left(i\right)-e_{1}\left(i-1\right)\\[5pt] e_{21}\left(i\right)=e_{2}\left(i\right)-e_{2}\left(i-1\right) \end{array}\right. \end{equation}

Then one can set the feedforward control law as:

(22) \begin{equation} \left(\theta _{12},\theta _{22}\right)=\left(\theta _{11},\theta _{21}\right)-k\left(e_{11},e_{21}\right) \end{equation}

Where $k$ is the feedforward coefficient.

5.4. Experimental study

As mentioned before, the average errors of the previous repeated experiments were sent into the feedforward control part to obtain corrected joint driving data. These data were then fed into the servo motor and reducer by a control card to control the manipulator. The actual trajectory data were collected by a laser tracker. The feedforward control effect with and without calibration is shown in Fig. 23, from which one can see that the trajectory errors at most points were further reduced after combining with the feedforward control. The working range was a circular area with a center of (0.3 m, 0 m) and a radius of 0.3 m.

Figure 23. Feedforward control effect with and without calibration.

However, after combining local calibration and feedforward control, the positioning errors of a very few points were worse than those of the original trajectory (Fig. 23). The spike to large values is at the position coordinate of (0.6, 0), which has reached the boundary of the working range. According to the singular configuration theory of manipulators, this kind of boundary position is prone to give rise to boundary singularity. Therefore, the spike to large values appears at the position coordinate of (0.6, 0), which is mainly the result of control failure caused by boundary singular configuration.

More experimental validations are implemented to illustrate the performances of the combined method of calibration approach and feedforward control within the workspace. The circular paths with the different radius and centers are used in the experiments, such as the radius of 0.2 m and 0.3 m with different centers. As a whole the effects can be improved by the feedforward control and the performance of the combined method will degrade for the circular trajectories with larger radius. However, the situations are not always like this. Because the certainty of error sources in the process, some experiments for the circular trajectory with the radius of 0.3 m will have better performance than that of the radius of 0.2 m. Similarly, the variation of the locations also will have the impact on the effects of the combined method at certain extent.

To further analyze the combined control effect of feedforward control and kinematic calibration, we suggest a concept of feedforward coefficient and calibration coefficient. This approach is similar to the proportional part of PID control. In addition, it is worth noting that the proposed feedforward control can only be applied to the specific trajectory of Fig. 23: Feedforward control effect with and without calibration.

5.5. Analysis of feedforward and calibration coefficients

As discussed above, the influence factors of the combined control effect can be divided into feedforward and calibration coefficients. To facilitate the experimental analysis, a circular trajectory with circle center of (0.3 m, 0 m) and radius of 0.3 m was also selected.

Figure 24. Trajectory errors with different feedforward coefficients. (a) Without calibration. (b) With local calibration.

It should be noticed that the feedforward coefficients used in this process are not normalized. When the feedforward coefficient is set to be 0, the system has no feedforward compensation, and the system has full feedforward compensation when the feedforward coefficient is equal to 1, as shown in Equation (22). Figure 24 shows the trajectory errors obtained for the manipulator with different feedforward coefficients and constant calibration coefficient. When the calibration coefficient was 0, the control effect of feedforward coefficient of 1 was much better than that a value of 0.5 (Fig. 24(a)). Both of these achieved a significant reduction of positioning errors compared with the original errors. As shown in Fig. 24(b) and after combining the local calibration, the control effect of nonzero feedforward coefficient was also better than that of the zero feedforward coefficient. However, it was different from the results in Fig. 24(a). It is noticed that a feedforward coefficient larger than 1 will result in overcompensation for the system. As shown in Fig. 24(b), the control effect of feedforward coefficient of 0.5 was better than that of feedforward coefficient of 1. It can be inferred that the effect with feedforward coefficient larger than 1 will be worse. The trajectory errors of combining the local calibration with feedforward coefficient larger than 1 will be larger than that with zero feedforward coefficient. It is the overcompensation for the control of the system.

However, the control effect of feedforward coefficient of 0.5 was better than that of feedforward coefficient of 1. This was different from the results shown in Fig. 24(a).

In conclusion, the feedforward control had a significant control effect on the trajectory errors. However, the optimal feedforward coefficient varied under different calibration coefficients. To further analyze the coupling effect of feedforward control and kinematics calibration, the test trajectory was also chosen as a circle center of (0.3 m, 0 m) and a radius of 0.3 m. We next experimented with different calibration conditions under feedforward coefficients of 0 and 1, and the experimental results are shown in Fig. 25. Without feedforward control, the trajectory errors of original, global, and local calibration are shown in Fig. 25(a). The trajectory errors of global calibration with feedforward control can be seen from Fig. 25(b). The trajectory errors of local calibration with feedforward control are shown in Fig. 25(c). With feedforward control, the trajectory errors of original, global, and local calibration are shown in Fig. 25(d), which shows that the superposition control effect was not significantly improved compared with the feedforward control.

The average error and standard deviation of different feedforward control and calibration methods are shown in Table VI. When the feedforward control is applied alone, the manipulator was able to reduce most positioning errors, with an average error of 0.2149 mm. After combining the calibration coefficients, the superposition control effects on average errors of 0.2061 mm and 0.2962 mm were not significantly improved compared with the feedforward control alone. Because there were measurement errors that could not be eliminated, the average error obtained (0.2061 mm) may be the best control effect of designed manipulator.

Table VI. Average error and standard deviation.

Figure 25. Trajectory errors of different feedforward control and calibration methods. (a) Without feedforward control (same as Fig. 17). (b) Global calibration with feedforward control. (c) Local calibration with feedforward control. (d) Feedforward control with different calibrations.

It is worth noting that comparing the accuracies of feedforward control and global calibration directly is not appropriate. The global calibration method applies to arbitrary trajectories in the entire workspace, while the feedforward control method is only suitable for specific trajectories.

As demonstrated in Fig. 25(d) and Table VI, pure feedforward control can achieve good performance for specific trajectories but cannot be applied to other trajectories. Similarly, the combination of global calibration and feedforward control method (shown in Table VI) is also limited to the applications of specific trajectories, while that without feedforward control method can be applied to arbitrary trajectories in the entire workspace. The effectiveness of global calibration has been validated through different trajectories, as shown in Fig. 14 and Table II.

It is a more convenient method for global calibration in practical industrial applications, because it can improve the control accuracy for all trajectories with only one calibration process. On the other hand, because of the limitations to specific trajectories, the feedforward control method has significant limitations in practical applications. And the feedforward control method designed in this paper is based on the kinematic calibration method, which is a supplementary method.

6. Conclusions

In this paper, the precise control of a designed three-degree-of-freedom manipulator was taken as the research object. Based on the kinematic model, an error model considering six error sources was established. The global sensitivity analysis of each error source was conducted to distinguish the primary and secondary sources. Then, the ant colony optimization was adopted to identify the actual structural parameters and joint output angles of the manipulator. Experimental results showed that after global calibration, the reduction percentage of the average error was 45.6%.

Considering that there were still some positioning errors after the kinematic calibration, a feedforward control method with different coefficients was next adopted to further improve accuracy. Experimental results showed that the after-feedforward control (without calibration), the percentage reduction of the average error was 77.9%. Then, we combined both the kinematics calibration and feedforward control and the resulting effectiveness was experimentally verified. In conclusion, the proposed feedforward control was limited in that it could be applied to the specific trajectory of a local-not whole workspace. That said, the global calibration method is applicable to the whole workspace. In the future, the machine vision based closed-loop feedback control method will be introduced to further improve the positioning accuracy of the manipulator.