Intelligent Trajectory Tracking Behavior of a Multi-Joint Robotic Arm via Genetic–Swarm Optimization for the Inverse Kinematic Solution

It is necessary to control the movement of a complex multi-joint structure such as a robotic arm in order to reach a target position accurately in various applications. In this paper, a hybrid optimal Genetic–Swarm solution for the Inverse Kinematic (IK) solution of a robotic arm is presented. Each joint is controlled by Proportional–Integral–Derivative (PID) controller optimized with the Genetic Algorithm (GA) and Particle Swarm Optimization (PSO), called Genetic–Swarm Optimization (GSO). GSO solves the IK of each joint while the dynamic model is determined by the Lagrangian. The tuning of the PID is defined as an optimization problem and is solved by PSO for the simulated model in a virtual environment. A Graphical User Interface has been developed as a front-end application. Based on the combination of hybrid optimal GSO and PID control, it is ascertained that the system works efficiently. Finally, we compare the hybrid optimal GSO with conventional optimization methods by statistic analysis.


Introduction
With the advancements in robotic technology, numerous types of robots have become involved in our daily life and help humans in many different areas. As one of the most common robots, the multi-joint manipulator robotic arm plays an important role in automotive, agriculture and bio-medical sectors due to its flexibility, robustness and accuracy [1][2][3].
The identification of the Inverse Kinematic (IK) plays an important role in the precision control of trajectory tracking [4,5]. Various IK solutions have been carried out for robotic arms [6,7]. For instance, Xu et al. [8] presented a combination brain of a computer interface and computer vision to move a robotic arm end-effector to a desired point by using a depth camera. A six degree of freedom (6DoF) robot with initialized commands from a user's brain signals combined with a point clouds model was verified with five healthy candidates without specific user training, showing acceptable accomplishment for complex tasks. Fang et al. [9] established a visual communication method using deep neural networks, in which the movements of a human arm were monitored and determined by the Denavit-Hartenberg (D-H) technique. Narayan et al. [10] presented a 5DoF robotic arm with a three-finger gripper and validated the IK in a simulation platform. Ye et al. [11] dealt with 5DoF manipulator forward-IK problems using Ferrari's and redundant Euler methods and validated them in a simulation model. Wei at al. [12] applied a neural network to a robotic arm and used environment feedback to reach a specific target point inspired by animal and human biological neural networks. They validated their approach using the penalty function to avoid the robot from reaching specific points. The results show the end-effector reached the target successfully. Ren et al. [13] developed generative neural networks to solve the IK for a robotic arm. They determined the IK by the D-H technique and Moveit, which is an application of a Robot Operating System (ROS) to control and monitor a robot.
Traditionally, the IK has been utilized to establish joint configurations of manipulators based on the end-effector position. However, the traditional IK methods cannot consider the continuity of configurations, collision avoidance and kinematic singularities that arises when attempting to follow the end-effector path [14]. In addition, solving IK problems is a difficult challenge because manipulators with more than 5DoF result in an infinite number of possible solutions for joint trajectories that determine the same position in the Cartesian space [15]. Traditional analytical solutions cannot directly calculate the one-tomany possible relationships in the Cartesian space. Therefore, evolutionary algorithms such as optimization methods are used to solve IK problems quantitatively [11,16]. For instance, Starke et al. [17] studied a mimetic evolutionary algorithm, which was a combination of the Genetic Algorithm (GA), Particle Swarm Optimization (PSO) and gradient-based optimization to address the IK solution for various industrial and anthropomorphic robots.
One of the approaches used in this paper is to combine GA and PSO in order to develop an optimal solution for the IK and Proportional-Integral-Derivative controller (PID) controller tuning. This optimization method has been presented in several works for various applications [18][19][20]. For instance, Dziwinski, et al. [21] presented a fuzzy-logic controller in which a combination of PSO and GA was used in parallel to improve PSO performance by adding crossover and mutation to the GA to avoid becoming trapped in local optima. Farand et al. [22] developed a combination of GA and PSO to reduce computational time and accuracy in comparison with other known methods such as GA and PSO for high-dimensional and complex functions.
The PID controller is one of the most common used classical control systems in different industries because of its flexibility, satisfactory results [23,24], ease of implementation in a control system and wide usage in industries [25,26]. In order to enhance the accuracy and robustness of classical PID control, one of the techniques is to combine it with optimization methods [27,28]. Belkadi et al. [29] worked on PSO with a random initial value to tune the parameters of the PID controller by minimizing the trajectory error. They verified their controller in a simulation model and compared it with conventional methods by numerical analysis. Phu et al. [30] used optimization with sliding mode control based on the Bolza-Meyer criterion to minimize the vibration effect. In another work, Suhaimin et al. [31] used a PID controller for a 5DoF robotic arm and controlled its joints for the point-to-point trajectory tracking of end-effectors.
The contribution of this paper is the development of an optimal hybrid IK and PID controller for joint trajectory tracking, using the Genetic-Swarm Optimization technique. The applicability of the proposed technique for the IK solution of end-effectors and steadystate error for control is compared with conventional optimization approaches such as the GA and PSO. In addition, a 5DoF robotic arm is selected for analysis due to its simple structure, flexible action, small volume, convenient operation and so on; this device is widely used in many fields and industries [32].
The rest of the paper is organized as follows: first, the kinematic and dynamic models of the 5DoF robotic arm are established using the D-H and Lagrangian method. Subsequently, an IK solution is determined by hybrid Genetic-Swarm Optimization (GSO) for angular trajectories of each joint reaching the target position. The joint angles determined by the IK are implemented in a closed-loop system using a PID controller, and the gains are tuned by the GA and PSO. The 3D models of the robotic arm are simulated in the Gazebo environment. Finally, a Graphical User Interface (GUI) is created to interact with the 3D model in the ROS environment.

