Precision Denavit–Hartenberg Parameter Calibration for Industrial Robots Using a Laser Tracker System and Intelligent Optimization Approaches

Precision object handling and manipulation require the accurate positioning of industrial robots. A common practice for performing end effector positioning is to read joint angles and use industrial robot forward kinematics (FKs). However, industrial robot FKs rely on the robot Denavit–Hartenberg (DH) parameter values, which include uncertainties. Sources of uncertainty associated with industrial robot FKs include mechanical wear, manufacturing and assembly tolerances, and robot calibration errors. It is therefore necessary to increase the accuracy of DH parameter values to reduce the impact of uncertainties on industrial robot FKs. In this paper, we use differential evolution, particle swarm optimization, an artificial bee colony, and a gravitational search algorithm to calibrate industrial robot DH parameters. A laser tracker system, Leica AT960-MR, is utilized to register accurate positional measurements. The nominal accuracy of this non-contact metrology equipment is less than 3 μm/m. Metaheuristic optimization approaches such as differential evolution, particle swarm optimization, an artificial bee colony and a gravitational search algorithm are used as optimization methods to perform the calibration using laser tracker position data. It is observed that, using the proposed approach with an artificial bee colony optimization algorithm, the accuracy of industrial robot FKs in terms of mean absolute errors of static and near-static motion over all three dimensions for the test data decreases from its measured value of 75.4 μm to 60.1 μm (a 20.3% improvement).


