Dynamic Model Identification for 6-DOF Industrial Robots

A complete and systematic procedure for the dynamical parameters identification of industrial robot manipulator is presented.The systemmodel of robot including joint frictionmodel is linear with respect to the dynamical parameters. Identification experiments are carried out for a 6-degree-of-freedom (DOF) ER-16 robot. Relevant data is sampled while the robot is tracking optimal trajectories that excite the system. The artificial bee colony algorithm is introduced to estimate the unknown parameters. And we validate the dynamical model according to torque prediction accuracy. All the results are presented to demonstrate the efficiency of our proposed identification algorithm and the accuracy of the identified robot model.


Introduction
In recent years, industrial robots have been greatly used as orienting devices in industry, especially in the shipbuilding, automotive, and aerospace manufacturing industries [1,2].Advanced control techniques for robots have become more and more affordable thanks to increasing power of computing resources and their dramatic cost reduction.However, the dynamical model of robot contains uncertainties in some parameters and many control methods are sensitive to their values especially in high speed operations.Hence, dynamical parameters identification approach has importance for developing model based controllers.
In terms of academic research, a standard robot identification procedure consists of dynamic modeling, excitation trajectory design, data collection, signal preprocess, parameter identification, and model validation [3].The parameter identification has attracted considerable attention from numerous researchers.Atkeson et al. [4] proposed the least square method to realize the estimation of dynamical parameters.Grotjahn et al. [5] used the two-step approach to perform the identification of robot dynamics.Gautier and Poignet [6] obtained a dynamical model of SCARA robot from experimental data with weighted least squares method.Behzad et al. [7] applied fractional subspace method to identify a robot model in simulation field.Recently, some intelligence computation algorithms have been reported as a useful tool in robot model identification.A traditional genetic algorithm (GA) was proposed to identify the autonomous underwater robot in [8].Liu et al. [9] introduced the improved genetic algorithm to obtain the space robot model.However, while dealing with complex and large-scale parameters identification problems, the GA algorithm would be stuck on local optimum.
Artificial bee colony algorithm (ABC) was first proposed by Karaboga in 2005 [10] and successfully applied to parameters identification of aerial robot [11].The ABC algorithm has been proved to possess a better performance in function optimization problems, compared with differential evolution algorithm (DE), particle swarm optimization algorithm (PSO), and GA algorithm [12].As we know, usual optimization algorithms conduct only one search operation in one iteration, but ABC algorithm can conduct both local search and global search in each iteration, and as a result the probability of finding the optimal parameters is significantly increased, which efficiently avoids local optimum to a large extent.In this paper, the ABC algorithm was introduced to conquer the parameters identification problem of the industrial robots.The identification experiment was implemented on 6-DOF ER-16 robot manipulator.

Dynamic Modeling
Since the -DOF industrial robot is represented by a kinematic chain of rigid bodies, the exhaustive description for its motion can be found in [13].The dynamic model of industrial robot is derived by the Newton-Euler or Lagrangian method: where  dyn is the -vector of actuator torques as well as the joint positions , velocities q , and accelerations q .M() is the  ×  inertia matrix, C(, q ) denotes the -vector including Coriolis and centrifugal forces, and G() is the -vector of gravity.
According to the modified Newton-Euler parameters [14] or the barycentric parameters [15], (1) can be rewritten as a linear form: where Φ dyn denotes the  × 10 observation or identification matrix, which depends only on the motion data. dyn is the barycentric parameter vector.This property considerably simplifies the parameters identification.Dynamic model of robots also contains the torques caused by joint frictions and inertias of actuator rotors apart from the effects of dynamic parameters in (2).The inertias of actuator rotors are generally provided by producers, and corresponding torques should be compensated for the dynamic equations.In fact, joint friction is a complex nonlinear model, especially during motion reversal.In order to simplify the model, the friction model consisting of only Coulomb and viscous friction [16] is given by where  fric is the friction torques and   ,  V , respectively, mean the Coulomb and viscous friction parameters.The integrated dynamic model of robots can be written as where   is actuator torques including  dyn and  fric .Φ  is the  × 12 observation matrix, and   is 12-vector of unknown dynamic parameters.In addition, the dynamic parameters of link  are governed by the form: where   ( = , , ) is the inertial tensor of link .Similarly,     denotes the first-order mass moment and   is the mass of link .
In general, the observation matrix Φ  in (4) is not a full rank; that is, not all dynamic parameters have an influence on the dynamic model.In order to obtain a set of minimum parameters, a case-by-case analysis method is adopted [17].Consequently, the dynamic model based on the basic dynamic parameters can be rewritten as where Φ is the  × ( + 2) observation matrix. is ( + 2)vector of dynamic parameters, including the basic parameters and the friction parameters. denotes the number of the minimum dynamical parameters.2 denotes the number of the friction parameters.