Dynamic and Kinematic Model
The robotic arm consists of a base, four links, a wrist and gripper that are connected to each other by joints in series. Figure 1 represents a 5DoF robotic arm, in which the coordinate systems of joints and global frames are presented. There are various methods used to determine the dynamic equation of robot manipulators, such as Newton-Euler, Kane and Hamilton approaches. In this work, an energy-based Lagrangian method has been adopted to determine the relation between the torque and angle of joints; one of the advantages of the Lagrangian method is that, unlike the Newton-Euler method, it is not necessary to determine internal forces between joints; therefore, it is quicker and easier to obtain the equation of motion [33]. The Lagrangian equation is given as follows: where B i is the joint friction coefficient; L is the Lagrangian function; T i is the torque of each link, with i = 1, 2, 3, 4, 5; θ i andθ i are the angular trajectory and velocity; and E p and E k are the total potential and kinetic energies, respectively. From [34], the equations of E p and E k are determined as follows: where m i and I i are the mass and inertia of each link; g is the gravity acceleration; and (ẋ di ,ẏ di ,ż di ) is the time derivative of the centroid position of each joint, where i = 1, 2, 3, 4, 5.
According to the geometric relation, the centroid position (x di , y di , z di ) of every linkage is written as where j X ∈ 3×1 is the position of joint (i − 1) th according to the reference frame; X d i ∈ 3×1 is the position of the centroid point of link i th relative to the reference frame; R z i ∈ 3 is the rotation matrix around z-axes according to the coordinate system placed in the i th joint; and i X d ∈ 3×1 represents the centroid position of the link i th regarding the coordinate system located in the joint i th . By substituting E k and E p in the Lagrangian function, the state space dynamic is determined in range of motion condition where while one joint is moving, the other ones are fixed, which is shown as follows: where θ ∈ 5×1 andθ ∈ 5×1 are the angular rotation and acceleration; τ ∈ 5×1 is the torque vector; and M ∈ 5 is a matrix containing mass and inertia elements, which is shown as follows: where e m i i = 1, 2, 3, 4, 5 represents the mass and inertia elements, expressed as follows: where l c i is the length of the centroid position for each link and l i is the length of every link. V ∈ 5 is the centrifugal, coriolis and friction matrix and G(θ) ∈ 5 represents the gravity matrix, expressed as follows: where I 5 ∈ 5 is the identity matrix and e g i shows the elements of mass and gravity matrices, represented as follows: where g represents gravitational acceleration. Table 1 illustrates the physical features of the robotic arm's links. In the next stage, the forward kinematic based on modified D-H (mD-H) algorithms has been developed to establish the relative position of the 5DoF robotic arm end-effector to its reference frame O [35]. Table 2 represents the parameters of the mD-H.
In Table 2, α i−1 , θ i , d i and a i−1 represent the twist angle, joint angle, link offset and link length, respectively. The mD-H homogeneous transformation is expressed as follows: The transformation matrix of the end-effector is the transformation matrix from the reference frame to the last frame, which is shown as follows: where 0 1 T, 1 2 T, 2 3 T, 3 4 T and 4 5 T are the transformation matrices of each joint to its previous joint. The transformation matrix from the reference frame to end-effector is as follows: where t 1,4 , t 2,4 and t 3,4 express the end-effector position relative to the reference frame, which is shown as follows:

