Performance Analysis of Numerical Integration Methods in the Trajectory Tracking Application of Redundant Robot Manipulators

Differential kinematic is one of the most important solution methods in robot kinematics. The main advantage of the differential kinematic method is that it can be easily implemented any kind of mechanisms. Also, an accurate and efficient kinematic based trajectory tracking application can be easily implemented by using this method. In differential kinematic method, we use Jacobian as a mapping operator in the velocity space. Inversion of Jacobian matrix transforms the desired trajectory velocities, which are the linear and angular velocities of the end effector, into the joint velocities. The joint velocities are required to be integrated to obtain the pose of the robot manipulator. This integration can be evaluated by using numerical integration methods, since the inverse kinematic equations are highly complex and nonlinear. Therefore, the performance of the trajectory tracking application of the robot manipulator is directly affected by the chosen numerical integration method. This paper compares the performances of numerical integration methods in the trajectory tracking application of redundant robot manipulators. Several widely used numerical integration methods are implemented into the trajectory tracking application of the 7-DOF redundant robot manipulator named PA-10 and simulation results are given.


Introduction
The configurations of redundant robot manipulators offer the potential to overcome many difficulties by increased manipulation ability and versatility [1,2].A desired trajectory can be tracked in many different configurations of redundant robot manipulators by using the extra manipulation ability.This feature can be easily implemented into the many robotic applications such as obstacle avoidance, singularity avoidance, complex manipulation etc. [3, 4 and 5].Redundant robots are also frequently used in complex industrial applications, service robots and humanoids [6].However the redundant robot manipulators have many advantageous, they require quite complex control structures and suffer from singularity problem.
A fundamental research task of redundant robot manipulation is to find out the appropriate way to control the system of redundant robot manipulator in the work space at any stage of the trajectory tracking.This control can be achieved by using dynamic or kinematic model based control.However, dynamic model based control gives us more realistic results than kinematic based control, it requires highly complex control structures.Therefore, it is generally used to control the few degrees of freedom (DOF) robot manipulators in the laboratory researches [7,8].Kinematic model based control can also be used to solve the problem of redundant robot manipulation.In the kinematic problem, the equation of motion is obtained without using the forces and torques which cause the motion and more elegant control structures can be obtained due to the nature of the kinematic problem.However kinematic model based control is not as realistic as dynamic model based control, it has a quite simple control structure and it can be easily implemented to the robotic problems which do not require force and torque controls.Therefore, kinematic model based control is widely used in many robotics applications such as industrial robotics [9,10].
Differential kinematic is one of the most important solution methods in robot kinematics [11,12].The main advantage of the differential kinematic is that it can be easily implemented any kind of mechanism.Also, an accurate and efficient kinematic based trajectory tracking applications can be easily implemented by using this method [13].Jacobian is used as a velocity mapping operator which transforms the joint velocities into the Cartesian linear and angular velocities of the end effector.A highly complex and nonlinear inverse kinematic problem of redundant robot manipulators can be solved numerically by just inversing the Jacobian matrix operator.However, differential kinematic based solutions can be easily implemented any kind of mechanisms, it has some disadvantages.The first one is that differential kinematic based solutions are locally linearized approximation of the inverse kinematic problem [14].Thus, we can only obtain the approximate solutions of the inverse kinematic problem by using this method.Although the solution results of this method are not real, the approximation results are generally quite sufficient for small joint velocities.The second disadvantage of this method is that it has a heavy computational load and big computational time because of numerical iterative approach [7].Requirement of the inversion of Jacobian matrix operator directly affects the computational efficiency in the inverse kinematics.Moreover, this requirement also causes the singularity problem.Singularity is one of the most important problems in robot kinematics [15].At or around the singular configurations of robot manipulator, Jacobian transformation generates high joint velocities which results in instability and large errors in the task space.Several studies have been published on the problem of kinematic singularity in the literature.In general, there are four main techniques to cope with the kinematic singularity problems.These are avoiding singular configuration method, robust inverse method, a normal form approach method and extended Jacobian method [16][17][18][19].However given techniques have some disadvantages which include computational load and errors.And the last disadvantage of the differential kinematics method is that, it requires numerical integration which suffers from numerical errors, to obtain the joint positions from the joint velocities [20].The numerical integration of joint velocities to compute joint positions causes a numerical drift which in turn corresponds to a task space error [21,22].An effective inversion of differential kinematics mappings can be realized by adopting the so-called closed-loop inverse kinematics algorithms which are based on the use of a feedback correction term on the task space error [23].However the drift-phenomena can be overcome by using the closed-loop inverse kinematic algorithm, the performance of the algorithm is still extremely affected by the chosen numerical integration method.
In numerical integration, we calculate the analytical integral approximately by using numerical techniques.This computation can also be called as quadrature [24].There are a wide range of quadrature methods available in the literature [25][26][27][28].Accuracy and performance are main requirements in the quadrature algorithms.In general, better accuracy can be obtained by increasing the computational complexity and computational load.Therefore, this case, in general, results as a trade-off between the accuracy and computational performance in the numerical integration algorithms.
In this paper, a performance analysis of the numerical integration methods in the trajectory tracking application of the redundant robot manipulators is presented in details.Three single-step numerical integration methods which are Euler integration, Runge-Kutta 2 and Runge-Kutta 4 and also three multi-step numerical integration methods which are Predictor & Corrector, Adams-Bashforth and Adams-Moulton methods are implemented into the differential kinematic based solution of the trajectory tracking application of the redundant robot manipulator.These methods are compared with respect to computational efficiency and accuracy.Simulation results of the trajectory tracking are given in section V.This paper is also included the differential robot kinematics in section II, numerical integration methods in section III, trajectory tracking algorithms in section IV. Conclusions and future works are drawn in the final section.

