Nonlinear Friction and Dynamical Identification for a Robot Manipulator with Improved Cuckoo Search Algorithm

,


Introduction
Modern industry is increasingly oriented towards the production of small batches of a large variety of products, asking for flexibility and automation in the manufacturing systems [1,2].Additionally, the increasing quality standards, international competition, and economic reasons put higher requirements on reliability and accuracy and especially on the speed of production processes.In this context, industrial serial robot manipulators have become an indispensable means of automation to increase productivity and flexibility of production units.Robots are programmed by teaching the sequence of the attitude and position which are necessary to execute the desired task.To reach a sufficient accuracy, this teaching is mostly done on-site and relies on a good repeatability, rather than on a good absolute accuracy.To improve the operation accuracy, a precise dynamical model is essential for accurate offline programming.
In the academic aspect, a typical manipulator identification procedure consists of dynamic modelling, excitation trajectory design, data collection, signal preprocess, parameter identification, and model validation [3].When a priori knowledge is available about the robot system, parametric models can be derived based on the laws of physics and mechanics resulting in a set of differential equations.The unknown dynamical parameters have a physical meaning and can be identified by several approaches.Atkeson et al. [4] used the least square method (LS) to implement the load estimation of dynamical parameters on a PUMA600 robot.According to weighted least squares method (WLS), Gautier and Poignet [5] proposed a dynamical identification approach only from the torque data, without other sensors.Grotjahn et al. [6] used the two-step method to execute the friction and rigid body identification of robot dynamics.Considering the effect of measurement noise, Olsen and Petersen [7] used the maximum likelihood estimating (MLE) method for parameters identification of an industrial robot with a statistical framework.Recently, some novel dynamical identification methods for robot manipulators have been reported using intelligence algorithms.For instance, Bingül and Karahan [8] integrated the particle swarm optimization (PSO) algorithm with LS algorithm to estimate the dynamical parameters of Staubli RX-60 robots.In the identification experiment, the velocity and acceleration are measured by three high-speed cameras, and the joint torques are measured by six load cell sensors.Without velocity and acceleration sensors, Ding et al. [9] used the motor current and joint positions to calculate the torques, velocity, and acceleration of joints and proposed an artificial bee colony algorithm (ABC) to obtain the unknown dynamical parameters.Nevertheless, when handling complicated and high-dimensional parameters identification problems, the flaw that the premature convergence can make those intelligence computation algorithms stuck in a local optimum.
As known, the friction is a major source of disturbances affecting motion quality.Therefore, it must be included as an additional component in robot modelling.In robot identification applications, a model including Coulomb and viscous friction is frequently applied [10,11].With such linear model, the parameters estimation is significantly simplified.However, this friction model is not capable of describing the experimentally measured friction characteristic, especially the static model at joint reversal.Aiming at these problems, we add a nonlinear friction model into the dynamical model of a 6-DOF industrial serial manipulator and use an improved cuckoo search (ICS) algorithm for dynamical identification.The idea of the method is to measure the positions and gravitational torques of different joints through designing Fourier series as excitation trajectories.The collected values are used to calculate the dynamical parameters based on ICS.In ICS, considering the outstanding performance of chaotic operator and emotion operator, these improved operators are used to enhance the performance of the classical CS.And the comparison of three different identification methods, CS, PSO, and ICS, illustrated the superiority of our proposed algorithm in the application of dynamical identification.
This article is organized as follows.The dynamical model with a nonlinear friction model of a robot manipulator is given in Section 2. For the unknown dynamical parameters, the improved cuckoo search algorithm is introduced to realize the parameters estimation in Section 3.Then, the design of excitation trajectory, data collection, and preprocessing are presented in Section 4. Later, an ER-16 robot is used as a test platform for identification experiment, and the results are analyzed in Section 5.In addition, the superiority of linear and nonlinear friction model has been compared through the model validation in Section 6.Finally, Section 7 discussed the key findings and prospective research target.