Optimal Inverse Kinematic
Since the number of variables is greater than the number of equations and the endeffector position is non-linear, the usage of traditional methods such as Gaussian elimination are not practical [36]. Thus, in this paper, the IK is defined as a mono-objective optimization problem. The desired position of the end-effector is set to be achieved by minimizing the objective function. In this study, the hybrid version of the GA and PSO, named GSO, is adopted to solve the IK problem, because GA is developed initially by random values due to its reliability and robust performance [37] and PSO is sufficient to find accurate results in a few iterations with low computational time. Subsequently, PSO is initialized by the results of the GA. In the GSO algorithm, the GA provides searching space and initial values for PSO to avoid becoming trapped in local optima.
The summation of squared error (SSE), which is a well-known statistic in multiple regression analyses [38], is chosen as an objective function because it shows the squared sum of residuals, which is the error between the measured and desired trajectory of the end-effector, and illustrates how close a regression line is to a set of residuals. The squaring is necessary to remove any negative signs. The objective function is given as follows: where e x , e y and e z are the errors, represented as follows: where x des , y des and z des are the desired positions of the endpoint regarding the reference frame.
x, y and z express the position of the endpoint, which are determined by the gene of the GA from Equations (18)- (20). The population structure of an iteration is represented in Figure 2. In each population, there is a gene which consists of each joint angle of the robotic arm, represented as x 1 , x 2 , x 3 and x 4 , which are θ 1 , θ 2 , θ 3 and θ 4 , respectively. In the robotic arm model, there are limitations for the angular trajectory of each joint, which create the searching space for the GA, which is as follows: θ 5 is set as 90 degrees and is not included in the design variables because it is assumed that the gripper is located at last link point down to grab the objects. The first iterations of the GA are set randomly within the searching space, as demonstrated in Equations (21)- (24). After the initialization, the objective function is calculated for each gene of the population for evaluation and sorted in ascending order. The next iterations are created by crossover and mutation. The crossover enhances the possibilities of finding the most optimum results by blending the previous iterations as children and parents using the uniform crossover operator. In addition, mutation is performed to maintain the diversity of the GA [39][40][41]. The algorithm is continued by the evaluation of each gene by determining the objective function followed by sorting in ascending order. This trend continues until the maximum iterations are reached. In the last iteration, because of the ascending sorting, the first gene is the result of GA and represents the optimum angles of joints which are needed to lead the endpoint to reach the desired position.
x ga is the output of GA which is used to create the range for the initial population of the PSO, which is shown as follows: where j stands for the number of particles in the first population and rand is the function used to generate a random value between x min and x max , which are the lower and higher bounds, given as follows [42]: x max = x ga + r where r ∈ 1×4 is a random vector between zero and one. After creating the particles of the first population, the objective function is determined for each particle to evaluate and sort in descending order. The particles of population for the next iteration are created as follows: where x i+1,j is the particle of the next iteration. v i+1,j ∈ 1×4 is a vector that represents the velocity and direction of each particle through the particle of the next iteration, which is shown as follows: where p best,i is called the best position, containing the particles that have the minimum objective function. g best is the global best, including the particles which are the minimum of the p best , which is shown as follows: where i and i max are the current and maximum number of iterations, respectively. In the first iteration, after evaluation, the minimum particle is saved as p best and g best , and the velocity is a zero vector.
In Equation (31), ω i is the inertia weight, where its adjustable value for each iteration is given by the following equation: where ω damp is the damping value for ω, set as 0.05, and c 1 and c 2 are coefficients of self and social recognition, respectively. The value of c 1 is greater than c 2 , and their summation should remain at 4 in all iterations [43].
where, a and b are the ascending and descending gains between zero and one, represented as follows:  After generating each population, the evaluation and sorting of its particles is developed. This trend is followed until the number of iterations is equal to i max . In this paper, the size of the population for GA and PSO is 40 and i max is 200 for each. Algorithm 1 and Figure 4 show the pseudo-code and flow chart of GSO. Create new iteration using crossover and mutation; 14: 15: Evaluate the population by determining the objective function; 16: 17: Sort the genes in ascending order; Create new population; 34: 35: Evaluate the particles of population; 36: 37: Set the minimum particle as the P best ; 38: 39: Set the minimum P best as g best ;