Differential Robot Kinematics
It is so hard, even impossible to find the analytical solutions of the inverse kinematic problem of the redundant robot manipulators except the limited special structures or very simple mechanisms.Therefore, differential kinematic based solution of the inverse kinematic problem of the redundant robot manipulators is widely used [29].Differential kinematic based solution can be formulated as follows, Let's p define the pose, which includes the position and orientation, of the end effector.It can be formulated in terms of the joint angles as follows, where indicates the joint angles and n is the number of the joints.If we differentiate it then, we get,   where is the Jacobian matrix and indicates the end effector's velocities.As it can be directly seen from the equation (2), Jacobian is a mapping operator which transforms the joint velocities into the Cartesian velocities.Joint velocities can be obtained by using the inversion of the Jacobian matrix operator as the following, where   g J q indicates generalized inverse of the Jacobian matrix.In the introduction part, we mentioned about that differential kinematic based solution is a locally linearized approximation of the inverse kinematics solution [30].It can be directly seen by using the Taylor expansion of the forward kinematics equation.It can be formulated as follows, If we assume that Δq is relatively small then, we can omit the high order terms of the Taylor expansion.Then, we get the approximate solution as the following, F q F q F q Δq F q q Δq q (5) As it can be directly seen from the equation ( 5) differential kinematic based solution is locally linearized approximation.This equation is admissible if the Δq is small enough and the higher order terms of the Taylor series can be neglected.Also the equation ( 5) equals to a simply first order Euler integration.The joint velocities should be integrated to obtain the joint positions by using the following equation, As it can be seen from the equation ( 6), the integrand is highly complex and nonlinear.Analytical solution of this integration, in general, is so hard, even impossible.Therefore, numerical integration methods should be used to solve this problem.The numerical integration solution techniques of the equation ( 6) will be discussed in the next section.
Note that, we derived the Jacobian operator analytically in this section.However, different methods can be used to find the Jacobian operator.Several methods to obtain the Jacobian operator can be found in [31][32][33] 3. Numerical Integration of The Joint Velocities Iterative solution methods are generally used to solve the inverse kinematic problem of the redundant robot manipulators because of highly complex and nonlinear inverse kinematic equations.The joint angles are obtained by numerically integrating the joint velocities.Therefore, the chosen numerical integration method extremely affects the computational efficiency and accuracy of the differential kinematic based trajectory tracking algorithm.
Here, several numerical integration methods are introduced.These integration methods can be divided into two main different approaches which are single-step and multi-step numerical integration methods.The formulations of these integration methods are as the following, [25][26][27] A. Single-Step Numerical Integration Methods

 Explicit Euler Integration Method:
Explicit Euler integration is the simplest numerical integration method.It is simply the first order Taylor expansion and formulated as follows where The strength of this method is that it can be easily implemented and also it has a very computationally light algorithm.However, the accuracy of this method is quite poor.It requires small sampling rates to obtain stable and accurate results.
 Runge-Kutta 2 Method: Runge-Kutta methods are one of the most widely used numerical integration methods.The formulation of the second order Runge-Kutta numerical integration method is as follows, where This method requires two calculations of the generalized inverse Jacobian for each step, so that the computational load of this method is higher than Explicit Euler integration.However it requires extra computation, it gives more accurate and stable results than Explicit Euler integration.

 Runge-Kutta 4 Method:
