Structural parameters identification for industrial robot using a hybrid algorithm

Abstract

To improve the precision and reduce the movement uncertainty of the industrial robot, a novel hybrid optimization algorithm which combines adaptive genetic algorithm with simulated annealing algorithm is proposed in this article. First, for the sake of increasing the global exploring ability of relevant individuals, the adaptive crossover and mutation operator are used in the phase of adaptive genetic algorithm. If the population optimized by adaptive genetic algorithm is trapped in the local optimal area and simultaneously meets the transformation rule, then it is consequently optimized by simulated annealing to enhance the population diversity and hunt for a better solution so that the probability of finding the global optimal solution is greatly increased. Then, corresponding experiments based on single point repeatability are conducted to acquire data and identify the structural parameters of the industrial robot. Moreover, the single point repeatability test and length test are all implemented at the same time to verify the effectiveness of the proposed method. At last, the result reveals that the proposed method is effective to identify the real structural parameters of the industrial robot, thus enormously decreasing the single point repeatability and length deviation at the same time, which extremely increases the precision and decreases the movement uncertainty of the industrial robot.

Keywords

Industrial robot single point repeatability parameters identification intelligent optimization algorithm

Introduction

With the rapid development of manufacturing and industrial applications, the industrial robot¹ is the universally used apparatuses in various fields, which makes a profound and critical impact in our production and life. For example, it is not strange to find that the industrial robot is widely used in the fields of measurement, assembly, detection, transportation, and so on.² Composed of a base, two arms, three pairs of rotational joints, a probe, and some other relevant components,² one kind of industrial robot illustrated in Figure 1 is a special instrument which makes it able to work expediently, move easily, and display data conveniently.

Figure 1.

The (a) physical structure and (b) coordinate system of the industrial robot.

Although this kind of industrial robot has many merits, but the coordinate of the probe tip modifies in a large scale if the structural parameters change in a small scale resulting from the fact that the coordinate of the probe tip is obtained by a nonlinear transformation function from all the rotational angles of the joints to the probe tip.² Based on this fact and analysis, it is obviously deduced that the coordinate error of the probe tip is greatly and mainly impacted by the structural parameters of the industrial robot.³ However, what should be noticed is that the real values of the structural parameters are not equal to the corresponding nominal ones in almost cases after the industrial robot has been utilized for a period of time as there exists the actual influential factors such as deviations of temperature and humidity, action of gravity, vibrations and interactions between various instruments, and so on.^3
–5 Hence, affirmation can be made that the real coordinate of the probe tip is not commonly equal to the one which is calculated using the nominal structural parameters. Therefore, it is extremely required for us to improve the accuracy of the industrial robot in various fields, to which a great number of researchers devoted much works and energies.

Nowadays, there are two main methods to improve the accuracy of the industrial robot.^6,7 The first one is to use higher precise components such as high precise coders, sensors, arms, and so on, which in turn dramatically increases the manufacturing cost of the industrial robot. The second one is to identify the structural parameters and compensate them again using a little auxiliary instruments or algorithms without any increased cost of components.^7,8

But, because of the nonlinear transformation function between the coordinate of the probe tip and the structural parameters used in the industrial robot,^2

–6 small deviations of the structural parameters may accumulate and lead to relative substantial coordinate errors of the probe tip. Thus, it is obvious that increasing only one or several components’ precision is not capable of effectively decreasing the whole measurement errors in a good degree. As a result, it is essential and necessary for us to execute another method, that is, identify the real structural parameters by other approaches. For example, Fang and Dang⁹ proposed a kinematic parameters identification method using quantum-behaved particle swarm optimization, which is able to simultaneously identify the open-chain series mechanism and parallel mechanism. Santolaria et al.¹⁰ identified the structural parameters of the industrial robot by utilizing four laser tracers to establish a geometric multilateration system and obtained data from laser tracks. Liu et al.¹¹ combined a geometric identification with the parameters identification to propose a new two-step calibration method. He et al.¹² used a laser tracker for spraying robot calibration, where the precision of the robot was improved by 60%. An elasto-geometrical calibration method was used by Kamali et al.¹³ in 2016, in which the geometric parameter errors and joint stiffness parameters were identified at the same time by using a laser tracker.

However, although the methods mentioned above give us a lot of selections and approaches to calibrate and compensate the structural parameters of the industrial robot, but they nearly have the same characteristic that almost all the methods previously mentioned rely on higher precise apparatus especially laser trackers^14,15 which are extreme expensive and affected easily by the external environment, resulting in a relative high cost and nondeterminacy. Moreover, by using lots of precise instruments, the calibration and compensation procedures of the methods are complicated and difficult for us to operate sometimes, which is their defect and bottleneck. To solve this problem and simplify the calibration procedure, a new method is proposed in this article based on single point repeatability by using a hybrid optimization algorithm of adaptive genetic algorithm (AGA)¹⁶ and simulated annealing (SA)^17,18 algorithm. Without using higher precise instruments or complex calibration procedures in the new proposed method, only a conical hole is fixed on a base, to which the probe tip is inserted and the arms of the industrial robot are rotated to a great number of different configurations to generate the acquired data set. Then the data set is utilized by the hybrid algorithm to get the optimized solution, that is, the optimal structural parameters of the industrial robot. And then, the identified parameters of the industrial robot are used to recalculate the single point repeatability and length deviation and some remarkable conclusions are presented at last.

To complete the calibration process, four steps are commonly carried out step by step^19,20: modeling, measurement, identification, and compensation. Modeling is to formulate the transformation function between the coordinate of the probe tip and the structural parameters of the industrial robot. As was presented by Roth et al.,¹⁹ a perfect kinematic model of the robot requires three conditions: (1) completeness, (2) continuousness, and (3) minimality. In recent years, it is known to us all that there are a lot of mathematical models used to calibrate the structural parameters of the industrial robot such as Denavit–Hartenberg (DH) model,^21,22 POE model,²³ CPC model,²⁴ MDH model,²⁵ S model,²⁶ and so on. However, what should be noted is that the DH model and MDH model are most widely utilized in robotic area on account that they have clear physical meaning and thus used conveniently. Therefore, without loss of generality, the MDH model is adopted and used in this article to formulate the transformation function between the coordinate of the probe tip and the structural parameters of the industrial robot.

