A Robotic Calibration Method Using a Model-Based Identification Technique and an Invasive Weed Optimization Neural Network Compensator

The study proposed a robotic calibration algorithm for improving the robot manipulator position precision. At first, the kinematic parameters as well as the compliance parameters of the robot can be identified together to improve its accuracy using the joint deflection model and the conventional kinematic model calibration technique. Then, an artificial neural network is constructed for further compensating the unmodeled errors. The invasive weed optimization is used to determine the parameters of the neural network. To show the advantages of the suggested technique, an HH800 robot is employed for the experimental study of the proposed algorithm. The improved position precision of the robot after the experiment firmly proves the practicability and positional precision of the proposed method over the other algorithms in comparison.


Introduction
Robot manipulators are widely used in industry to attain many duties such as welding, painting, pick and place task, etc. The construction of robot manipulators is characterized using their kinematic model parameters. However, experience has shown that industrial manipulators have much greater repeatability than accuracy [1]. Therefore, the utility of robot manipulators would be significantly enhanced if they were made to be as accurate as they are repeatable. In producing and assembly, many errors arise that could not be taken into account by the nominal geometric model. For that reason, there is a demand to create model-based robotic calibrations that depend on an error model that symbolizes the connection between the errors of geometric parameters and the end effector positioning errors. The kinematic calibration method has been widely researched by numerous studies [2][3][4][5][6][7][8][9][10][11][12][13] to describe the proper geometric model. Denavit, Hartenberg et al. proposed one of the most fundamental calibration methods that is widely used. The DH model is based on homogenous transformation matrices. The model is a description of the kinematic relations between the links of a kinematic chain connected by 1 degree-of-freedom lower pair joints [6][7][8]. There are also other methods such as the zero-reference position model that is proposed by Gupta et al., which has been widely used by many studies [9,14,15]. Another fundamental model-based calibration method is the complete and parametrically continuous (CPC) model, which relies on the singularity-free line representation by Zhuang et al. [10,11]. Park and Okamura [12] employed the product of exponentials (POE) error model to robot calibration. This model has been adopted by several researchers [4,13]. There have been some attempts to replace the least square estimation with the Kalman filter, particle filter, fuzzy theory, etc. [16][17][18]. However, those attempts did not seem to be very effective.
In addition, the kinematic errors are not the only source of robot manipulator errors. Whitney et al. [19] used a PUMA 560 robot to determine that the most significant error sources for the robot were several nongeometric errors. Judd and Kasinski [20] studied an AID-900 robot and found that the geometric errors were responsible for approximately 95% of the measured error. Jean et al. [21] and Becque [22] reported that flexibility in joints and in links is responsible for 8-10% of the total position and orientation errors. Therefore, more recently, researchers have been devoted to the identification of compliance errors as well as kinematic errors in industrial robots [5,20,23]. Jian Zhou and Hee-Jun Kang [5] proposed a method for simultaneously identifying the kinematic and compliance parameter of the robot. The model-based calibration has been widely applied due to its fast computing and knowledge of error sources [5,24,25]. However, the accuracy of the model-based calibration method is dependent on the model accuracy. It is practically impossible to create a model that contains all the sources committing to the end effector errors. Some error sources such as gear backlash, temperature variation, and other errors are difficult to model correctly and completely.
To archive further accuracy, nongeometric calibration has been widely studied for compensating for the sources of errors that could not be taken into account by geometric calibration [17,21,[26][27][28][29][30][31][32][33]. The conventional back-propagation neural network (BPNN) [34] is widely adopted by researchers [35] for compensating the unmodeled errors to increasing the precision of the robot. In robotic calibration processing, BPNN is usually employed to construct the relationship between the end effector position and the corresponding joint angle configuration [36][37][38][39]. However, the conventional BPNN has some drawbacks such as getting stuck in local minima and slowing convergence [40]. To overcome these drawbacks, some heuristic algorithms have been used for training the network [41][42][43][44][45]. One among them is the Invasive Weed Optimization (IWO) algorithm.
In 2016, Mehrabian et al. presented the Invasive Weed Optimization (IWO) algorithm [46]. This method is motivated by the spreading of weeds to find a suitable place for expanding and breeding. The technique is characterized by fast convergence for global optimization. Based on the properties of IWO, this algorithm is employed for optimizing the parameters of a neural network to reduce its cons such as a high dependence on input data. Furthermore, the IWO prevents the neural network from falling to local minima.
This work proposed a new calibration method that includes the model-based calibration technique and unmodeled calibration method. At first, the kinematic and compliance parameters are simultaneously determined by a model-based calibration method [5]. Then, the method needs to compensate for some unmodeled errors that cannot be ignored such as friction, mechanical transmission error, and thermal expansion. A suggested neural network optimizing by IWO is used to compensate the residual positional errors. It should be noted here that the proposed method is a combination of model-based and artificial neural network (ANN) methods that used the IWO technique to determine the weight and bias. Meanwhile, most of the ANN-based technique is applied after the kinematic calibration. This calibration method simultaneously calibrates both the kinematic errors and joint compliances. After the simultaneous calibration, the IWO-NN-based compensation is accomplished for the unmodeled nongeometric errors. By using the IWO neural network, the proposed calibration method seems to reach the global minima easily. Therefore, the IWO neural network can be said to have better convergence capability than the traditional backpropagation NN. Finally, a HH800 robot is employed for the experimental study of the proposed algorithm to compare with four other calibration methods, including the conventional kinematic calibration method (KM), as well as the simultaneous identification of joint compliance and kinematic parameters method (SKCM), and the combination of an NN compensator and SKCM method (NN-SKCM). The advantages of the method have been shown: the enhanced position accuracy of the manipulator after the calibration confirms the feasibility and greater positional accuracy over the other calibration methods. Additionally, the adopted IWO neural network has better convergence capability than the back-propagation neural network in this calibration process. This advantage makes the proposed method more feasible in real offline programming environments.

Kinematic Model of the HH800 Robot
The HH800 robot [36] model is given in Table 1 and Figure 1. The homogeneous transformation matrix from the robot base frame to the end effector frame can be computed by:  The HH800 robot [36] model is given in Table 1 and Figure 1. The homogeneous transformation atrix from the robot base frame to the end effector frame can be computed by:

D-H Parameters of the Main Open Chain
The homogenous transformation matrix from the robot 6th frame to the end effector frame: The passive joint position 3 is formed from the joints 2 and 3 as follows [36,47]:  The homogenous transformation matrix from the robot 6th frame to the end effector frame:

D-H Parameters of the Main Open Chain
The passive joint position θ 3p is formed from the joints θ 2 and θ 3 as follows [36,47]:

Simultaneous Joint Stiffness and Kinematic Parameters
In a robot static configuration, a robot joint torque causes a twist deformation about a rotation shaft (that represents the entire drive train from the motor to the associated robot link). Therefore, the shaft can be considered as a torsional spring in the compliance modeling. This study investigates the characteristics of a torsion spring because they are related to modeling of rotational joint compliance.
The characteristics of torsion springs are basically presented by non-linear functions: for example, τ = k 1 * ∆θ c + k 2 * (∆τ c ) 3 , where τ is the spring torque, ∆θ c is the spring rotational deformation, and k 1 and k 2 are the coefficients. When the robot joint deformation is small, the linear part becomes dominant. Now, we can assume that the functional relationship between the joint torque and its deformation is linear in this calibration process.
Assuming that the compliance of the robot links is small and can be neglected, therefore, the elastic errors are mostly dependent on the compliance of joints under the link self-gravity and external payload. Accepting that the robot joint deformations can be expressed by a linear function of the joints torque, the deformation of the ith joint can be expressed by the effective torques τ i : where k i is the joint stiffness value of the ith joint. The Cartesian position errors due to small joint deflections can be modeled as: The compliance vector is denoted by C = c 1 c 2 · · · c n T . The deflection vector is At the equilibrium position, τ = diag(τ 1 , τ 2 , · · · , τ n ) is the effective torque of robot joints. The real position vector of the robot end-effector can be expressed as: where P kin is the end effector position that is calculated by the kinematic parameters. ∆P kin , ∆P c , and ∆P extra are the position errors due to the kinematic parameter errors, joint deflections, and the residual errors due to the unmodeled sources, respectively. Jian Zhou and Hee-Jun Kang [5] presented a method for simultaneous identifying the kinematic and compliance parameter of a robot: where ∆X is a (3 × 1) vector of three position errors of the robot end-effector. J is a (3 × p) matrix that relates the column vectors ∆X and ∆P kin (p = 27 is the total number of kinematic parameters). By using the least-square method, the kinematic and joint compliance parameter can be computed at the same time.

IWO-NN Errors Compensator Technique
In Equation (6), the position errors ∆P extra are still large after applying the algorithm of simultaneous identification of joint compliance and kinematic parameters. This occurred because of some unmodeled sources such as friction, thermal extension, mechanical transmitting errors, which are hardly considered in geometric calibration. To decrease these unmodeled errors, a nongeometric compensator should be carried out for compensating. In this study, an IWO-NN errors compensator is taken to eliminate the residual positional errors. The method is fully described as follows.
First, the robot geometry error ∆φ and joint compliance parameter C are identified together (Equation (7)) using the least squares method. The total positional error vector is denoted by ∆P (3 × 1): where P m , P c indicate the positions of the end effector by measuring and computing. ∆P kin , ∆P c are the kinematic error and compliance error that are calculated by applying the method in Ref. [5].
The residual position error after the modeled-based calibration process is ∆P r . Here, the IWO-NN is hired to compensate the residual error ∆P r .
The IWO-NN containing six inputs representing the joint configurations is θ n = [θ 1 , θ 2 , · · · , θ 6 ]. There are five nodes in the hidden layer ( Figure 2). The tag-sigmoid is selected as the activate function for the hidden layer as follows [36]: Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 13 Three outputs of the NN have the linear activation function. The outputs of the NN are used for compensating the residual errors ∆ . The error in the output layer is where is the output of the neural network. The mean square error is calculated as follows: where m = 3 represents the three dimensions of the end-effector position. The IWO is a recently proposed population-based heuristic optimization algorithm that mimics the spreading of weeds to find a suitable place for expanding and breeding [46]. IWO has four main stages: (i) initialization, (ii) reproduction, (iii) space distribution, and (iv) ranking and selection.
(i) Initialization: The population of solution is generated. The seeds are randomly distributed over all the space problem. The number of seeds is chosen.
(ii) Reproduction: In this step, every seed reproduces the next generation. The number of every seed's heir depends on its fitness in the community. If the seed has a good fitness, it will produce more seeds than another seed that has a worst fitness in the society.
(iii) Space Distribution: This is the progress by which the offspring are randomly spread over all the searching space. The spreading process is made to distribute the offspring close to their parents' location. Therefore, the seeds with lower fitness in the population are deleted over time. The distribution function is described as below: where is the standard deviation (SD) at the present step, is the maximum number of The error in the output layer is e = ∆P r − P nn (11) where P nn is the output of the neural network. The mean square error is calculated as follows: where m = 3 represents the three dimensions of the end-effector position. The IWO is a recently proposed population-based heuristic optimization algorithm that mimics the spreading of weeds to find a suitable place for expanding and breeding [46]. IWO has four main stages: (i) initialization, (ii) reproduction, (iii) space distribution, and (iv) ranking and selection.
(i) Initialization: The population of solution is generated. The seeds are randomly distributed over all the space problem. The number of seeds is chosen.
(ii) Reproduction: In this step, every seed reproduces the next generation. The number of every seed's heir depends on its fitness in the community. If the seed has a good fitness, it will produce more seeds than another seed that has a worst fitness in the society.
(iii) Space Distribution: This is the progress by which the offspring are randomly spread over all the searching space. The spreading process is made to distribute the offspring close to their parents' location. Therefore, the seeds with lower fitness in the population are deleted over time. The distribution function is described as below: where σ iter is the standard deviation (SD) at the present step, iter max is the maximum number of iterations before stopping the algorithm, σ initial is the initial value of the SD, σ f inal is the final value of the SD, and n is the non-linear modulation index.
(iv) Ranking and Selection: When the maximum number of populations in the colony is hit, all the seeds are evaluated, including new seeds and their producers. The unqualified fitness seeds will be eliminated.
In this work, the IWO is introduced for optimizing the weights and biases of an NN. The optimizing process is described by the flowchart in Figure 3.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 13 member of the population is = 53 according to the number of the weights and biases of the fully connected NN, which has six inputs, three outputs, and five hidden nodes.
The loss function is used for evaluating the fitness of member of the population: Here, is the matrix of mean square errors (Equation (12)) of the ith set of NN respectively to 50 inputs and outputs of the calibration data. 1 = 1; 2 = 0.7 are the weights that represent the contribution of mean and max function in forming the cost function of each solution i. By using Equation (14), the number of seeds for member of the population can be calculated: where is the function that rounds to the nearest integer less than or equal to the input number. is the ratio calculated by the following equation: The seed of member can be calculated by = + * (17) Assume that X i = {x i1 , x i2 , · · · , x in } is the ith member that contains the weights and biases of the NN. An initial population including 10 members is randomly spread over the searching space. The maximum population size is 50. The minimum number of seeds is S min = 0, and the maximum number of seeds is S max = 5. The variance reduction exponent is n = 2 (Equation (13)). The initial value of the SD is σ initial = 0.05, while the final value of the SD is σ f inal = 0.0005. The dimension of member X i of the population is i = 53 according to the number of the weights and biases of the fully connected NN, which has six inputs, three outputs, and five hidden nodes.
The loss function is used for evaluating the fitness of member X i of the population: Here, E i is the matrix of mean square errors (Equation (12)) of the ith set of NN respectively to 50 inputs and outputs of the calibration data. l 1 = 1; l 2 = 0.7 are the weights that represent the contribution of mean and max function in forming the cost function of each solution i.
By using Equation (14), the number of seeds for member X i of the population can be calculated: where f loor is the function that rounds to the nearest integer less than or equal to the input number. r i is the ratio calculated by the following equation: The seed S j of member X i can be calculated by using σ iter (Equation (15)). K is a matrix that has the same dimension as X i . Every element of K is a random number between 0 and 1. For speeding the process, the values of elements of S j are bounded. The lower and upper bounds are −1.7 and 1.7, respectively. After this process, all the seed S j of member X i are merged into the population and considered members of the population. If the number of members overcomes the maximum population size, the worst members having the highest valued loss function are eliminated. Then, the processing is repeated again until it reaches the maximum iteration or one of the members archives the desired loss function value. Figure 4 shows the flowchart of the proposed method.

Experiment and Validation Results
An HH800 robot is employed to examine the proposed method to clarify the effectiveness and practicability of it. The proposed method is examined in contrast with the other four methods to show its advantages. To demonstrate the effectiveness of the proposed method (IWO-NN-SKCM) in

Experiment and Validation Results
An HH800 robot is employed to examine the proposed method to clarify the effectiveness and practicability of it. The proposed method is examined in contrast with the other four methods to show its advantages. To demonstrate the effectiveness of the proposed method (IWO-NN-SKCM) in improving the robot position accuracy, four calibration methods are carried and compared in both the calibration and validation process. The conventional kinematic calibration method (KM) [6], as well as the simultaneous identification of joint compliance and kinematic parameters method (SKCM) [5], and the combination of NN compensator and SKCM method (NN-SKCM) [36] are used in this experiment calibration.

Experimental Calibration Results
In this experiment, a robot calibration system consists of a Hyundai HH800 robot (6 degree of freedom (d.o.f)) that has one closed-loop actuating mechanism, a 3D point sensing device (API Laser Tracker, measurement accuracy of 0.01 mm/m, repeatability of +/−0.006 mm/m), and an accompanying laser reflector. The reflector is fixed at a particular location of the robot end-effector. The system is arranged as shown in Figure 5.

Experiment and Validation Results
An HH800 robot is employed to examine the proposed method to clarify the effectiveness and practicability of it. The proposed method is examined in contrast with the other four methods to show its advantages. To demonstrate the effectiveness of the proposed method (IWO-NN-SKCM) in improving the robot position accuracy, four calibration methods are carried and compared in both the calibration and validation process. The conventional kinematic calibration method (KM) [6], as well as the simultaneous identification of joint compliance and kinematic parameters method (SKCM) [5], and the combination of NN compensator and SKCM method (NN-SKCM) [36] are used in this experiment calibration.

Experimental Calibration Results
In this experiment, a robot calibration system consists of a Hyundai HH800 robot (6 degree of freedom (d.o.f)) that has one closed-loop actuating mechanism, a 3D point sensing device (API Laser Tracker, measurement accuracy of 0.01 mm/m, repeatability of +/−0.006 mm/m), and an accompanying laser reflector. The reflector is fixed at a particular location of the robot end-effector. The system is arranged as shown in Figure 5.

Laser Tracker
Reflector Location HH800 Robot  In order to acquire suitable measurement data for robot parameter identification and IWO-NN training, the robot moves its end-effector to positions such that they entirely cover the workspace. The three-dimensional coordinates of the end points are measured by the Laser Tracker and saved in a computer. At the same time, the associated robot joint readings also are recorded.
These measurements are grouped as follows: A set of 40 robot configurations (Q1) is used for parameter identification and collected over all the workspaces. By using the SKCM method, the kinematic and joint compliance parameters are identified. These parameters are presented in Tables 2 and 3, including four joint compliance parameters and 29 kinematic parameters. Another set (Q2) of 50 robot configurations is also randomly selected to determine the weights and biases of the neural network that has five hidden nodes, six inputs, and three outputs. The reason why we use a different set (Q2) from the first set (Q1) is that the neural network compensator has more general error compensation capability over the entire robot workspace. A set of 50 arbitrary endpoints (Q3) for robot accuracy validation is collected over all the workspaces.  Table 3. Identified D-H parameters of Hyundai robot HH800 ("-": unavailable, "X": unidentifiable). In this calibration progress, the HH800 robot is used to apply four different calibration algorithms. Table 4 and Figure 6 demonstrate the results of these calibration methods. Figure 6 provides a visual result of the absolute position errors of each calibration pose using four calibration methods, including the conventional kinematic calibration method (KM), as well as the simultaneous identification of joint compliance and kinematic parameters method (SKCM), and the combination of NN compensator and SKCM method (NN-SKCM). It shows that the position errors generated by the proposed method are the lowest and better converging in comparison with the other methods. By using the IWO-NN SKCM, the mean of errors is reduced from 0.5961 to 0.3450 mm by as much as 42.12% compared to the KM method. This mean of errors is also decreased as much as 31.52% in comparison to the SKCM method (0.5038 mm). In comparison to the NN-SKCM method, the mean of errors result is declined 15.92% (0.4103 mm). The IWO-NN is better than NN-SKCM in reducing the maximum of absolute position errors (41.88%). The proposed method also generates the lowest maximum position error (0.6374 mm) and standard deviation (0.1624 mm).

Experimental Validation Results
The validation process is examined to illustrate the general capability of the suggested calibration technique over the entire robot workspace. In this progress, the set Q3 includes 50 robot configurations of the manipulator that are not used in the calibration process, and it is taken to examine the four different methods. The results of the validation process are demonstrated in Figure 7 and Table 5. Figure 7 shows the residual error of 50 poses using four methods in the validation process. It is clear to see from the figure that the proposed method is the best over both the positions used in the calibration process and the general position in the overall workspace.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 13 Figure 6. Absolute position error of the HH800 robot after calibration.

Experimental Validation Results
The validation process is examined to illustrate the general capability of the suggested calibration technique over the entire robot workspace. In this progress, the set Q3 includes 50 robot configurations of the manipulator that are not used in the calibration process, and it is taken to examine the four different methods. The results of the validation process are demonstrated in Figure  7 and Table 5. Figure 7 shows the residual error of 50 poses using four methods in the validation process. It is clear to see from the figure that the proposed method is the best over both the positions used in the calibration process and the general position in the overall workspace.
By using the proposed method, the mean of errors is reduced from 0.5981 to 0.4662 mm by as much as 22.07% better in comparison to the KM method. The mean of errors is also decreased by as much as 26.72% in comparison to the SKCM method (0.5458 mm), and there is a 10.14% decline from the mean of errors by using the NN-SKCM method (0.5187 mm). The proposed method also generates the lowest maximum position error (0.6538 mm) and standard deviation (0.1333 mm).

Advantages of the IWO-NN Compensator
In this work, the robot configurations are divided into several datasets that are compensated by  By using the proposed method, the mean of errors is reduced from 0.5981 to 0.4662 mm by as much as 22.07% better in comparison to the KM method. The mean of errors is also decreased by as much as 26.72% in comparison to the SKCM method (0.5458 mm), and there is a 10.14% decline from the mean of errors by using the NN-SKCM method (0.5187 mm). The proposed method also generates the lowest maximum position error (0.6538 mm) and standard deviation (0.1333 mm).

Advantages of the IWO-NN Compensator
In this work, the robot configurations are divided into several datasets that are compensated by the NN and the proposed IWO-NN. From this experimental application, the NN appear to be usually fallen into the local minima, and a reinstalling of random weights and bias is required to reach the global minima. On the other hand, the IWO-NN seems to archive the global minima quite efficiently. For that reason, the proposed IWO-NN can be said to have better convergence capability. Moreover, an accuracy consistency of experimental calibration (0.345 mm) and validation (0.4662 mm) confirms the abilities of the suggested calibration method.

Conclusions
In this study, a new calibration method with an error compensating IWO neural network is proposed for enhancing the robot positional accuracy of the industrial manipulators. By combining the joint deflection model with the conventional kinematic model of a manipulator, the geometric errors and joint deflection errors can be simultaneously considered to increase its positional accuracy. Then, a neural network is designed to additionally compensate the unmodeled errors, such as friction, mechanical transmission error, thermal expansion specially, and nongeometric errors. The teaching-learning-based optimization (IWO) method is employed to optimize the weights and biases of the neural network.
Real experimental studies are carried out on the HH800 manipulator to show the efficiency of the proposed method. The advantages of the method had been shown, such as the enhanced position accuracy of the manipulator after the calibration confirms the feasibility and its greater positional accuracy over the other calibration methods. Additionally, the adopted IWO neural network has better convergence capability than the back-propagation neural network in this calibration process. From this implementation experience, while the back-propagation neural network seems to be easily got into the local minima and to need reiteration by randomly resetting the weights and biases to reach the global minima, the IWO-NN seems to quite easily reach the global minima. This advantage allows that the proposed method is more feasible in a real offline programming environment. However, a heuristic optimization method is used to determine the weights and bias of the NN. Therefore, the proposed method takes time for computing.
In the future, the work could be expanded by modeling the relationship of the joint deflections of the robot and the effective torques as a polynomial function to increase the precision of modeling the robot's joints compliance.