Control System and Tuning
A closed-loop control system is developed for each joint, and its parameters are adjusted by GSO to converge by adjusting the required torque toward each joint. The desired angular trajectory is determined by the IK. Figure 5 demonstrates the control system of the robotic arm.
In Figure 5, J 1 (t), J 2 (t), J 3 (t) and J 4 (t) are plants of each joint; (x t , y t , z t ) is the desired position; θ d 1 , θ d 2 , θ d 3 and θ d 4 are the desired angular trajectory determined by the IK; θ a 1 , θ a 2 , θ a 3 and θ a 4 are the actual angular trajectory for each joint; and e 1 , e 2 , e 3 and e 4 are the trajectory errors that are the difference between the desired and actual angular trajectory, given as follows: The PID controllers C 1 (s)-C 4 (s) for each joint are given as follows: C i (t) = K P e i (t) + K I e i dt + K D de i dt (40) The tuning of the PID controller is assumed to be an optimization problem, and its parameters are defined as design variables. The tuning processes are carried out by the GSO algorithm, in which the GA starts to optimize the design variables based on random initial parameters, and subsequently the algorithm is continued by PSO based on the output of the GA. Initial parameters of the first population are randomly chosen between 0 and 1, given as follows: where x is the particle of the PSO and gene of the GA in each population and j and j max are the current and maximum number of populations. After setting the initial values for particles, an evaluation is carried out based on the objective function of tuning, which is the absolute steady-state error: where θ act and θ des are the actual and desired joint trajectory, respectively. θ act is measured from a simulation model in real-time, and θ des is developed by the optimal IK. After evaluation and sorting in descending order, the particles of the population for next iteration are generated by mutation and crossover. Whenever the number of iterations meets half of the maximum iterations, the algorithm is switched to PSO. The searching space of PSO is limited around the results of the GA to lead the algorithm toward global optima, as follows: The next populations of PSO are created by the particles of the next iteration. The algorithm continues until the number of iterations reaches the maximum. The output is an optimal set for PID parameters. Figure 5. Block diagram of control system for each joint.