Dynamic Modelling
According to the literature [12], a -DOF serial manipulator is described as a kinematic chain of several rigid bodies.Hence, we can utilize the Newton-Euler method to deduce the dynamical model of the manipulator: where joint torque , joint position q, joint velocity q , and joint acceleration q are -dimensional vector.  denotes the -dimensional joint friction vector.M(q) represents  ×  inertial matrix, C(q, q ) is a -dimensional vector including Coriolis and centrifugal forces, and G(q) is -dimensional gravity vector.
Equation (1) except for friction torques can be rewritten as a linear form if using the modified Newton-Euler parameters [13] or the barycentric parameters [11]: where Φ  only contains the motion data, which can be treated as ×10 identification matrix or observation matrix.  is the barycentric parameter vector.This conversion vastly reduces the complexity of parameters identification.In addition, the dynamical parameters of link  are governed by the form where   ( = , , ) is the inertial tensor of the link .
Similarly,   is the mass of the link  and     is the inertia moment.
Generally, the identification matrix Φ  in (2) is not full of rank, that is, not all dynamical parameters give contribution to the joint torques.In the literature [14], some methods like case-by-case analysis or singular value decomposition are adopted to eliminate the redundant parameters.And the barycentric parameter vector   can be replaced by a vector of minimal barycentric parameters  ∈ R  with  < 13.Hence, (2) can be transformed into another form  = Φ (q, q , q ) , where Φ denotes the × observation or identification matrix and  is the number of minimal barycentric parameters.Except for the dynamical parameters in (2), there also exists friction torques and extra torques caused by inertias of motor rotors.In general, the inertias of motor rotors are provided by manufacturers, and corresponding torques can easily be compensated to the dynamical model.As for joint friction model, it is regarded as a complex nonlinear model.To simplify the model, the Coulomb and viscous friction were used to describe the friction model.But the researchers found [15] that the friction torques of some joints exceeded their full speed range, and the simple friction model could not cover the characteristics, especially at motion reversal.A better description of the joint friction characteristics may be based on the following nonlinear equations: where  0 is the zero drift error of friction torque,   is the Coulomb friction coefficient,  V is the viscous friction coefficient, and   and   are the experiential friction coefficients.
It should be noted that this model has a discontinuity at zero velocity.In summary, using the nonlinear friction model yields the whole dynamical model of the 6-DOF robot as where   denotes the 5 friction torques vector and 5 is the number of friction coefficients.Obviously, the classical least square method could not solve the above nonlinear equation.Hence, applying an intelligence algorithm for solving this problem may be a feasible method.

Introduction to Chaos
Theory.Chaos theory is epitomised by the so-called "buttery" detailed by Lorenz [16].He discovered that tiny changes in an initial state would make a radically different final result and typically rendering long-term prediction impossible.Chaotic map has ergodic and stochastic properties, which is regarded as a bounded nonlinear system with deterministic dynamical behavior.Moreover, it has a very sensitive dependence on initial conditions.There are many forms about chaotic map, such as tent map, Gauss map, logistic map, and tent map [17].Considering the high robustness and stability of tent map, we choose it to generate the chaotic sequence.The description of the tent map is written as where  is iterations and  denotes a positive number.The value of   is updated with the initial condition   ∈ (0, 1).A tent map for  = 0.65 after 500 iterations is shown in Figure 1.
As the value of  increases,  +1 gets a new value.In this paper, the parameter  is chosen as 0.6 after many trials.The local search can be described as the following formulation, which gains an insight into CS search mechanisms:

Basic
where    and    are two different random solutions, (⋅) represents a Heaviside function, the symbol "⊗" represents point-to-point multiplication,  is a random parameter which satisfied a uniform distribution,  is step size, and  is a scaling coefficient over zero.
For the global search, the process is controlled by a Lévy flight behavior.The mathematical formula for this behavior is described as follows: where (, ) is a Lévy flight function that complicates integration,  represents the power coefficient, and The basic steps of CS are presented as follows.
Step 1 (parameters initialization).The initial positions of nests are randomly generated from where  ∈ {1, . . ., },  ∈ {1, . . ., },  is the position size,  is the dimensionality of unknown dynamical parameters, and   and   are the lower and upper limits in the solution space.
Step 2 (objective value calculated).In this paper, CS searched for the dynamical identification by minimizing the objective function written as follows: where   is the objective function,  is the data length, and   ( = 1, 2, 3) denotes joint torques measured from the first three joints.Similarly,   ( = 1, 2, 3) denotes predicted torques calculated by the identified model.  ( = 1, 2, 3) is a weight coefficient between 0 and 1.
After generating a group of new nest positions, the superior positions are retained to the next generation by comparing the objective value between the original positions and the postupdate nest positions.
Step 4. A random number of uniform distributions  ∈ [0, 1] is compared with   .If  >   , then the located nest is changed randomly.Otherwise, the nest remains unchanged.
Step 5.The termination condition of the CS is determined.If the condition is satisfied, then the optimal solution becomes an output.Otherwise, Step 3 is repeated.