The main contributions of this article can be briefly summarized as follows:

A new computation mechanism combining AGA with SA to formulate simulated annealing-based adaptive genetic algorithm (AGA-SA) is introduced in this article. Through this integration, the hybrid algorithm AGA-SA not only overcomes the weakness by SA that the AGA is easy to be trapped in the local optimal area but also strengthens the ability of local searching by AGA. That is to say, the AGA-SA is not only able to jump out of the local optimal area when it is trapped in it but also maintains strong local searching ability at the same time, which resulting in a higher probability of not only finding the preciser solution of the structural parameters but also increasing the precision and decreasing the movement uncertainty of the industrial robot.

The kinematic equation is mathematically analyzed in detail and conclusions are derived that it is impossible for us to precisely measure the real coordinates of the probe tip of the industrial robot, which reveals that it consequently nearly impossible for us to directly calculate the real structural parameters by just acquiring the joint angles and coordinate of the industrial robot. Hence, the algorithm AGA-SA and method of single point repeatability are naturally introduced and integrated to optimize the structural parameters and acquire their real ones.

Different from other methods, it is not necessary for us to set precise initial parameters of the structural parameters to be calibrated in the proposed calibration method, that is, the parameters are set arbitrarily as long as they are not too big or small, which simplifies the computation and specification process of the parameters.

The article is organized as follows. The second section formulates the mathematical model of the transformation function from the joint angles to the probe tip. The third section shows the relevant principles and analyses of the scheme. Then, the hybrid algorithm combing AGA and SA is introduced in the fourth section. Moreover, experiments including single point repeatability test and length deviation test are conducted in the fifth section. At last, remarkable conclusions are given in the sixth section.

Kinematic modeling and analysis

Kinematic modeling

The aim of kinematic modeling is to construct the mathematical relationship between the coordinate of the probe tip and the structural parameters of the industrial robot. According to MDH model,²⁵ the coordinate of the probe tip can be computed by equation

\begin{array}{l} P = [\begin{matrix} x \\ y \\ z \\ 1 \end{matrix}] = (\prod_{i = 1}^{6} T_{i - 1, i}) \cdot [\begin{matrix} l_{x} \\ l_{y} \\ l_{z} \\ 1 \end{matrix}] = (\prod_{i = 1}^{6} [\begin{matrix} cos θ_{i} \\ sin θ_{i} cos α_{i - 1} \\ sin θ_{i} sin α_{i - 1} \\ 0 \end{matrix} \begin{matrix} - sin θ_{i} & 0 & a_{i - 1} \\ cos θ_{i} cos α_{i - 1} & - sin α_{i - 1} & - d_{i} sin α_{i - 1} \\ cos θ_{i} sin α_{i - 1} & cos α_{i - 1} & d_{i} cos α_{i - 1} \\ 0 & 0 & 1 \end{matrix}]) \cdot [\begin{matrix} l_{x} \\ l_{y} \\ l_{z} \\ 1 \end{matrix}] \end{array}

where $θ_{i}$ is the joint angle, d_i is the linkage length, $a_{i - 1}$ is the joint offset, $α_{i - 1}$ is the torsion angle, and l_x , l_y , and l_z are the components of the probe length in x direction, y direction, and z direction, respectively. And, x, y, and z are the coordinates of the probe tip in x direction, y direction, and z direction, respectively.

It should be noted that there is a deviation between the zero position of the coder and that of the joint, so the joint angle $θ_{i}$ is expressed by equation

θ_{i} = θ_{0, i} + θ_{s, i}

where $θ_{0, i}$ is the deviation of the zero position and $θ_{s, i}$ is the rotational angle of the coder. Note that $θ_{0, i}$ is a constant and a structural parameter. The nominal structural parameters of the industrial robot used in this article are shown in Table 1.

Table 1.

The nominal structural parameters of the industrial robot.

Joint no.	$θ_{0, i} (^{\circ})$	$d_{i} (m m)$	$a_{i - 1} (m m)$	$α_{i - 1} (^{\circ})$
1	0	147.5	0	0
2	0	156	0	−90
3	0	112	395	0
4	0	127.5	375	0
5	0	130	0	−90
6	0	0	0	90
$P r o b e b i a s (m m)$	$l_{x} = 0, l_{y} = 0, l_{z} = 47$

Kinematic analysis

Set $X = (θ_{0, 1}, θ_{0, 2}, θ_{0, 3}, θ_{0, 4}, θ_{0, 5}, θ_{0, 6}, d_{1}, d_{2}, d_{3}, d_{4}, d_{5}, d_{6}, a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, α_{0}, α_{1}, α_{2}, α_{3}, α_{4}, α_{5}, l_{x}, l_{y}, l_{z})$ , then the coordinate $(x, y, z)$ of the probe tip can be rewritten as the structural parameter X and the readings of the coder by equation

[\begin{matrix} x \\ y \\ z \end{matrix}] = [\begin{matrix} f_{x} (θ_{s, 1}, θ_{s, 2}, θ_{s, 3}, θ_{s, 4}, θ_{s, 5}, θ_{s, 6}, X) \\ f_{y} (θ_{s, 1}, θ_{s, 2}, θ_{s, 3}, θ_{s, 4}, θ_{s, 5}, θ_{s, 6}, X) \\ f_{z} (θ_{s, 1}, θ_{s, 2}, θ_{s, 3}, θ_{s, 4}, θ_{s, 5}, θ_{s, 6}, X) \end{matrix}]

The total differential with respect to every structural parameter of equation (3) is formulated by equation