Introduction
Industrial robots are vital pieces of equipment for completing a modern manufacturing process. An industrial robot integrates several advanced technologies including mechanics, electronics, interfaces, and control in its structure [1]. To perform precision positioning [2], assembly, peg-in-hole [3][4][5][6][7], object handling, and manipulation on a factory floor [8], it is required to integrate accurate robot forward kinematics (FK) into its control methodology [9][10][11]. There exist some uncertainties in industrial robot FKs, which decrease its operational accuracy. Zero offset [12], irregularities in industrial robot geometry, robot manufacturing tolerances, assembly tolerances, structural deformations, and environmental factors are major sources of uncertainty in industrial robot FKs. Therefore, to increase the precision of an industrial robot FK, it is required to perform a calibration procedure [13,14]. The FK calibration is the procedure of improving the precision of the industrial robot Denavit-Hartenberg (DH) parameters to improve the robot FK model to have better positional accuracies [15,16]. The calibration of industrial robots is performed in multiple levels.
Three different calibration levels have already been identified for industrial robots [17]. The first calibration level (level I) is the one associated with robot joint angle encoders. The purpose of calibration in this level is to find a correcting relationship between the industrial robot joint angle readings from the installed joint angle transducer and the actual joint To perform the calibration of industrial robots, four optimization algorithms of DE, PSO, ABC, and GSA are utilized. The experimental setup in this paper consists of a laser tracker and collaborative robot. The laser tracker used in this study is a Leica AT960, a contactless metrology device capable of performing the 3D absolute position measurement of its retroreflector using laser interferometry principles. The positional precision of Leica AT960 is 3 µm/m (https://www.hexagonmi.com/-/media/Hexagon%20MI%20Legacy/ m1/metrology/general/brochures/Leica%20AT960%20brochure_en.ashx (accessed on 1 May 2022)), and its maximum measurement volume is up to 40 m. The absolute position measurements from the laser tracker are used to estimate the robot DH parameters. The joint angle readings from the industrial robot are used within its FK to calculate the industrial robot position. However, the positions calculated from the joint angles of the robot contain uncertainty because of the uncertainty in the industrial robot DH parameters. A cost function is defined based on the error between the FK robot positions and the measured position from a laser tracker. The data points selected for this experiment correspond to the static and near static motion of the industrial robot. Data are split to train and test to validate the accuracy of the results. The cost function is then optimized using the four optimization algorithms of DE, PSO, ABC, and GSA. DE belongs to the evolutionary optimization category of optimization, PSO and ABC belong to the swarm intelligence optimization category, and GSA is a physics-inspired optimization algorithm. The results of these optimizations give the precise DH parameters of the industrial robot. The comparison made between the results obtained using these algorithms reveals that ABC outperforms DE, PSO, and GSA in terms of the mean square value of the calibration error. It is further observed that, by using ABC, it is possible to decrease the mean absolute error of industrial robot positions for the test data from 75.4 to 60.1 µm, which is a 20.3% improvement in accuracy. The results reveal that ABC results in a superior performance compared to DE, PSO, and GSA for industrial robot DH parameter estimation purposes. This paper is organized as follows: in Section 2, the overall methodology including an industrial robot FK, intelligent optimization methods, and the proposed calibration approach are introduced. The experiment setup used for measurements is presented in Section 3. Experimental results are presented in Section 4. Section 5 concludes the paper.

Methodology
The overall calibration process, which implements metaheuristic optimization approaches to perform calibration, is presented in this section (see Figure 1). Each element in this flowchart is presented in detail in the next few sections. Although joint angle readings from an industrial robot can be used to find its end effector position, uncertainties associated with industrial robot DH parameters, obtained from its geometry, result in positional inaccuracies. The manufacturing tolerance leads to differences in DH parameter values from one robot to another of a similar build. In this paper, a laser tracker system, Leica AT960-MR, is utilized to increase the accuracy of DH parameter values.

Data Preprocessing
This laser tracker is capable of measuring accurate non-contact robotic positions at a maximum distance of 60 m. In this work, the laser tracker is used within 2.8 m from the robot base. It uses laser interferometry principles and internal precise encoders for positioning. The 3D position readings from the robot are used for DH parameter calibrations using metaheuristic optimization approaches. The position reading from the industrial robot is conducted at a sampling frequency of 125 Hz, and the sampling frequency of the laser tracker system is 10 Hz. To perform calibration, it is required to synchronize the measurements from the industrial robot joint angle encoders and resample them at the same time instances as those for the laser tracker system. It is further required to extract samples occurring at a linear velocity of less than 2 mm/s using the nominal DH parameters of UR5 to have a quasi-static calibration for the robot. Overall, we have 209 data points which are split into train and test data using a 70/30 ratio. The resulting synchronized

Data Preprocessing
This laser tracker is capable of measuring accurate non-contact robotic maximum distance of 60 m. In this work, the laser tracker is used within 2 robot base. It uses laser interferometry principles and internal precise enco tioning. The 3D position readings from the robot are used for DH paramet using metaheuristic optimization approaches. The position reading from robot is conducted at a sampling frequency of 125 Hz, and the sampling fre laser tracker system is 10 Hz. To perform calibration, it is required to sy measurements from the industrial robot joint angle encoders and resamp same time instances as those for the laser tracker system. It is further requ samples occurring at a linear velocity of less than 2 mm/s using the nomina ters of UR5 to have a quasi-static calibration for the robot. Overall, we have 2 which are split into train and test data using a 70/30 ratio. The resulting sync values are then used with a metaheuristic optimization approach to calibra rameters of industrial robots. The details of the overall calibration approach in this section.

FK Model of UR5
The mathematical function operating on industrial robot joint angles to find the Cartesian coordinates of an industrial robot within a 3D workspace is called FKs. The link transformation matrix from the link i − 1 to the link i using the DH parameters of the robot, as in Table 1, depends on the corresponding joint angle of the industrial robot and its other DH parameters [35,36] (see Figure 2).
Sensors 2023, 23, x FOR PEER REVIEW 6 of 21 where is the distance between the center of the retroreflector and the center of the robot end-effector, which is measured as equal to 0.1695 . Furthermore, to conduct the calibration, the end-effector of the robot is considered facing down, with its TCP axis-rotation vector being equal to ( 0 0).

Formulating the Estimation Problem as a Cost Function
It is required to define the industrial robot DH parameter estimation problem as a cost function to be optimized. This is a function of the industrial robot DH parameter and returns a quadratic function of the industrial robot FK error. Using each set of DH parameters given by each of the three optimization algorithms, DE, PSO, and GSA, the industrial robot FK is constructed. The positions given by the industrial robot FK are within its coordinate, which is the fixed coordinate attached to its base. The change in the coordinate The end effector positions in all three dimensions in the industrial robot coordinate are obtained as follows. Although the values of the FK parameters are unknown and will be estimated, their numerical values, according to the robot manufacturer, are as follows (https://www.universalrobots.com/articles/ur/application-installation/dh-parameters-for-calculations-of-kinema tics-and-dynamics/ (accessed on 1 May 2022)).
where d is the distance between the center of the retroreflector and the center of the robot end-effector, which is measured as equal to 0.1695m. Furthermore, to conduct the calibration, the end-effector of the robot is considered facing down, with its TCP axis-rotation vector being equal to π 0 0 .

Formulating the Estimation Problem as a Cost Function
It is required to define the industrial robot DH parameter estimation problem as a cost function to be optimized. This is a function of the industrial robot DH parameter and returns a quadratic function of the industrial robot FK error. Using each set of DH parameters given by each of the three optimization algorithms, DE, PSO, and GSA, the industrial robot FK is constructed. The positions given by the industrial robot FK are within its coordinate, which is the fixed coordinate attached to its base. The change in the coordinate is required to have the industrial robot positions in the laser tracker coordinates.
where x rl , y rl , and z rl are the robot end-effector positions using a laser tracker in a laser tracker coordinate, and Tr rl is the transformation matrix from the robot base coordinate system to the coordinate system of the laser tracker. The transformation matrix Tr rl can be easily calculated using a least squares algorithm [38]. Using the calculated transformation matrix Tr rl , the end effector positions in laser tracker coordinates are calculated as follows.
where x rl , y rl , and z rl represent the robot end effector position in the x, y, and z axes using robot joint encoders in laser tracker coordinates. The position errors (er) can be calculated for each point as the difference between the robot end effector position using joint encoder data in the laser coordinate and the robot end effector measurements from the laser tracker.
The overall mean value of position errors for all measured points associated with each of these industrial robot FKs needs to be calculated using (11). The mean absolute position errors over all measured points are used to construct the cost function for each object in the intelligent optimization algorithm.
where N m is the total number of measured points. The metaheuristic optimization procedure is then used to find the optimal values of the industrial robot DH parameters.

Metaheuristic Optimization Algorithms
Metaheuristic approaches are generally designed to solve complex optimization problems where the cost function is not explicitly given and/or suffers from multiple local minima. In this paper, four metaheuristic optimization algorithms of DE, PSO, GSA, and ABC are briefly explained, and they are used to estimate industrial robots' DH parameters.

Differential Evolution
Evolutionary algorithms are frequently used as metaheuristic optimization approaches to solve complex optimization problems. These algorithms are inspired by the natural selection principles and the Darwinian theory of evolution. Individual members of the evolutionary algorithm represent solutions for the optimization problem. In each generation of these algorithms, operators such as crossover, mutation, and selection are applied to the individuals to generate the next generation. The selection operator selects the individuals for crossover or mutations using their fitness function and some random operators. Crossover is an operator that is applied to two or more selected individuals to generate new generations as their mathematical combination. Mutation, on the other hand, is an operator that acts as a single individual to generate a new individual. The mutation operator avoids premature convergence by mutating individuals.
Candidate solutions are presented by a six-dimensional vector Individuals representing solutions to the optimization problem need to be uniformly distributed between the permeable minimum and permeable maximum value of the solutions within each dimension. The mutation operator acts on the best individual by adding a term using two other randomly selected individuals from the population as follows [32]: where V i DE is the mutant vector, X DE,best is the best individual member of DE, and the parameter F is the scale factor, which is selected as equal to 0.5 in this paper. Furthermore, t refers to the number of generations. The indices r j 1 and r i 2 are mutually exclusive integers randomly generated within the range [1, N], where N is the total number of individuals.
After performing the mutation operation, the crossover operation is applied to each individual and its mutant one [32].
where X i,j DE is the result of the crossover. The selection operator performs a comparison between the cost function associated with X i DE and U i DE , and the result of the selection procedure is used to perform the crossover and mutation to generate the next generation [32].
The mutation, crossover, and selection steps are repeated iteratively until the optimization termination condition of the maximum number of iterations is met. The maximum number of iterations considered in this paper is equal to 300.

Particle Swarm Optimization
The particle swarm optimization algorithm is based on swarm intelligence. This algorithm is developed by observing the behavior of fish and birds in a swarm [28]. The solutions to the optimization problem are presented as a position vector in a swarm. The velocity vector determines the change in the swarm position vector for the next iteration. The updates in the velocity vector are performed using the best personal experience of each swarm member and the best experience over the whole swarm. A term preserving the inertia of movement exists in the velocity vector update equation to add more exploration features to the swarm. To apply PSO to calibrate industrial robot FKs, the position vector in PSO is associated with the DH parameters as follows [28]: where X i refers to the solutions within PSO and i refers to the i-th particle within the swarm. The positions in the next generation of PSO using its current position vector and velocity vector are updated as follows [28]: where t refers to the current iteration, pbest i (t) presents the personal best experience of the i-th particle, gbest i (t) represents the overall best experience within the swarm, 0 < c 1 , c 2 are the two positive constants, and r 1 , r 2 are two uniform and random numbers from the interval of [0, 1]. The parameter c 1 is the coefficient associated with the best personal experiment of the particles in the swarm, and the parameter c 2 is the coefficient associated with the best global experiment of the particles within the swarm. The parameter w is the inertia weight, which makes the swarm follow its previous search direction. The stability criteria for PSO require the following condition to be valid for its parameters [39,40].
It is further observed in [41] that while a large value for w improves exploration, a small value guarantees a good exploitation capability for PSO.

Artificial Bee Colony
The ABC algorithm is designed to imitate the principle of foraging, which explains how bees work together to search to find the best hive and maximize the colony yield [33]. In an ABC algorithm, employed bees, onlookers, and scout bees are the three groups of bees which exchange information to perform optimization [34]. The main step of this algorithm is summarized as follows.
calculate the fitness associated with each member of the population 3.
repeat the following loop: a. produce a new set of industrial robot DH parameters as the solutions for the optimization problem using the employed bee using v ij = x ij + ϕ ij x ij − x kj , k 1, . . . , SN, j 1, . . . , 6, where ϕ ij [0, 1] is a uniform random number b.
calculate the fitness function associated with each solution f (x i ) c.
for each solution, calculate its selection probability value as follows: produce the new solutions v i for the onlookers from the solutions x i selected depending on p i and evaluate them e.
apply a greedy selection process for onlookers Other than the swarm intelligence and evolutionary optimization algorithms, there exists a third class of intelligent optimization approaches inspired by physics. GSA is an algorithm in this category inspired by the Newtonian gravitational forces between objects. In this algorithm, each solution is defined as a position of an object. The object positions are updated using velocity and acceleration vectors. Each object in this optimization algorithm benefits from a mass property. The larger mass value is assigned to the solution with a better cost function. The acceleration vector at each iteration is updated such that the objects with smaller weights are accelerated towards the ones with larger weights. This makes the objects with a smaller mass value scan the solution space towards the objects with a larger mass value. After a few iterations, all solutions converge to the object with a larger mass value, which will conclude the optimization. The mathematical description of the algorithm is explained within the next paragraphs.
Each object in this algorithm benefits from several properties of mass, positions, velocity, and accelerations. Solutions in a six-dimensional solution space are represented by object positions.
The mass value corresponding to the ith particle at iteration number t is called the non-normalized mass value and is represented by m i (t) [13].
where f worst (t) is the overall worst fitness function value and f best (t) represents the overall best fitness function value. Therefore, m i GSA (t) satisfies m i GSA (t) ∈ [0, 1], with the mass value corresponding to the best solution being equal to one and the mass value corresponding to the worst solution being equal to zero. The parameters f worst (t) and f best (t) are updated at every iteration as follows.
where x j i is the jth component of the particle position and N is the total number of particles. Each particle mass is updated and normalized at tth [13].
Sensors 2023, 23, 5368 10 of 21 The overall gravitational force (F i GSA (t)) acting on the ith particle is calculated using the gravitational law of force [13].
where k b is the number of selected best solutions, ||. || represents the Euclidean norm, ε is a small value added to prevent singulairty, r p is the power value for the Euclidean distance between the two particles, G(t) is the gravitational constant, and r j ∈ [0, 1] is a uniform and random value. The gravitational constant is updated at each iteration using the following equation [13].
where G 0 has a constant real value and t max is the maximum value of the iterations of the algorithm. The acceleration term for each object is calculated according to Newton's second law of motion by dividing the applied force to the ith mass by its mass value [13].
where A i GSA (t) ∈ R d is the d-dimensional acceleration of the particles. The velocity value corresponding to each object is updated using the acceleration term and velocity vector [13].

Hardware Setup for the Experiment
The hardware setup for performing the experiment is explained under this section. To perform the calibration test, we need the industrial robot and independent calibration equipment, as shown in Figure 3. A detailed explanation of the hardware required to perform the calibration is explained in this section.

Industrial Robot: UR5
The Universal Robots-UR5 is a six-degrees-of-freedom (DOF) industrial robot capable of handling a 5 kg load with an angular velocity of 180°/s and an angular acceleration of 180°/s . The total reach of this robot without its grippers is 850 mm. Joint angle encoders and current measurements are performed in each joint. Joint angle measurements are used as the input to the robot FK to estimate its positions and orientations. However, the nominal position repeatability of UR5 is 0.1 mm. To collect joint angle values from UR5, its main controller is connected to a PC via LAN connectivity and a hub. The

Industrial Robot: UR5
The Universal Robots-UR5 is a six-degrees-of-freedom (DOF) industrial robot capable of handling a 5 kg load with an angular velocity of 180 • /s and an angular acceleration of 180 • /s 2 . The total reach of this robot without its grippers is 850 mm. Joint angle encoders and current measurements are performed in each joint. Joint angle measurements are used as the input to the robot FK to estimate its positions and orientations. However, the nominal position repeatability of UR5 is 0.1 mm. To collect joint angle values from UR5, its main controller is connected to a PC via LAN connectivity and a hub. The PC used for this experiment operates under Linux 18.04 OS, and the software interface for the robot is provided by ROS Melodic. The ROS driver that provides the connectivity is available through a GitHub webpage (https://github.com/UniversalRobots/Universal _Robots_ROS_Driver (accessed on 9 April 2023)). This ROS driver publishes a range of rostopics including joint angle values, joint angular velocities, and motor currents over time @125 Hz. In total, 38 waypoints are programmed for the UR5, which travels them linearly in 600 s. Figure 4 illustrates the robot joint angle values over time for all six joints of UR5. The position data gathered from the robot are resampled @10 Hz to match the laser tracker frequency.
ation of 180°/s . The total reach of this robot without its grippers is 85 encoders and current measurements are performed in each joint. Join ments are used as the input to the robot FK to estimate its positions However, the nominal position repeatability of UR5 is 0.1 mm. To collect from UR5, its main controller is connected to a PC via LAN connectivi PC used for this experiment operates under Linux 18.04 OS, and the sof the robot is provided by ROS Melodic. The ROS driver that provides t available through a GitHub webpage (https://github.com/UniversalRob bots_ROS_Driver (accessed on 9 April 2023)). This ROS driver publish topics including joint angle values, joint angular velocities, and motor c @125 Hz. In total, 38 waypoints are programmed for the UR5, which tra in 600 s. Figure 4 illustrates the robot joint angle values over time for al The position data gathered from the robot are resampled @10 Hz to matc frequency.

Laser Tracker System
The laser tracker used in this experiment is AT960-MR, which is a portable dynamic 6DOF laser measurement system manufactured by Hexagon metrology GMBH, Wetzlar, Germany. This metrology equipment benefits from a single-class II laser source and is an IEC certified IP54 which guarantees ingress protection against dust and other contaminants for the unit. It benefits from Wifi connectivity, and in this work, it is set up to collect 10 data samples per second. Distance measurements are performed by using a retroreflector as the measurement target mounted on the UR5 end effector. This type of measurement system is a widely used one for precisely inspecting critical distances, locations, and surfaces [42] (see Figure 5). The target for the laser tracker is a precision Leica 1.5 red ring reflector detectable through the laser tracker at a maximum distance of 60 m. In this work, the laser tracker is used within 2.8 m from the robot base. The reflector used in this experiment is using the principle of the corner cube. To reflect the beam, three plane mirrors at right angles to one another are used. The measurement point is the center of the reflector. To perform measurements, it is required to have a clear line of sight between the laser source and the retroreflector as the laser target. the laser tracker is used within 2.8 m from the robot base. The reflector use iment is using the principle of the corner cube. To reflect the beam, three p right angles to one another are used. The measurement point is the center To perform measurements, it is required to have a clear line of sight bet source and the retroreflector as the laser target.  Table 2.  Figure 5. UR5 with a retroreflector mounted on it as the target for the laser tracker.
Further specifications and environmental conditions of the laser tracker are presented in Table 2.

Data Gathering
To perform level II calibration on the industrial robot using the laser tracker system, the data flow graph as demonstrated in Figure 1 is followed. While data are collected from the laser tracker using Wi-Fi, the LAN network is used for collecting data from UR5. The ROS melodic package under Linux 18.04 is used to collect joint data including joint angle values, joint angular velocities, and joint efforts from UR5 @~125 Hz. The laser tracker data are collected using Spatial Analyzer 2019 software under a Windows 10 operating system (see Figure 6). It is required to shift, synchronize, and resample joint angle data for UR5 to match laser tracker data samples. Furthermore, to extract a near static data sample, the angular velocities less than 2 mm/s are considered. The total number of samples satisfying this angular velocity constraint is 209 points, from which 70% is used for training and the rest is used for testing purposes. The ROS melodic package under Linux 18.04 is used to collect joint data including joint angle values, joint angular velocities, and joint efforts from UR5 @~125 Hz. The laser tracker data are collected using Spatial Analyzer software under the Windows 10 operating system (see Figure 7). The connectivity required to gather the data from the laser tracker using SA and the data from UR5 using ROS melodic are illustrated in Figure 8.
the laser tracker using Wi-Fi, the LAN network is used for collecting data from U ROS melodic package under Linux 18.04 is used to collect joint data including jo values, joint angular velocities, and joint efforts from UR5 @ ~125 Hz. The lase data are collected using Spatial Analyzer 2019 software under a Windows 10 o system (see Figure 6). It is required to shift, synchronize, and resample joint angle UR5 to match laser tracker data samples. Furthermore, to extract a near static data the angular velocities less than 2 mm/s are considered. The total number of samp fying this angular velocity constraint is 209 points, from which 70% is used for and the rest is used for testing purposes. The ROS melodic package under Linux used to collect joint data including joint angle values, joint angular velocities, a efforts from UR5 @ ~125 Hz. The laser tracker data are collected using Spatial A software under the Windows 10 operating system (see Figure 7). The connect quired to gather the data from the laser tracker using SA and the data from UR ROS melodic are illustrated in Figure 8.   values, joint angular velocities, and joint efforts from UR5 @ ~125 Hz. Th data are collected using Spatial Analyzer 2019 software under a Window system (see Figure 6). It is required to shift, synchronize, and resample join UR5 to match laser tracker data samples. Furthermore, to extract a near stat the angular velocities less than 2 mm/s are considered. The total number of fying this angular velocity constraint is 209 points, from which 70% is us and the rest is used for testing purposes. The ROS melodic package under used to collect joint data including joint angle values, joint angular veloc efforts from UR5 @ ~125 Hz. The laser tracker data are collected using Sp software under the Windows 10 operating system (see Figure 7). The c quired to gather the data from the laser tracker using SA and the data fr ROS melodic are illustrated in Figure 8.

Performance Measurement
There exist different performance indexes for validating the accuracy models. In this paper, the mean absolute error (MAE) and standard deviati performance measurement evaluation. These performance indexes are defi ically as follows.
where is the number of samples, presents the standard deviations sion, represents the data measured by the laser tracker in each of the t dimensions, and represents the measurements by the FK of the robot coordinates in each of the three positional dimensions: , , and .

Results
The parameters chosen for the three intelligent optimization algorithm GSA, and ABC are provided in Table 3. To have a better comparison, the p for all three algorithms and the number of iterations are kept the same. Fi strates the convergence graph for all these optimization algorithms versu number. It is observed from the figure that DE converges faster than the o tion algorithms of ABC, GSA, and PSO.

Performance Measurement
There exist different performance indexes for validating the accuracy performance of models. In this paper, the mean absolute error (MAE) and standard deviation are used for performance measurement evaluation. These performance indexes are defined mathematically as follows.
where N is the number of samples, σ i presents the standard deviations in each dimension, A ik represents the data measured by the laser tracker in each of the three positional dimensions, and F ik represents the measurements by the FK of the robot in laser tracker coordinates in each of the three positional dimensions: x, y, and z.

Results
The parameters chosen for the three intelligent optimization algorithms of DE, PSO, GSA, and ABC are provided in Table 3. To have a better comparison, the population size for all three algorithms and the number of iterations are kept the same. Figure 9 demonstrates the convergence graph for all these optimization algorithms versus the iteration number. It is observed from the figure that DE converges faster than the other optimization algorithms of ABC, GSA, and PSO.    Table 4 demonstrate the improvement made using the proposed calibration method. The MAE values associated with the calibrated FK using ABC in all three dimensions are reduced with respect to the uncalibrated version. The mean absolute value of error is reduced from 75.4 to 60.1 µm for the calibrated FK using ABC, which is a 20.3% improvement. It is further observed that the standard deviation of error is reduced from 100.6 to 76.1 µm. The improvement for the variance of error is 24.4%. The error trend associated with the industrial robot FK tuned by ABC as well as that of the original industrial robot FK are depicted in Figures 13-15

Conclusions
The sources of uncertainty in industrial robot FKs include ma bly tolerances, dimension measurement errors, and environmenta deals with identifying more precise robot DH parameters to incr racies of industrial robots. The proposed approach is a paramet using a laser tracker. Therefore, not only is the proposed approac positioning robots, but the parameter values obtained using this a ical meanings. To perform the calibration, a laser tracker system metrology device is used. The laser tracker used in this experime

Conclusions
The sources of uncertainty in industrial robot FKs include manufacturing and assembly tolerances, dimension measurement errors, and environmental conditions. This paper deals with identifying more precise robot DH parameters to increase the positional accuracies of industrial robots. The proposed approach is a parametric calibration approach using a laser tracker. Therefore, not only is the proposed approach a precise approach for positioning robots, but the parameter values obtained using this approach also have physical meanings. To perform the calibration, a laser tracker system which is a non-contact metrology device is used. The laser tracker used in this experiment is a Leica AT960, with precision up to 3 µm/m, performing measurements 2.8 m away from the robot base. The industrial robot used in this experiment is a UR5, an industrial robot manufactured by Universal Robots.
Using the four optimization algorithms of DE, PSO, ABC, and GSA, the tuning of industrial robot DH parameters is performed. It is observed that ABC outperforms DE, PSO, and GSA in terms of increasing the precision of robots. It is further observed that using the calibration approach proposed in this paper benefits from the ABC optimization algorithm results in decreasing the mean absolute value of error for the test data from 75.4 to 60.1 µm, which is a 20.3% improvement.

Future Works
In this paper, it is shown that the more accurate DH parameters can improve the accuracy of industrial robot FKs. The further accuracy improvement of industrial robots in future work requires online precision position feedback from measurement equipment and advanced computational inverse kinematic algorithms such as damped least squares. To increase the accuracy and convergence speed of the damped least squares algorithm, the DH parameter values obtained in this paper will be used. Furthermore, since the continuous monitoring of calibration and its validity over time is required, a mechanism for calibration monitoring is introduced to keep the position error within a valid interval. Such calibration monitoring system requires the continuous inspection and evaluation of the robot movements and positioning. This proposed approach will be implemented on more industrial robots to test the efficacy, performance, and implementability of the proposed methodology over a large number of industrial robots. The generalization of the proposed approach over a large number of industrial robots makes it possible to have statistical analysis of it and investigate the expected performance improvement using this approach in a more general form. Moreover, the calibration of industrial robots when the robot is operated at high speeds will also be considered. Contact approaches for industrial robot calibration [43,44] are usually lower-cost approaches for estimating more precise DH parameters of industrial robots. Motion sensors are another type of sensor which have already contributed to the precise calibration of industrial robots [45,46]. These types of sensors may perform orientation calibration with higher performance. The fusion of contact sensors for position measurements and inertia measurement sensors for orientation calibration would be considered in a future study.