The formulation of the fourth order Runge-Kutta numerical integration method is as follows, where This method requires four calculations of the generalized inverse Jacobian for each step, so that the computational load of this method is higher than Runge-Kutta 2. This extra computation improves the numerical integration results and the solutions, which are more accurate and stable than Runge-Kutta 2 based solutions, can be derived by using this method.As the order of the Runge-Kutta numerical integration method increases, more accurate and stable results are obtained by using more computationally complex algorithms.

 Euler Trapezoidal Predictor & Corrector Method:
Euler Trapezoidal Predictor & Corrector method is an algorithm that proceeds in two steps.First, the prediction step calculates a rough approximation of the desired quantity.In this step simple Explicit Euler integration is used.Second, the corrector step refines the initial approximation using another means.In this step trapezoidal integration is evaluated by using the predicted and current joint velocities.The formulation of this method is as follows, where   ˆk t q  is the predicted joint velocity vector in which Euler Trapezoidal Predictor & Corrector method also requires two computation of the generalized inverse of Jacobian operator so that the computational load increases.It gives more accurate and stable results than Explicit Euler Integration.

 Adams-Bashforth Method (Fourth Order):
Adams-Bashforth is a widely used multi-step explicit numerical integration method.It can be formulated as follows, If If If where Several different order of this method can be used to obtain the numerical integration.Here, we used the fourth order Adams-Bashforth numerical integration algorithm which is the most widely used one.This method requires three step backward values of joint velocities and one computation of generalized inverse of Jacobian operator for each step so that it is a computationally light numerical integration method.As the order of the method increases, more accurate and stable results are obtained.

 Adams-Moulton Method (Fourth Order):
Adams-Moulten is also a widely used multi-step implicit numerical integration method.Here, Adams-Bashforth algorithm is used in the numerical integration of the predicted joint velocities and Adams-Moulton algorithm is used in the numerical integration of the corrected joint velocities.It can be formulated as follows, If If where   ˆk t q  is the predicted joint velocities in which (24) Here, we used the fourth order Adams-Moulten numerical integration algorithm which is the most widely used one.This method requires two step backward values of joint velocities and one step forward predicted joint velocities.It also requires two computations of generalized inverse of Jacobian operator for each step so that computationally load increases.This extra computation improves the numerical integration results and the solutions which are more accurate than Adams-Bashforth based solutions, can be derived by using this method.As the order of the method increases, more accurate and stable results are obtained.

Trajectory Tracking Application
The trajectory tracking application of the redundant robot manipulator is implemented by using the following two simulink block diagrams which are shown in figures 1 and 2. The first one shows us the trajectory tracking application by using the explicit numerical integration methods which are Explicit Euler Integration, Runge-Kutta 2, Runge-Kutta 4, and Adams-Bashforth methods.In this application, a desired trajectory is generated for the end effector of the robot arm in the Desired Trajectory block and it is transferred to the Jacobian block.In the Jacobian block, the joint velocities are obtained by using the velocity mapping.Inputs of the Jacobian block are tip point velocities, error and positions of the robot manipulators.The tip point velocity input, which are linear and angular velocities of the robot manipulator, is the desired trajectory.The error input is used to obtain the closed loop inverse kinematic solution, which can be shown in figure 1. Closed loop inverse kinematics is used to cope with the drift phenomena [12].And the position input is used to obtain the Jacobian iteratively in each step of the simulation.Joint velocities are generated in the Jacobian block and then, they are transferred to the Numerical Integration block.In the Numerical Integration block, explicit numerical integration methods, which are Explicit Euler Integration, Runge-Kutta 2, Runge-Kutta 4, and Adams-Bashforth, are used to obtain the joint angles.The obtained joint angles are transferred to the Forward Kinematics block.In the Forward Kinematics block, we obtain the pose of the robot manipulator's and the positions of the each robot manipulator's joints.The positions of each robot manipulator's joints are required to obtain the Jacobian operator iteratively.Therefore, these joint positions are transferred to the Jacobian block.Also, the pose of the robot manipulator is obtained in the Forward Kinematics block.The pose of the robot arm is required to obtain the Jacobian operator iteratively and also the closed-loop kinematic structure.Therefore, the pose of the robot arm is also transferred to the Jacobian block.The error of the trajectory tracking is calculated by using the following formulas, where