Basic Principles of Identification Algorithm.
In order to introduce the search mechanism of ABC algorithm, we should define three essential components: employed bees, unemployed bees, and food source [11].And the unemployed bees are divided into following bees and scout bees.The population of the colony bees is   , the number of employed bees is   , and the number of unemployed bees is   , which satisfies the relation   = 2  = 2  .We also define  as the dimension of solution vector, that is, the number of the unknown parameters.ABC algorithm treats each unknown parameter as a food source.The detailed procedure of executing the proposed algorithm is described as follows.
Step 1. Randomly initialize a set of possible solutions ( 1 , . . .,    ), and the particular solution   can be governed by where  ∈ {1, . . ., } denotes the th dimension of the solution vector.  min and   max mean the lower and upper bounds, respectively.
Step 2. Apply a specific function to calculate the fitness of the solution   according to the following equations and select the top   best solutions as the number of the employed bees: where fit  is the fitness function,   is the objective function,  is the data length, and   ( = 1, 2, 3) is the vector of the actual torques data from the first three joints.Similarly,   ( = 1, 2, 3) is the vector of the predicted data from the identified model.  ( = 1, 2, 3) is a weight coefficient between 0 and 1.
Step 3.Each employed bee searches new solution in the neighborhood of the current position vector in the th iteration as follows: where  ∈ {1, . . ., },  ̸ = , both  and  are randomly generated, and    is a random parameter in the range from −1 to 1. Next, we apply the greedy selection equation (11) to choose the better solution between V   and    into the next generation: Step 4. Each following bee selects an employed bee to trace according to the parameter of probability value.The formula of the probability method is described as Step 5.The following bee searches in the neighborhood of the selected employed bee's position to find new solutions.
Update the current solution according to their fitness.
Step 6.If the search time trial is larger than the predetermined threshold limit and the optimal value cannot be improved, then the location vector can be reinitialized randomly by scout bees according to the following equation: Step 7. Output the best solution parameters achieved at the present time, and go back to Step 3 until termination criterion  max is met.
The detailed procedure of ABC algorithm for parameters identification can be also depicted in Figure 1.

Excitation Trajectory.
When designing an identification experiment for the robot, it is necessary to design proper excitation trajectories to ensure the accuracy of estimation in presence of disturbances [18].In this work, a finite Fourier series is adopted as excitation trajectories, that is, a finite sum of harmonic sine and cosine functions.The trajectories for joint  of a robot are designed as , cos (  ) , (14) where  , is the offset term and   is the fundamental pulsation of the Fourier series.This Fourier series specifies a periodic function with period   = 2/  .Each Fourier Initialization N s , T max , limit solution and generate Trial > limit?
T > T max ?
Figure 1: Sketch of the identification algorithm.
series contains 2 + 1 parameters, and  , ,  , are the amplitudes of the sine and cosine functions.
The noise immunity and convergence rate of an identification experiment depend directly upon the constraints of the excitation trajectories.It is important to emphasize that the configurations for which measurements are taken must correspond to a well-conditioned reduced observation matrix since the constraints represent some limits for input/output.In the literature, the constraints of the excitation trajectories can be described as min cond (Φ) where  min and  max are the lower and upper of the joint positions, q max and q max are the upper of velocities and accelerations,  is optimal trajectory parameters,   is the available workspace of robot, and  max is the maximum joint torque.

Preprocessing of Measured Data.
The measured torques are obtained through collecting the data of motor current, which is described as follows: where  is motor current. is just coefficient.Since there are measurement noises in experiments, it is necessary to preprocess the collection data before identification.In order to remove outliers and attenuate the effect of interference signal, a five-spot triple smoothing method is adopted to smooth the raw data according to the following equations: where  = 3, . . .,  − 2, Y = [ 1 ,  2 , . . .,   ] is the measured raw data, and Y = [ 1 , . . .,   ]  is the data for identification after preprocessing.The more the number of using ( 17) is, the smoother the curves will be.It should be noted that excessively using (17) to smooth the raw data can lead to the error of the parameters identification increasing.
In addition, the velocities and accelerations of joints cannot be measured directly.However, these pieces of information are usually obtained by joint positions, and numerical differentiation for joint positions can amplify the measurement noise and decrease accuracy of the velocities and accelerations.An analytical approach is adopted to overcome the aforementioned difficulty, which was proposed in [19] and used successfully in [20].The average joint positions are approximated as finite Fourier series through the linear least square technique and the joint velocities and accelerations can be estimated by the derivatives of the obtained finite Fourier series.