Improved Operator.
In general, the parameter  in Lévy flight is the key factor to affect the convergence of CS.Due to the infinite variance and mean, the classical CS algorithm may have a premature search process.To overcome this problem, the tent map is used to generate a chaotic sequence for parameter .In that case, the algorithm searches the new position in the neighborhood of the current optimal position.Meanwhile, a new emotional acceptance criterion is used to prevent the algorithm from getting trapped into local optima.
In the ICS, two cuckoos' emotions (positive and negative) can be described as follows: where  represents the function of cuckoo's emotion,  is a constant,  is the stimulus function,  0 is a stimulus threshold, and  is the objective function.It should be noted that  is selected as 1 and  =   .

Design of Exciting Trajectories and Data Preprocess
It is essential to consider these conditions before designing an identification experiment, that is, (1) whether the excitation trajectory is sufficient to provide fast and accurate parameters estimation and (2) whether the processing of the experimental data is simple and yields stable and accurate results.In fact, the imprecise modelling and measurement noise are assumed to be an additive frequently distributed zero-mean stochastic disturbance, causing bias errors and uncertainty in the dynamical identification.To reduce the effect of the disturbance, an appropriate excitation trajectory must be designed carefully.The choice of parameterization for the excitation trajectory is a very important issue.It directly determines the number of parameters in the optimization problem and the effort needed to calculate velocity and acceleration from the joints positions measurements.The first considerations on finding excitation trajectory for the dynamical identification of manipulator were proposed by Armstrong [19].Then, Grotjahn and Daemi [20] proposed an interpolated trajectory consisting of two parts, part I overcame the given boundary conditions and part II overcame the homogeneous boundary conditions.Furthermore, Gautier [21] used fifth-order polynomials to obtain smooth joint trajectories to be executed by the manipulator, and the polynomial coefficients are fixed by imposing continuity constraints between the trajectory segments.In this paper, we adopt a finite Fourier series which was proposed by Swevers et al. [22] as excitation trajectories.The trajectory for joint  of a manipulator is designed as , cos (  ) , (13) where  ,0 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 series contains 2 + 1 parameters, and  , ,  , are the amplitudes of the trigonometric functions.It should be noted that although this trajectory is hard to implement on an industrial setup, this parameterization has several advantages.It is possible to calculate the velocity and acceleration in the frequency domain, which avoids phase distortions.
The antinoise ability and convergence speed of an identification experiment are related to the constraints of the excitation trajectory.It should be noted that the configurations of measurements must correspond to a good conditional simplified identification matrix, since the corresponding input/output represents some limitations.In literature [23], the constraints of an excitation trajectory are given 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 velocity and acceleration,  is optimal trajectory parameters,   is the available workspace of robot, and  max is the maximum joint torque.
When the robot joints repeatedly track the excitation trajectories with the PID controllers, motor current and joints positions can be sampled in the time domain.The motor current can be transformed into joint torques with a simple torque constant.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-point triple smoothing method is adopted to smooth the raw data according to the following equation [24]: where  = 3, . . .,  − 2,  = [ 1 ,  2 , . . .,   ] is the sampled raw data, and  = [ 1 , . . .,   ]  is the data for identification after preprocessing.The more the number of using ( 15) is, the smoother the curves will be.It should be noted that excessive use (15) to smooth the raw data can lead to the error of the parameters identification enlarging.Due to the shortage of speed sensor and acceleration sensor, the velocity and acceleration of joints cannot be measured directly.The calculation of the joint velocities and acceleration can be performed by analytical differentiation of the measured angles [25].The velocity and acceleration are given as follows: For the analytical differentiation, the measured encoder readings are first approximated with the least square method, as a finite sum of trigonometric functions.This approximation corresponds to frequency domain differentiation combined with frequency domain windowing.Firstly, the discrete Fourier transform of the preprocessing encoder readings is calculated and the excited frequency curves are selected by frequency domain windowing.Then, the selected curves are multiplied with the frequency response of a pure single and double differentiator, that is, multiplied with  and − 2 , with  being the frequency in radians per second.The obtained frequency spectra are then transformed back into time domain using the inverse discrete Fourier transform, yielding joint velocity and acceleration.The application of a Fourier series gives a trajectory which is continuously differentiable up to any order.And it can avoid the excitation of unknown dynamic effects.It is clear that the finite Fourier series has several advantages over the classical excitation trajectories.