[\begin{matrix} d x \\ d y \\ d z \end{matrix}] = [\begin{matrix} \sum_{i = 1}^{6} \frac{\partial f_{x}}{\partial θ_{0, i}} \cdot d θ_{0, i} + \sum_{i = 1}^{6} \frac{\partial f_{x}}{\partial d_{i}} \cdot d d_{i} + \sum_{i = 1}^{6} \frac{\partial f_{x}}{\partial a_{i - 1}} \cdot d a_{i - 1} \\ \sum_{i = 1}^{6} \frac{\partial f_{y}}{\partial θ_{0, i}} \cdot d θ_{0, i} + \sum_{i = 1}^{6} \frac{\partial f_{y}}{\partial d_{i}} \cdot d d_{i} + \sum_{i = 1}^{6} \frac{\partial f_{y}}{\partial a_{i - 1}} \cdot d a_{i - 1} \\ \sum_{i = 1}^{6} \frac{\partial f_{z}}{\partial θ_{0, i}} \cdot d θ_{0, i} + \sum_{i = 1}^{6} \frac{\partial f_{z}}{\partial d_{i}} \cdot d d_{i} + \sum_{i = 1}^{6} \frac{\partial f_{z}}{\partial a_{i - 1}} \cdot d a_{i - 1} \end{matrix} \begin{matrix} + \sum_{i = 1}^{6} \frac{\partial f_{x}}{\partial α_{i - 1}} \cdot d α_{i - 1} + \frac{\partial f_{x}}{\partial l_{x}} + \frac{\partial f_{x}}{\partial l_{y}} + \frac{\partial f_{x}}{\partial l_{z}} \\ + \sum_{i = 1}^{6} \frac{\partial f_{y}}{\partial α_{i - 1}} \cdot d α_{i - 1} + \frac{\partial f_{y}}{\partial l_{x}} + \frac{\partial f_{y}}{\partial l_{y}} + \frac{\partial f_{z}}{\partial l_{z}} \\ + \sum_{i = 1}^{6} \frac{\partial f_{z}}{\partial α_{i - 1}} \cdot d α_{i - 1} + \frac{\partial f_{z}}{\partial l_{x}} + \frac{\partial f_{z}}{\partial l_{y}} + \frac{\partial f_{z}}{\partial l_{z}} \end{matrix}]

If the structural parameters of the industrial robot change in an extreme small scale, then the differentials of the structural parameters and the coordinates of the probe tip are referred as equal to their corresponding deviations, that is

\begin{array}{l} [d θ_{0, 1}, d θ_{0, 2}, d θ_{0, 3}, d θ_{0, 4}, d θ_{0, 5}, d θ_{0, 6}, d d_{1}, d d_{2}, d d_{3}, \\ d d_{4}, d d_{5}, d d_{6}, d a_{0}, d a_{1}, d a_{2}, d a_{3}, d a_{4}, d a_{5}, d α_{0}, d α_{1}, \\ d α_{2}, d α_{3}, d α_{4}, d α_{5}, d l_{x}, d l_{y}, d l_{z}] = \\ [Δ θ_{0, 1}, Δ θ_{0, 2}, Δ θ_{0, 3}, Δ θ_{0, 4}, Δ θ_{0, 5}, Δ θ_{0, 6}, Δ d_{1}, Δ d_{2}, \\ Δ d_{3}, Δ d_{4}, Δ d_{5}, Δ d_{6}, Δ a_{0}, Δ a_{1}, Δ a_{2}, Δ a_{3}, Δ a_{4}, Δ a_{5}, \\ Δ α_{0}, Δ α_{1}, Δ α_{2}, Δ α_{3}, Δ α_{4}, Δ α_{5}, Δ l_{x}, Δ l_{y}, Δ l_{z}] \end{array}

and

[d x d y d z] = [Δ x Δ y Δ z]

On the condition that the nominal structural parameters are given beforehand and 27 structural parameters need to be calibrated, if 9 groups of real coordinates of the probe tip and corresponding joint angles are obtained by some instruments, then there are 27 groups of equations derived by equations (4), (5), and (6), which can help to not only directly calculate the deviations of the structural parameters but also compensate the nominal structural parameters by these deviations in turn.

Principle and analysis

Analysis

Based on the analyses previously mentioned, the coordinates of the probe tip can be obtained and the structural parameters can be calibrated by equation (3). Specifically speaking, for a given X, equations can be formulated from equation (4) by using higher precise instruments to acquire N groups of different points ${}^{k}P ({}^{k}x, {}^{k}y, {}^{k}z) (k = 1, 2, \dots, N)$ which are mathematically expressed by equation

[\begin{matrix} {}^{k}x \\ {}^{k}y \\ {}^{k}z \end{matrix}] = [\begin{matrix} f_{x} ({}^{k}θ_{s, 1}, {}^{k}θ_{s, 2}, {}^{k}θ_{s, 3}, {}^{k}θ_{s, 4}, {}^{k}θ_{s, 5}, {}^{k}θ_{s, 6}, X) \\ f_{y} ({}^{k}θ_{s, 1}, {}^{k}θ_{s, 2}, {}^{k}θ_{s, 3}, {}^{k}θ_{s, 4}, {}^{k}θ_{s, 5}, {}^{k}θ_{s, 6}, X) \\ f_{z} ({}^{k}θ_{s, 1}, {}^{k}θ_{s, 2}, {}^{k}θ_{s, 3}, {}^{k}θ_{s, 4}, {}^{k}θ_{s, 5}, {}^{k}θ_{s, 6}, X) \end{matrix}]

where $k = 1, 2, \dots, N$ . Here, index k indicates the k groups of coordinates of the probe tip. By solving equation (7), the real structural parameters are obtained and can consequently be used as the new nominal parameters of the industrial robot.

However, the fact that the general coordinates ${}^{k}P ({}^{k}x, {}^{k}y, {}^{k}z) (k = 1, 2, \dots, N)$ in equation (7) are measured by higher precise instruments makes it difficult even impossible for us to implement the process in real fields. Therefore, it is extreme necessary for us to use other ways to settle this problem. One method to solve this difficulty and simplify this process is to use the single point repeatability calibration. For this process, a single point is measured for a great number of times with different poses for obtaining the corresponding joint angles and calibrating the real structural parameters, which is the main principle of this method.

Objective function

As mentioned above, after fixing the conical hole on the base, the probe tip is inserted to the conical hole still and the arms are rotated to N kinds of different configurations for obtaining N groups of different joint angles ${}^{k}θ_{i} (i = 1, 2, ..., 6; k = 1, 2, ..., N)$ . Suppose an arbitrarily given structural parameter X is known, then a series of points ${}^{k}P ({}^{k}x, {}^{k}y, {}^{k}z) (k = 1, 2, ..., N)$ can be computed by substituting ${}^{k}θ_{i} (i = 1, 2, ..., 6; k = 1, 2, ..., N)$ and X into equation (7) and the deviations between the measured coordinates and the real coordinates are thus calculated by equation