Results and Discussion
The optimal IK and PSO tuning of controllers was applied in the 3D simulation of a robotic arm in a 3D environment to simulate robots integrated with ROS [45]. A GUI was programmed by using Python to run the algorithms and communicate with the simulation model. In addition, by providing a camera in the Gazebo environment, it was possible to monitor the results visually and numerically, as shown in Figure 6. The desired position for the three different algorithms-i.e., GA, PSO and GSO-could be selected according to the IK method, and the actual position of end-effector, error of the actual trajectory and desired trajectory of each joint could be monitored. In addition, the controller parameters could be tuned in real time and observed. Table 3 compares the optimal results for GA, PSO and GSO for the IK solution, while f obj is the SSE for various sets of optimization parameters.
Various sets of parameters were defined to observe the influences of changes in parameters on the optimization algorithms. The mean of f obj for GSO in set 3 gso has the lowest f obj of 7.9×10 −15 , and the maximum of f obj is 4.99×10 −14 , which is the nearest value to its mean compared to other results. This causes the lowest variance of all tests. The mean of the PSO results is lower than GA, while GSO shows the minimum results, which represents a significant improvement for the results obtained by the GSO algorithm. This is due to the hybrid of the GA and PSO algorithms in series; creating the initial values of PSO based on results of the GA increases accuracy compared to using each algorithm individually. In addition, by increasing the number of iterations and particles of PSO, GSO and PSO algorithms show improvements in their results. The H-values were measured by the Kruskal-Wallis method and were 0.03, 16.07 and 15.84 for the GA, PSO and GSO respectively. The test was calculated with the assumption of α = 0.05; therefore, the critical value for this test was 5.99. Since PSO and GSO had greater H-values than critical points, there were significant differences among the groups of f obg calculated by PSO and GSO. Table 4 represents the computational time in seconds for GA, PSO and GSO, while the parameters are determined as set 3 ga , set 3 pso and set 3 gso , respectively. The mean computational time of PSO was less than GA and GSO by 5.82 s and 1.37 s, respectively. Although the computational time of PSO was the lowest, the combination of GA and PSO reduced the computational time consumption significantly compared to GA by 4.45 s. By considering the value of f obj of GSO in Table 3 and the computational time, the usage of GSO for IK solution can be seen to have resulted in significant improvements in accuracy and computational time consumption.
In order to test the IK results solved by GSO in the robotic arm model in the Gazebo environment, nine desired position coordinates were expressed. Table 5 illustrates the coordinates of the nine target positions and the angles for each joint.  Figure 7 shows the objective function f obj for three ways of tuning PID parameters with 400 iterations, and the parameters of GSO were the same as set 3 gso in Table 3. From the results, it can be observed that GSO converges faster than GA and PSO, because by establishing the angle of joints for the desired points, the angular trajectories are exported to the control system and then tuned by GSO. The performance of the closed-loop control system is validated for each joint, in which the angular trajectories solved by the IK are set as desired (θ d i ). Table 6 represents the tuned PID parameters.   Figure 8 and Table 7 compare the angular trajectories and average errors tuned by the GA, PSO and GSO while the end-effector moved to the desired positions and its PID parameters were tuned by the GA, PSO and GSO. Figures 9 and 10 show the error and angular velocity for each joint while the end-effector tracked points A, B, C, and D. In Figures 8 and 9, overshoot can be observed when the joints are changing trajectories. In Figure 10, there are fluctuations when there is a change in position of the joints and they are tracking each point. The rest of the graph levels off at zero. This shows the stability of the control system, because while joints do not move and are in a stable situation, the joints do not shake. In Table 7, AE GA , AE PSO and AE GSO are averages of the SSE for GA, PSO and GSO. The actual trajectory of each joint shows the significant effects of GSO compared to PSO and the GA due to the initialization of PSO by results of the GA for GSO, which allows the algorithm to find the global optima more accurately. From AE GA , AE PSO and AE GSO , it can be concluded that the GSO resulted in a lower value because of its efficiency in finding precise results and lower chance of becoming trapped in local optima.

Conclusions
This paper presented a hybrid optimal IK solution of a 5DoF robotic arm to determine the joint trajectories based on its end-effector position. The Denavit-Hartenberg method was used to establish the kinematic and was solved by GSO, which is a combination of the GA and PSO. The trajectories were implemented with PID control as desired trajectories for each joint. The tuning of PID parameters was presented as optimization problem and carried out by GSO. A GUI was created to operate and visualize the performance of the robotic arm in a virtual environment. This method addresses the issue of finding one-to-many possible solutions of IK for a 5DoF robotic arm to control each joint efficiently and precisely to determine end-effector positions in Cartesian space.
The results show that GSO has a lower average error for each joint than PSO and GA. For instance, for joint 4, the average error for one specific path for GSO is 19.33% and 22.7% less than PSO and the GA, respectively. These results show that initialization particles for PSO by using the GA can give more accurate results and avoid algorithms becoming trapped in local optima. In addition, the mean computational time for GSO is lower than the GA by 4.45 s and higher than PSO by 1.37 s. Therefore, GSO is a sufficient algorithm for IK solution for robotic arms.
This method can be used for any robotic arms to control their end-effector. The limitation of this work is that we did not apply the hybrid proposed method to a real robotic arm. In addition, the target position of the end-effector was chosen by the user. Therefore, in future work, this position can be issued by sensors such as depth camera or tag marker measurement algorithms. Furthermore, the proposed IK and control system developed by the GSO algorithm was validated in the 3D simulation environment of Gazebo, and the effect of sensor noises was not considered; this can be covered in the future work. In addition, the proposed method can be tested for applications in which the accurate position of end-effectors is needed, such as welding, material handling and thermal spraying or any other industrial applications.