Simulation Results
A redundant anthropomorphic robot arm structure that is shown in figure 3 is used for the simulation studies.Anthropomorphic robot arm structure is widely used in many robotics applications such as industrial applications (i.e.welding assembly etc.) and non-industrial applications (i.e.rehabilitation, human motion assistance, etc.) [34,35].
The robot structure which is shown in figure 3 is obtained by using the Mitsubishi PA-10 anthropomorphic robot arm.Mitsubishi PA-10 robot arm features an articulated arm with 7-DOF for high flexibility.It spreads a wide range area in many robot applications.Simulation studies of the trajectory tracking application are performed by using Matlab and the animation applications are performed by using virtual reality toolbox (VRML) of Matlab which can be seen in figure 3.As it can be seen from the figures 4, 5 and 6 a satisfactory trajectory tracking application is implemented for the desired circular path.The circular path is tracked precisely as shown in figure 4. Also, both of the joint positions and velocities are smooth during the trajectory tracking application as shown in figures 5 and 6.In this trajectory tracking application, Runge-Kutta 4 is used for the real time numerical integration with 100 t   ms sampling rate.There are two main requirements, which are accuracy and computational efficiency, in the real time numerical integration applications.The computational efficiency performances of the proposed numerical integration methods in the trajectory tracking application can be seen in the figures 7 and 8.As it can be directly seen from the figures 7 and 8 the most computationally efficient method is Explicit Euler Integration and the least computationally efficient method is Runge-Kutta 4. The computational efficiency results of the numerical integration methods are quite similar.However, these small performance differences of the numerical integration methods may drastically affect the performance of the real time application.The computation time is evaluated by using Matlab's tic-toc commands.Running environment is as follows, Accuracy is the other important requirement in the numerical integration applications.The accuracies of the proposed numerical integration methods are analyzed in the trajectory tracking application and the simulation results are given in the figures 9,10,11,12,13

 Explicit Euler Integration