{}^{k}E = \sqrt{{({}^{k}x - x)}^{2} + {({}^{k}y - y)}^{2} + {({}^{k}z - z)}^{2}}

where x, y, and z are the real coordinates of the probe tip in direction of x, y, and z, respectively. Note that it is not easy for us to measure the real coordinate $(x, y, z)$ of the probe tip especially in real fields, hence the coordinate $(x, y, z)$ is replaced by its mean value expressed by equation

{\begin{matrix} x = \frac{1}{N} \sum_{i = 1}^{N} {}^{k}x \\ y = \frac{1}{N} \sum_{i = 1}^{N} {}^{k}y \\ z = \frac{1}{N} \sum_{i = 1}^{N} {}^{k}z \end{matrix}

After given ${}^{k}E$ , then their average is

\bar{E} = \frac{1}{N} \sum_{k = 1}^{N} {}^{k}E

However, the average deviation $\bar{E}$ cannot indicate the distribution regularity of this group of data, so the standard deviation $σ$ is introduced by equation

σ = \sqrt{\frac{1}{N - 1} \sum_{i = 1}^{N} {({}^{k}E - \bar{E})}^{2}}

Consequently, the objective function for X based on “3 sigma principle” is

F (X) = \bar{E} + 3 σ

Let $X^{*}$ be the real structural parameter of the industrial robot, then

F (X^{*}) = \bar{E^{*}} + 3 σ^{*} = 0

Equation (12) is utilized as the objective function to search for the global optimal solution which is accepted as the real structural parameter of the industrial robot and compensates them in turn. Note that the closer X approaches to $X^{*}$ , the smaller $F (X)$ is. Therefore, the main work thereafter is to obtain the different joint angles and use them to search for the global optimal solution by usage of optimization algorithm.

Simulated annealing-based adaptive genetic algorithm

Adaptive genetic algorithm

Genetic algorithm (GA),²⁷ a natural bioinspired intelligent optimization algorithm simulating the biological mechanism of selection, crossover, and mutation, is a widely used evolutionary optimization method in various fields, which obeys the principle of “survival of the fitness” proposed by Charles Robert Darwin. In GA, the individuals’ survival abilities are indicated by their own fitness values, that is, the individuals with larger fitness values have better survival abilities and thus have larger probabilities of executing the genetic operation and giving birth to child. Note that the individuals with larger fitness values have better survival abilities while the aim of this problem in this article is to search for the best solution with the minimal objective value, hence the fitness function takes the form of

f i t n e s s (X) = \frac{1}{F (X) + r e a l lim}

where $f i t n e s s (\cdot)$ is the fitness value of the corresponding variable, $r e a l lim$ is an extreme small real constant which ensures that the denominator is not equal to zero when $F (X)$ is equal or extremely approaches to zero.

However, although GA has a good global searching ability, it can be trapped in a local optimal area with a certain probability especially when the optimized problem is extreme nonlinear such as the problem in this article. Therefore, the AGA¹⁶ is used in this article to increase the probability of hunting for the global optimum. In AGA, the probability of crossover and mutation are all self-adaptive dynamically changing by taking the form of

P_{c} = {\begin{cases} k_{1} (f_{max} - f^{'}) / (f_{max} - \bar{f}), f^{'} \geq \bar{f} \\ k_{3}, f^{'} < \bar{f} \end{cases}

and

P_{m} = {\begin{cases} k_{2} (f_{max} - f) / (f_{max} - \bar{f}), f \geq \bar{f} \\ k_{4}, f < \bar{f} \end{cases}

where $k_{1}, k_{2}, k_{3}, k_{4} \leq 1$ are constants, $f_{max}$ , $\bar{f}$ , $f^{'}$ , f are the maximal fitness value, the mean fitness value, the larger fitness value of the two crossover individuals, and the mutation fitness value, respectively. In this article, $k_{1}, k_{2}, k_{3}, k_{4}$ are set to be 1, 0.5, 1, and 0.5, respectively.

For clarification, the details about AGA are specified as follows.

Ⅰ.

Selection. In selection phase, the individuals are selected by roulette wheel²⁸.

Ⅱ.

Crossover. Let $X_{i} (t)$ and $X_{j} (t)$ be any two arbitrary individuals in the tth generation, then their corresponding newly generated individuals in the $(t + 1)$ th generation after executed the operation of crossover are updated by equation

X_{i} (t + 1) = r X_{i} (t) + (1 - r) X_{j} (t)

and

X_{j} (t + 1) = r X_{j} (t) + (1 - r) X_{i} (t)

where r is a constant given beforehand.

Ⅲ.Mutation. Mutation operation is executed successively as the following steps:

①Generate the mutation position $u = r a n d i (D)$ , where D is the dimension of the variable X and $r a n d i$ is a system function which can randomly selects an integer from integer array varying from 1 to D wherein this article $D = 27$ .

②Define variable $a = X_l o w e r limit (u)$ , $b = X_u p p e r limit (u)$ , where $X_l o w e r limit$ and $X_u p p e r limit$ are the minimal limit vector and maximal limit vector of variable X, respectively.

③Let $s (t) = 1 - {(\frac{t}{G})}^{b}$ , where t is the current generation number, G and b are the maximal generation number and a constant, respectively, which are all given beforehand.

④Let $X_{i, u} (t)$ and $X_{i, u} (t + 1)$ be the values in the $u$ th position of the tth and $(t + 1)$ th generation, respectively, of any arbitrary individual i; $X_{i, u} {(t + 1)}_{l o w e r}$ and $X_{i, u} {(t + 1)}_{u p p e r}$ be the lower value and upper value of $X_{i, u} (t + 1)$ , respectively, then

X_{i, u} {(t + 1)}_{l o w e r} = X_{i, u} (t) - s (t) (X_{i, u} (t) - a)

and

X_{i, u} {(t + 1)}_{u p p e r} = X_{i, u} (t) + s (t) (b - X_{i, u} (t))

⑤At last, the value $X_{i, u} (t + 1)$ in the $u$ th position of $X_{i} (t + 1)$ is updated by equation

X_{i, u} (t + 1) = X_{i, u} {(t + 1)}_{l o w e r} + r a n d \cdot (X_{i, u} {(t + 1)}_{u p p e r} - X_{i, u} {(t + 1)}_{l o w e r})

where $r a n d$ is a random function arbitrarily generating a number from 0 to 1.

SA algorithm and AGA-SA

SA^17,18 algorithm, which simulates the annealing process of the metal, has a good characteristic of escaping from the local optimal area when the best solution found so far is trapped in the local optimal region. Motivated by this trait, if the solution optimized by AGA is trapped in the local optimal region, then it is optimized by SA through a transformation rule to help the current best individual to escape from the local optimal region and wander widely in a broader space. Let t, $X_{b e s t} (t)$ and $f i t n e s s (X_{b e s t} (t))$ be the current generation number, the best individual found so far until tth generation, and the fitness value of the best individual $X_{b e s t} (t)$ , respectively, then the transformation rule takes the form of

{\begin{cases} | f i t n e s s (X_{b e s t} (t)) - f i t n e s s (X_{b e s t} (t - K)) | \leq ε, \\ f i t n e s s (X_{b e s t} (t)) \leq δ, t > K \geq 1 \end{cases}

where K is a positive integer, $ε$ and $δ$ are transformational thresholds. Note that K, $ε$ , and $δ$ are constants and set by ourselves beforehand. That is to say, when the solution $X_{b e s t} (t)$ optimized by AGA meets equation (22), then the optimization process enters the SA process, that is, the best individual $X_{b e s t} (t)$ optimized by AGA is optimized by SA furtherly. Note that if the best solution $X_{b e s t} (t)$ found so far is not the real parameter until AGA-SA runs a great number of generations and the best solution in the last generation is nearly kept still, then the whole optimization process cannot find the global optimum. Thus, a termination condition is set to finish the unnecessary redundant computation of AGA-SA if the optimized result meets the accuracy demand. In this article, the termination condition is set as

F (X_{b e s t} (t)) \leq η

f i t n e s s (X_{b e s t} (t)) \geq τ

t \geq G

where t is the generation number, G is the maximal generation number, $τ$ and $η$ are thresholds which determine the ending of the whole optimization process.

The relative parameters are set as follows:

①Population number $N = 60$ ; $r = 0.6$ ; $b = 3$ ; $η = 0.1$ .

② $τ = 80$ ; G = 35,000; $K = 25$ ; $ε = 0.1$ ; $δ = 2$ ; $T_{0} = 10^{4}$ ; $J = 30$ ; $T_{1} = 10^{- 3}$ ; $μ = 0.9$ .

To make comparison, the AGA, SA, and AGA-SA are respectively utilized as the algorithm to get the structural parameters of the industrial robot and their corresponding obtained optimized structural parameters are respectively implanted into the corresponding controller to replace the corresponding old ones. Moreover, the industrial robot with the old structural parameters and the corresponding new structural parameters is reused again to conduct the experiments of single point repeatability and length deviation test and the relevant outcomes are respectively presented thereafter. The flow chart of AGA, SA, and AGA-SA for calibrating the structural parameters of the industrial robot are shown in Figures 2, 3, and 4, respectively.

Figure 2.

The flow chart of AGA. AGA: adaptive genetic algorithm.

Figure 3.

The flow chart of SA. SA: simulated annealing.

Figure 4.

The flow chart of AGA-SA. AGA-SA: simulated annealing-based adaptive genetic algorithm.

Experiments

Data acquisition

As is illustrated in Figure 5, after the conical hole and the base of the industrial robot are simultaneously fixed on the platform firstly, then the arms of the industrial robot are manipulated to make the probe tip inserted to the conical hole and fitted with it. Moreover, the arms are rotated again to N groups of different poses and the readings ${}^{k}θ_{i} (i = 1, 2, \dots, 6; k = 1, 2, \dots, N)$ of the encoder at each position are recorded from the man–machine interface illustrated in Figure 6 while maintaining the probe tip fitted with the conical hole during the whole process. The data set ${}^{k}θ_{i} (i = 1, 2, \dots, 6; k = 1, 2, \dots, N)$ is then acquired and utilized thereafter. Note that the top and bottom diameter of the conical hole are a little bigger and smaller than that of the probe tip, respectively. Therefore, the probe tip can be trapped and fitted in the conical hole to maintain the spatial position of the center of the probe tip unchanged when the arms are rotated to different poses.

Figure 5.

The experimental test of the single point repeatability. (a) and (b) are global figure of the single point repeatability test and partial enlarged figure of A, respectively.

Figure 6.

The man–machine interface.

Parameters identification

After obtaining N groups of joint angles, they are used as the data set by AGA, SA, and AGA-SA, respectively, to optimize the structural parameters. If a certain best individual $X_{b e s t} (t)$ satisfies equation (23) or equation (24), then it is acknowledged that the algorithm successfully finds the optimal solution $X_{b e s t} (t)$ denoted by $X^{o p t i m u m}$ which is accepted as the real parameter of the industrial robot. Note that $X^{o p t i m u m}$ which is derived by relevant algorithm and test data is usually infinitely close but not equal to $X^{*}$ due to the existence of errors of real engineering problems.

Other critical problems are the range settings of the various variables of the industrial robot and their boundary treatments when the variables exceed the corresponding boundaries. In this article, because the values of the length parameters are nonnegative numbers and the nominal values of $a_{i - 1} (i = 1, 2, 5, 6)$ , d ₆, and $α_{i - 1} (i = 1, 3, 4)$ are all 0, so the deviation ranges of $a_{i - 1} (i = 1, 2, 5, 6)$ and d ₆ are set to be (0, 10); of $a_{i - 1} (i = 3, 4)$ and $d_{i} (i = 1, 2, ..., 5)$ are set to be ±10; of $α_{i - 1} (i = 1, 2, ..., 6)$ is set to be ±10°; of $θ_{0, i} (i = 1, 2, ..., 6)$ is set to be ±1°; and of l_x , l_y , and l_z are set to be ±10, respectively.

When a certain variable y among the structural parameters surpasses the corresponding limit, then it is newly generated randomly by a generator, which is expressed by Criterion A illustrated as:

Criterion A

If y exceeds the corresponding constraint $y = (u p lim i t - l o w e r lim i t) . * r a n d + y$ end

where $u p limit$ and $l o w e r limit$ are the upper boundary and lower boundary of variable y, respectively.

Note that the closer the any two arbitrary individuals X_i and X_j are, the smaller of the gap between their fitness values are. That is to say, if ${| X_{i} - X_{j} |}_{2} < {| X_{i} - X_{k} |}_{2}$ , then $| f i t n e s s (X_{i}) - f i t n e s s (X_{j}) | < | f i t n e s s (X_{i}) - f i t n e s s (X_{k}) |$ , where ${| \cdot |}_{2}$ indicates 2-norm of the corresponding variable. Based on this fact, the term $f i t n e s s_{max} (t) - f i t n e s s_{a v e r a g e} (t)$ is introduced in this article to indicate the individual density of the population in the tth generation, where $f i t n e s s_{max} (t)$ and $f i t n e s s_{a v e r a g e} (t)$ are the maximal fitness value and average fitness value in the tth generation, respectively. The changing curves of the objective function and $f i t n e s s_{max} (t) - f i t n e s s_{a v e r a g e} (t)$ optimized by AGA-SA are shown in Figures 7 and 8, respectively. The calibrated structural parameters of the industrial robot optimized by AGA, SA, and AGA-SA are shown in Tables 2, 3, and 4, respectively.

Figure 7.

Evolutionary curve of the objective function optimized by AGA-SA. AGA-SA: simulated annealing-based adaptive genetic algorithm.

Figure 8.

Evolutionary curve of $f i t n e s s_{max} (t) - f i t n e s s_{a v e r a g e} (t)$ optimized by AGA-SA. AGA-SA: simulated annealing-based adaptive genetic algorithm.

Table 2.

The calibrated structural parameters of the industrial robot with respect to AGA.

Joint no.	$θ_{0, i} (^{\circ})$	$d_{i} (m m)$	$a_{i - 1} (m m)$	$α_{i - 1} (^{\circ})$
1	0.0715	147.4267	0.2543	0.1486
2	0.0368	156.7328	0.0863	−89.5327
3	0.1287	112.0745	395.4387	−0.4005
4	−0.0845	127.7328	374.8456	0.9362
5	−0.0046	129.5328	0.3275	−90.3645
6	0.0403	0.1053	0.0743	90.1563
Probe bias (mm)	$l_{x} = 0.3456$ , $l_{y} = 0.5974$ , $l_{z} = 46.2367$

AGA: adaptive genetic algorithm.

Table 3.

The calibrated structural parameters of the industrial robot with respect to SA.

Joint no.	$θ_{0, i} (^{\circ})$	$d_{i} (m m)$	$a_{i - 1} (m m)$	$α_{i - 1} (^{\circ})$
1	−0.0847	147.4407	0.3063	0.1732
2	−0.0375	156.4382	0.1453	−90.4328
3	0.1386	111.0653	394.4367	0.2873
4	0.0763	127.1432	374.7864	0.8563
5	0.0028	130.3296	0.2946	−89.5368
6	0.0354	0.0872	0.0943	89.6847
Probe bias (mm)	$l_{x} = 0.3653$ , $l_{y} = 0.6831$ , $l_{z} = 46.1532$

SA: simulated annealing.

Table 4.

The calibrated structural parameters of the industrial robot with respect to AGA-SA.

Joint no.	$θ_{0, i} (^{\circ})$	$d_{i} (m m)$	$a_{i - 1} (m m)$	$α_{i - 1} (^{\circ})$
1	−0.0632	147.5328	0.2987	0.1572
2	0.0482	155.6032	0.0954	−89.4632
3	0.1054	112.0564	394.6532	−0.3045
4	−0.0963	127.2465	375.1678	0.8023
5	−0.0032	129.6423	0.2658	−90.2865
6	0.0354	0.0984	0.0845	90.2789
Probe bias (mm)	$l_{x} = 0.2384$ , $l_{y} = - 0.6523$ , $l_{z} = 47.8205$

AGA-SA: simulated annealing-based adaptive genetic algorithm.

Verification of the calibrated parameters

Because it is difficult for us to measure the real structural parameters of the industrial robot, thus the single point repeatability test illustrated in Figure 5 and the length test illustrated in Figure 9 are simultaneously used in this phase to verify the effectiveness and correctness of the method used in this article. As is shown in Figure 5, the single point is tested for many times and repeatability results with respect to AGA, SA, and AGA-SA are shown in Figure 10.

Figure 9.

The experimental test of the standard distance L between the (a) left and (b) right standard balls.

Figure 10.

The changing curves of the single point repeatability. (a) is the changing curve of single point repeatability before calibration and (b), (c), and (d) are the changing curves of single point repeatability after calibration with respect to AGA, SA, and AGA-SA, respectively. AGA: adaptive genetic algorithm; SA: simulated annealing; AGA-SA: simulated annealing-based adaptive genetic algorithm.

As is shown in Figure 9, the two standard balls are attached to the fixture plate with a standard length L between the two ball centers $O_{1}$ and $O_{2}$ and the fixture plate is fixed by a vise on the platform. In this test, the arms of the industrial robot are rotated to make the ball of the probe tip attached to the standard balls in different positions to obtain different groups of joint angles, with which the length L between the two standard balls is computed through the transformation function. For the sake of narrative, the error $Δ L$ between the standard length L and the tested length $L^{'}$ is defined as

Δ L = | L - L^{'} |

Let ${L'}_{1}$ and ${L'}_{2}$ be the tested length using the nominal structural parameters and identified structural parameters of the industrial robot, respectively, then

Δ L_{1} = | L - {L'}_{1} |

and

Δ L_{2} = | L - {L'}_{2} |

For comparison, the standard length L is tested for many times and the corresponding error curves are shown in Figure 11. Note that the radiuses of the two standard balls are all R and the length between the two centers $O_{1}$ and $O_{2}$ of the corresponding standard balls is a fixed distance L, which is illustrated in Figure 9(b).

Figure 11.

The changing curve of the deviation of the standard length L between the two standard balls.

Discussion

After the experiments are completed and verifications are accomplished which are illustrated by the tables and figures, observations can be made and conclusions can be drawn consequently, which are divided into Part 1 and Part 2, respectively.

Part 1. The observations and conclusions about single point repeatability

Before calibration, it is observed from Figure 10(a) that:

①The maximal (Max), minimal (Min), mean (Mea), and standard deviations (Std) of the single point repeatability before calibration are 4.5898 mm, 0.3409 mm, 2.3534 mm, and 0.9331 mm, respectively.

②The Max, Min, Mea, and Std of the single point repeatability before calibration in x direction are 3.3673 mm, −3.1302 mm, 0.020606 mm, and 2.0215 mm, respectively; in y direction are 4.1769 mm, −1.3273 mm, 0.44338 mm, and 1.2210 mm, respectively; and in z direction are 1.7511 mm, −1.4132 mm, −0.17053 mm, and 0.9722 mm, respectively.

However, after optimized by AGA, Figure 10(b) tells us that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration are 0.2344 mm, 0.0630 mm, 0.1438 mm, and 0.0384 mm, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration in x direction are 0.1395 mm, −0.1436 mm, 0.0164 mm, and 0.0853 mm, respectively; in y direction are 0.1490 mm, −0.1482 mm, −0.0101 mm, and 0.0890 mm, respectively; and in z direction are 0.1381 mm, −0.1468 mm, −0.0190 mm, and 0.0815 mm, respectively.

After optimized by SA, Figure 10(c) tells us that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration are 0.4224 mm, 0.0675 mm, 0.2433 mm, and 0.0950 mm, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration in x direction are 0.2850 mm, −0.1865 mm, 0.0133 mm, and 0.1371 mm, respectively; in y direction are 0.2999 mm, −0.1924 mm, 0.0645 mm, and 0.1501 mm, respectively; and in z direction are 0.2758 mm, −0.1920 mm, 0.0551 mm, and 0.1337 mm, respectively.

After optimized by AGA-SA, it is derived from Figure10(d) that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration are 0.1125 mm, 0.0063 mm, 0.0833 mm, and 0.0189 mm, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration in x direction are 0.0794 mm, −0.0737 mm, 0.0017 mm, and 0.0516 mm, respectively; in y direction are 0.0698 mm, −0.0724 mm, −0.0075 mm, and 0.0454 mm, respectively; and in z direction are 0.0787 mm, −0.0766 mm, −0.0039 mm, and 0.0513 mm, respectively.

From the results presented above, analyses can be conducted and some conclusions can be acquired about single point repeatability as follows:

After optimized by AGA, it is concluded that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration have been decreased by 95%, 81%, 94%, and 95%, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration in x direction are decreased by 96%, 95%, 20%, and 96%, respectively; in y direction are decreased by 96%, 89%, 77%, and 92%, respectively; and in z direction are decreased by 92%, 89%, 89%, and 92%, respectively.

After optimized by SA, conclusions are made that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration have been decreased by 91%, 80%, 90%, and 89%, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration in x direction are decreased by 92%, 94%, 35%, and 93%, respectively; in y direction are decreased by 93%, 86%, 85%, and 88%, respectively; and in z direction are decreased by 84%, 86%, 68%, and 86%, respectively.

After optimized by AGA-SA, it is observed that:

①The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration have been decreased by 98%, 98%, 96%, and 98%, respectively.

②The Max, Min, Mea, and Std of the single point repeatability after calibration compared with those before calibration in x direction are decreased by 98%, 98%, 92%, and 98%, respectively; in y direction are decreased by 98%, 95%, 83%, and 96%, respectively; and in z direction are decreased by 96%, 95%, 98%, and 95%, respectively.

Part 2. The observations and conclusions about length error

Before calibration, it is observed from Figure 11 that the Max, Min, Mea, and Std of the length error $Δ L_{1}$ are 0.3757 mm, −0.2773 mm, 0.0478 mm, and 0.2024 mm, respectively.

After optimized by AGA, SA, and AGA-SA, it is derived from Figure 11 that the Max, Min, Mea, and Std of $Δ L_{2}$ after optimized by AGA are 0.0985 mm, −0.0986 mm, 0.0077 mm, and 0.0575 mm, respectively; after optimized by SA are 0.1988 mm, −0.1940 mm, −0.0093 mm, and 0.1171 mm, respectively; and after optimized by AGA-SA are 0.0493 mm, −0.0499 mm, −0.0080 mm, and 0.0289 mm, respectively.

From the results presented above, analyses can be conducted and some conclusions can be acquired that the Max, Min, Mea, and Std of $Δ L_{2}$ optimized by AGA with respect to $Δ L_{1}$ are decreased by 74%, 64%, 84%, and 72%, respectively; optimized by SA with respect to $Δ L_{1}$ are decreased by 44%, 30%, 81%, and 42%, respectively; and optimized by AGA-SA with respect to $Δ L_{1}$ are decreased by 86%, 82%, 83%, and 86%, respectively.

Overall, from the observations and analysis presented above, it is evidently derived that:

It is nearly impossible for us to conduct the method of measuring the real coordinates of the probe tip of the industrial robot and directly computing the real structural parameters using the tested data in the real fields, thus a method that combines testing the rotational angles of the probe tip of the industrial robot and utilizing the hybrid algorithm AGA-SA is proposed to simplify the process and resolve the contradictions.

The Max, Min, Mea, and Std of the single point repeatability and length error, including their three components in the corresponding directions, are all decreased with a big scale after optimized by corresponding algorithms (AGA, SA, and AGA-SA). However, nearly all the Max, Min, Mea, and Std of the single point repeatability and length error are decreased with a greater scale after optimized by AGA-SA than those optimized by AGA or SA at the same time.

Through this analysis, it is shown that the precision and relevant errors of the industrial robot are respectively increased and decreased drastically with a bigger corresponding increasing and decreasing magnitude after optimized by AGA-SA than those optimized by AGA and SA, which indicates that this proposed method and algorithm are effective and more suitable for us to calibrate the real structural parameters of the industrial robot than AGA and SA simultaneously.

Conclusions

A kinematic model of the industrial robot is established based on MDH model firstly. Then, the kinematic model and error formula are analyzed and results are derived that: (1) it is extreme difficult even impossible for us to measure the real coordinates of the industrial robot and; (2) only using a simple conventional algorithm is easy to be trapped in the local optimal region and thus incapable of searching for the real structural parameters of the industrial robot. Moreover, based on these facts and to solve the problem presented previously, the algorithm of AGA and SA are combined to formulate a hybrid algorithm AGA-SA which simultaneously has a good local searing ability and global exploring capacity. And then, the experiments of the single point repeatability and the length test are conducted combining the hybrid algorithm to calibrate the structural parameters and verify the effectiveness of the proposed method. The result shows that the single point repeatability and the length deviation after calibrated by AGA-SA are all decreased substantially and with a bigger magnitude compared with those calibrated by AGA or SA, which suggests that this method is effective and more suitable for us to greatly increase the precision and decrease the movement uncertainty of the industrial robot than AGA and SA.

The proposed modeling and identification method can be commonly expanded to the calibration of a serials of robots, where the identification method is required to increase the accuracy of the corresponding robots. This will be furtherly studied in our future work.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Kejin Liu

References

Zurong

Jie

. Review of performance testing of high precision reducers for industrial robots. Measurement 2021; 183: 109794.

Zhibin

. An overview of calibration technology of industrial robots. IEEE/CAA J Autom Sinica 2021; 8: 27–40.

Gao

Sun

Jing

, et al. Structural parameter identification for 6 DOF industrial robots. Mech Syst Signal Process 2017; 113: 145–155.

Gao

Zhao

. Decoupling of kinematic parameter identification for articulated arm coordinate measuring machines. IEEE Access 2018; 6: 1.

Arai

AWAS

. A kinematic calibration method for industrial robots using autonomous visual measurement. CIRP Ann Manuf Technol 2006; 55: 1–6.

Gao

Zhang

San

, et al. Kinematic calibration for industrial robots using articulated arm coordinate machines. Int J Model Identif Contr 2019; 31: 16.

Guo

, et al. A self-calibration method for articulated arm coordinate measuring machines. T Can Soc Mech Eng 2016; 40(4): 645–655.

Filion

Joubair

Tahan

, et al. Robot calibration using a portable photogrammetry system. Robot Cim-Int Manuf 2018; 49: 77–87.

Fang

Dang

. Kinematic calibration method of robots based on quantum-behaved particle swarm optimization. J Mech Eng 2016; 52: 23–30.

10.

Santolaria

Aguilar

Yagüe

, et al. Kinematic parameter estimation technique for calibration and repeatability improvement of articulated arm coordinate measuring machines. Precision Eng 2008; 32: 251–268.

11.

Liu

Zhao

Xie

, et al. Kinematic calibration and compensation for a robot with structural deformation. Robot 2015; 37(3): 376–384.

12.

H-L

Ling

S-X

, et al. Spraying robot kinematical parameter calibration based on laser tracker. Contr Eng China 2016; 31: 16–26.

13.

Kamali

Joubair

Bonev

, et al. Elasto-geometrical calibration of an industrial robot under multidirectional external loads using a laser tracker. In: 2016 IEEE international conference on robotics and automation (ICRA), Stockholm, Sweden, 16–21 May 2016. IEEE.

14.

Huangfu

Hang

Cheng

, et al. Research on robot calibration based on laser tracker. Singapore: Springer, 2016.

15.

Nubiola

Bonev

. Absolute calibration of an ABB IRB 1600 robot using a laser tracker. Robot Cim-Int Manuf 2013; 29: 236–245.

16.

Srinivas

Patnaik

. Adaptive probabilities of crossover and mutation in genetic algorithms. IEEE Trans Syst Man Cyber 2002; 24: 656–667.

17.

Correia

Moita

Moleiro

, et al. Optimization of metal–ceramic functionally graded plates using the simulated annealing algorithm. Appl Sci 2021; 11: 729.

18.

Fermani

Rossit

Toncovich

. A simulated annealing algorithm for solving a routing problem in the context of municipal solid waste collection. Commun Comput Informat Sci 2021; 1048: 63–76.

19.

Roth

Mooring

Ravani

. An overview of robot calibration. Informat Technol J 1987; 3: 377–385.

20.

Majarena

Santolaria

Samper

, et al. An overview of kinematic and calibration models using internal/external sensors or constraints to improve the behavior of spatial parallel mechanisms. Sensors (Basel, Switzerland) 2010; 10: 10256–10297.

21.

Wang

Chen

Shao

, et al. Analysis of flexible supported industrial robot on terminal accuracy. Int J Adv Robot Syst 2018; 15: 255901690.

22.

Piotrowski

Barylski

. Modelling a 6-DOF manipulator using Matlab software. Comput Sci 2014; 34: 45–55.

23.

Luo

Gao

Huang

, et al. An improved minimal error model for the robotic kinematic calibration based on the POE formula. Robotica 2021; 39: 1–20.

24.

Zhuang

Roth

Hamano

. A complete and parametrically continuous kinematic model for robot manipulators. IEEE Trans Robot Automat 2002; 8: 451–463.

25.

Gao

Liu

, et al. Kinematic calibration of industrial robots based on distance information using a hybrid identification method. Complexity 2021; 2021: 1–10.

26.

Stone

Sanderson

AC.

Statistical performance evaluation of the S-model arm signature identification technique. In: Proceedings of the 1988 IEEE international conference on robotics and automation, Philadelphia, PA, USA, 24–29 April 1988. IEEE.

27.

Park

Song

. Meta-models and genetic algorithm application to approximate optimization with discrete variables for fire resistance design of A60 class bulkhead penetration piece. Appl Sci 2021; 11: 2972.

28.

Vishnoi

Jain

Sharma

. An efficient nuclei segmentation method based on roulette wheel whale optimization and fuzzy clustering. Evolution Intellig 2019; 14: 1367–1378.