Parameters Identification
To test the effectiveness and viability of our proposed method, the identification procedure was implemented on the first three joints of the ER-16 6-DOF industrial robot manipulator without payload, as shown in Figure 2. It should be noted that the values of the dynamical 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∼6 joints [26].The motion constraints of the ER-16 are given in Table 1.
We use 5th order Fourier series as the excitation trajectories in the experiment.The fundament pulsation is 0.05 Hz, resulting in a period of 25 s.The data is sampled with 1 Hz.The excitation trajectories with 0.25 Hz bandwidth are shown in Figure 3, containing 11 optimal parameters in each joint.These parameters are listed in Appendix A. Figure 4 shows the trajectory of the effect center point of the ER-16 in the workspace.Identification procedures discussed above are carried out with ICS algorithm in MATLAB 2016b programming environment on an Intel Core i7-3770 PC running Windows 7. No commercial tools are used.In order to test the performance of the proposed algorithm, the classical CS algorithm and PSO algorithm [27] are employed as the comparison.The initial parameters of the three algorithms are listed in Table 2.All the algorithms are run three times and the maximal iteration is 200.The best results would be retained.Theoretically, those results are large probability global optimum.Additionally, the search scope of the unknown parameters is also listed in Appendix B.
Figure 5 shows the evolution curves of the above three algorithms regarding (11).The figure demonstrates that the objective value increases as the generation iterates with time, gradually converging to an optimal result.Compared with CS and PSO, ICS achieves a better result with the higher objective value after 50 iterations.The optimal objective value of ICS is 29.0753, whereas those of CS and PSO are 30.2513and 31.8865,respectively.Obviously, ICS has found the optimum with the objective value equal to 29.0753 which is a 4.04% or 9.67% increase with respect to the CS or PSO.It can be seen from the result that the chaotic operator can not only avoid the traditional CS being trapped in a local optimum but also improve its robustness and efficiency.
The dynamical model of the first three joints contains 15 barycentric parameters and 15 friction parameters.These  5.0000 parameters identified by our proposed algorithm are listed in Table 3. Figure 6 compares the measured torques for the excitation trajectory with the predicted torques based on the identified dynamical parameters.The results show that the identified data by applying the three algorithms have the same trend as the measured data.Nonetheless, it turns out that the predicted torques generated by ICS approximate the actual test torques best.To verify the precision of the identified model by the above three algorithms, the correlation coefficient between the measured torques   and predicted torques   , defined as the normalized cross-covariance function, is used 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.From Table 4, we can clearly see that the correlation coefficient of the identification results produced by ICS is better than those of other algorithms.It indicates that the identified procedure based on ICS has a higher identification precision due to the strong search ability of ICS.

Validation Experiment
According to literature [9], the peaks in the prediction error would abruptly increase when the joints reverse.And the In order to validate the accuracy of the dynamical model with the nonlinear friction model, we used the linear friction model as a comparison.Similarly, the dynamical parameters of the two models are identified by the ICS algorithm.The excitation trajectories and simulation conditions are set the same as those in Section 5.The comparison of the measured torques and the predicted torques based on different friction models is given in Figure 7.As we have seen, both the predicted torques from linear and nonlinear dynamical models can match the measured torques well.Nonetheless, the nonlinear model describes the friction behavior better for the points of joint inversion, obviously in Figure 8.The nonlinear friction model can reduce the error peak, which is beneficial to the design of the controllers.

Conclusions
In this paper, an improved cuckoo search algorithm has been proposed to solve the dynamical parameters identification for a robot manipulator with nonlinear friction property.The dynamical model has been established by Newton-Euler method and processed into a linear form.Then, a nonlinear friction model is added to the joint model for realizable friction compensation at motion reversal.Based on the experimental data collected from identification experiment, we use a novel identification algorithm ICS including CS methods, the chaotic operator, and emotion operator to identify the unknown parameters of the robot model.Compared to other two identification algorithms, CS and PSO, the model generated by using our proposed algorithm matches the actual torques better.What is more, ICS has a fast convergence speed and strong search ability.Furthermore, a linear friction model is used as a comparison to test the effect of the nonlinear friction model in describing the friction characteristic.The results show that the nonlinear friction model can restrain the saltation at motion reversal, resulting in higher identification accuracy.In the future, we will attempt to study the controller design based on the identification dynamical model.

Figure 3 :
Figure 3: Optimized exciting trajectories in joint space.

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

Figure 7 :
Figure 7: Comparison of the different friction models.

Table 1 :
Motion constraints of robot.

Table 4 :
Comparison of the correlation coefficient.