Figure 9.Total orientation and position errors (radian and meter) of the end effector for the Euler integration method and the sampling rates are (a)     Adams-Bashforth Methods (Fourth Order) Figure 13.Total orientation and position errors (radian and meter) of the end effector for the fourth order Adams-Bashforth method and the sampling rates are (a) As it can be seen from the figures 9, 10, 11, 12, 13 and 14, the most accurate method is Runge-Kutta 4 and the least accurate method is Explicit Euler Integration.Explicit Euler Integration based solution gives fairly poor accuracy results in the trajectory tracking application.The results of this method can be seen in the figure 9. Small sampling rates, which increase the computational loads of the trajectory tracking algorithm, should be used to improve the accuracy of the numerical integration method and also to avoid the instability.As it can be seen from the figure 9 (a), Explicit Euler Integration based solution makes the system instable and errors get higher when the sampling rate is 1 t   sec.The Euler Trapezoidal Predictor & Corrector method also gives poor accuracy results in the trajectory tracking application.The result of this method can be seen in the figure 12.However, it has poor accuracy as Explicit Euler Integration; it is stable when the sampling rate is 1 t   sec.The other predictor & corrector based numerical integration method is Adams-Moulten.This method uses Adams-Bashforth algorithm for the prediction and Adams-Moulten algorithm for the correction.As it can be seen from the figures 13 and 14, both of these methods give quite accurate and stable results in the trajectory tracking application.However, Adams-Moulten based solution is more accurate and stable than Adams-Bashforth based solution.At last, Runge-Kutta based numerical integration results can be seen in the figures 10 and 11.As it can be seen from the figures, Runge-Kutta based numerical integration based solution gives very accurate and stable results.Also, the fourth order Runge-Kutta method (Runge-Kutta gives the most stable and accurate numerical integration results.Runge-Kutta based numerical integration gives quite accurate results when the sampling rates are high.For instance, position errors of the end effector are about At or around the singular configurations of robot manipulator, Jacobian transformation generates high joint velocities which results in instability and large errors in the task space.The main reason of the large errors in the task space is that the differential kinematic based solution is the locally linearized approximation of the inverse kinematic problem and it is valid only for relatively small joint velocities.This case is explained clearly in the section of differential robot kinematics.Therefore, whatever the numerical integration method used task space errors will be increased at or around the singular configurations of the robot arm.However, accuracy results will be different in terms of the used numerical integration methods.The position and orientation errors around the singular configurations of robot arm can be seen in the figures 15 and 16.Figures 15 and 16 show the trajectory tracking errors around the singular configurations of robot arm using single and multi steps numerical integration methods respectively.As it can be seen from the figures 15 and 16, Runge-Kutta 4 numerical integration method gives the most accurate and stable results around the singular configurations of robot arm.However, the errors are still getting high around the singular configurations of robot arm even if Runge-Kutta 4 is used.Since, singularity does not depend on the used numerical integration method.It actually based on the configurations of the robot manipulator.Therefore, we should use different techniques such as avoiding singular configuration, robust inverse etc., to cope with the kinematic singularity problems.The results of the avoiding singular configuration method are shown in the figures 9, 10, 11, 12, 13 and14.In this method, a trajectory, which does not suffer from kinematic singularity problems, is generated and the kinematic problem is solved by using this singularity free trajectory.Also, there are several robust inverse kinematic algorithms to cope with the kinematics singularity problem [36,37].Damped least squares method is one of the most widely used robust inverse kinematics algorithm [38,39].If we apply the robust inverse kinematics algorithm using the damped least squares method then, we get the errors around the singular configurations of robot arm as shown in the figures 17 and 18.

 Euler Trapezoidal Predictor & Corrector Method
As it can be seen from the figures 17 and 18, robust kinematics algorithm using the damped least squares method gives satisfactory trajectory tracking results around the singular configurations of robot arm.Moreover, Runge-Kutta 4 based numerical integration methods gives more accurate results around the singular configurations of robot arm.

Conclusion
In this paper, we analyzed the performance of numerical integration methods in the trajectory tracking application of redundant robot manipulators.The performance of the trajectory tracking algorithm is drastically affected by the chosen numerical integration method.For instance, more accurate and more computationally efficient trajectory tracking algorithm can be obtained by changing the numerical integration methods.Even, the trajectory tracking algorithm may become instable due to the chosen numerical integration method.Here, we compared six different numerical integration methods with respect to computational efficiency and accuracy.Among these methods, Runge-Kutta and Adams methods give satisfactory results.When we compare the Runge-Kutta and Adams methods, Runge-Kutta based algorithms give more accurate and stable results however; they require extra computation.Thus, the Adams methods are more computationally efficient than Runge-Kutta methods.In the trajectory tracking application, Runge-Kutta based algorithms give quite satisfactory results when the sampling rates are high.As the sampling rates increase, computational load of the trajectory tracking algorithm decreases.However Runge-Kutta based algorithms require extra computations and have high computational load, the satisfactory results at high sampling rates may reduce even eliminate the computational efficiency disadvantages' of this method.

Figure1.
Figure1.Simulink Block Diagram of Trajectory Tracking Simulation Application

Figure 3 .
Figure 3. PA-10 Robot arm animation in virtual reality toolbox (VRML) A circular trajectory tracking application is implemented by using the proposed kinematic control structure.Figures 4, 5 and 6 show the trajectory tracking application results.The desired and actual paths are shown in the figure 4. Joints positions and velocities are given in the figures 5 and 6.

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Path tracking on the coordinates x and y

Figure 7 .
Figure 7. Simulation times of the single-step numerical integration based solutions (sec)

Figure 8 .
Figure 8. Simulation times of the multi-step numerical integration based solutions (sec)

Figure 10 .
Figure 10.Total orientation and position errors (radian and meter) of the end effector for the Runge-Kutta 2 and the sampling rates are (a)

Figure 11 .
Figure 11.Total orientation and position errors (radian and meter) of the end effector for the Runge-Kutta 4 and the sampling rates are (a)

Figure 12 .
Figure 12.Total orientation and position errors (radian and meter) of the end effector for the Euler Trapezoidal Predictor & Corrector method and the sampling rates are (a)

Figure 14 .
Figure14.Total orientation and position errors (radian and meter) of the end effector for the fourth order Adams-Moulten method and the sampling rates are (a)1 t   s (b) 100 t   ms (c) 10 t   ms (d) 1 t   ms

10 
for Runge-Kutta 4 when the sampling rate is 100 ms .High sampling rates improve the computational performance of the algorithm.

Figure 15 .Figure 16 .Figure 17 .Figure 18 .
Figure 15.Total orientation and position errors (radian and meter) of the end effector around the singular configurations of robot arm for the Explicit Euler, Runge-Kutta 2 and Runge-Kutta 4 numerical integration methods.Sampling rate is 100 t   ms