Parameters
Identification.An experiment is conducted to test the proposed identification algorithm.The ER-16, shown in Figure 2, is a 6-DOF industrial robot manipulator without payload.And the link frame of the robot is depicted in Figure 3.The geometric parameters of the ER-16 robot are given in Table 1.Only the first joints are considered here.A fundamental pulsation of 0.05 Hz is selected for the excitation trajectories, resulting in a period of 25 s.As shown in Figure 4, the commanded trajectories are five-term    Fourier series, involving 11 optimal trajectory parameters for each joint which are listed in the Appendices, and a 0.25 Hz bandwidth.The 3D visualization of this optimized trajectory in the workspace of the robot is shown in Figure 5.The total measured time is 25 s, corresponding to 1 period of the excitation trajectory.The data is sampled with 1 kHz.Identification procedures are carried out with ABC algorithm in Matlab 2013b programming environment on an Intel Core i7-3770 PC running Windows 7. No commercial tools are used.According to [12], the performance of ABC algorithm is relative to the population size of colony bees.As the population size increases, the algorithm produces better results.However, after a sufficient value for colony size, any increment in the value does not improve the performance of ABC algorithm.And the control parameter limit is based on location vector reinitialized frequency.As the value of limit approaches infinity, the total number of location vectors reinitialized goes to zero.After many trials, in this paper, we set the parameters of ABC algorithm as follows:   = 30, limit = 15, and  max = 50.It should be noted that when the predetermined iterations exceed 50, the ABC algorithm converged and the objective value could not be obviously improved.And the search scope of the unknown parameters is also listed in the Appendices.The objective value of optimization process for the parameters is shown in Figure 6.The result shows that the convergence speed of our algorithm is fast and the final objective value which is calculated by ( 9) is 0.3182.
The robot dynamic model for the first joints contains 21 parameters, 15 base parameters, and 6 friction parameters.The parameters identified by our proposed algorithm are listed in Table 2.It should be noted that the values of the dynamic parameters of the first three joints are much bigger than those of the other three joints.It is reasonable to ignore the effect of the torques caused by the 4, 5, and 6 joints.Figure 7 compares the measured torques for the excitation trajectories with the predicted torques based on the identified dynamical model.Although the results show that the predicted error is slightly big during velocity reversal, the predicted torques have the same trend as the measured torques.It indicates that ABC algorithm has a strong ability to find the optimal parameters.
To verify the precision of the identified model by ABC algorithm, the correlation coefficient between the measured torques   and predicted torques   , defined as the normalized cross-covariance function, is applied to estimate how well the identified model can reproduce the measured torques, and the function is defined as where  = (1/)∑  =1   and   = (1/)∑  =1   .The closer the correlation coefficient is to unity, the better the identified model is.While the coefficient is close to zero, the identified model is poor.As a result, the correlation coefficients of the predicted first joints are 0.9533, 0.9856, and 0.9801, respectively.It indicates that the identified parameters have satisfactory precision.

Model Validation.
Since our aim for current stage is to investigate the validity of the model calculated by our proposed method, we focus on the validation test.Obviously, the appropriate validation test is to use the identified model in application and evaluate its success.As shown in Figure 8, three-term Fourier series are chosen as the excitation trajectories the same as the aforementioned excitation trajectories.Corresponding optimal trajectory parameters for each joint are listed in the Appendices.And the comparison of the measured and predicted torques is shown in Figure 9.It indicates that the model we obtain is capable of accurately predicting the actuator torque data.In addition, the correlation coefficients of the predicted first joints are 0.9272, 0.9534, and 0.9606 according to (18).The validation test not only shows very good results but also demonstrates that our proposed identification method is reliable enough.

Conclusion
In this paper a systematic procedure for the dynamical parameters identification of a 6-DOF industrial robot has been presented.We design optimal periodic excitation trajectories to integrate the identification experiment, data collection, and signal preprocess.All the unknown parameters are well identified by ABC algorithm.When comparing the measured torques and the predicting torques, we conclude that our proposed method can accurately estimate the robot dynamical parameters.Further, a model validation has been carried out to verify the validity of the identified model.The results of this paper are useful for researchers and manufactures of industrial robots.

A. Optimal Parameters of Excitation Trajectories
The optimal trajectory parameters for the five-term Fourier series are listed as follows: The optimal trajectory parameters for the three-term Fourier series are listed as follows:

B. Search Scope of Dynamical Parameters
After many trials, the probable suboptimal or optimal search scope of the unknown parameters is listed in Table 3.

Figure 7 :
Figure 7: Comparison of the measured torques and predicted torques.