Meta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking Capability Under Joint Limits

MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX ISSN: 1803-3814 (Printed), 2571-3701 (Online) https://doi.org/10.13164/mendel.2022.1.041 Meta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking Capability Under Joint Limits Ganesan Kanagaraj1,  , S A R Sheik Masthan2, Vincent F Yu3 1,2Department of Mechatronics Engineering, Thiagarajar College of Engineering, Anna University, Madurai, Tamil Nadu, India 3Department of Industrial Management, National Taiwan University of Science and Technology, Taipei, Taiwan [email protected]1,  ,sa[email protected]2,[email protected]3 Abstract In robot-assisted manufacturing or assembly, following a predefined path became a critical aspect. In general, inverse kinematics offers the solution to control the movement of manipulator while following the trajectory. The main problem with the inverse kinematics approach is that inverse kinematics is computationally complex. For a redundant manipulator, this complexity is further increased. Instead of employing inverse kinematics, the complexity can be reduced by using a heuristic algorithm. Therefore, a heuristic-based approach can be used to solve the inverse kinematics of the robot manipulator end effector, guaranteeing that the desired paths are accurately followed. This paper compares the performance of four such heuristic-based approaches to solving the inverse kinematics problem. They are Bat Algorithm (BAT), Gravitational Search Algorithm (GSA), Particle Swarm Optimization (PSO), and Whale Optimization Algorithm (WOA). The performance of these algorithms is evaluated based on their ability to accurately follow a predefined trajectory. Extensive simulations show that BAT and GSA outperform PSO and WOA in all aspects considered in this work related to inverse kinematic problems. Keywords: Inverse Kinematics, Redundant Robot Manipulator, Path Tracking, Meta-Heuristic Algorithm, BAT, Particle Swarm optimization, Gravitational Search, Whale Optimization. Received: 05 April 2022 Accepted: 15 June 2022 Online: 24 June 2022 Published: 30 June 2022 1 Introduction Recent increase in the use of robot manipulators in industries has attracted a great deal of attention to control the robot manipulator. Their wide applications in industrial systems include welding, painting, assembly, etc. A welding torch or a paint sprayer will be attached as a tool in the end effector of the robot manipulator. This tool has to accurately follow a predefined path in order to perform a particular task. Robot manipulators are made of several connected links. These links are to be controlled accurately in order to follow the given reference trajectories precisely. However, it is typically not tractable to control them accurately because of their high nonlinearity and unmodeled uncertainties. Thus, many works have been carried out to resolve such difficulties. Several methodologies have been used for solving the inverse kinematics problem. Quaternion transformation approach was proposed to solve the inverse kinematic problem wherein solution for a general 7-link 7R mechanism is presented [17,18,19]. [26] proposed a detailed derivation of inverse kinematics using exponential rotational matrices by breaking the 6R-chain in the middle to form two open 3R-chains. Later [53] used the same quaternion approach instead of the regular Newton-Euler and Lagrange method as it simplified the modelling of the kinematics and dynamics of rigid multi-body systems. Similar approaches were followed for trajectory tracking also. Jacobian matrix-adaption method was employed to overcome the two major limitations in Jacobian-matrix-pseudo-inverse (JMPI) for tracking control of robot manipulators [7]. A decentralized control strategy with finite-time convergence is developed for the trajectory tracking of a space manipulator in [41]. In this work, the robot manipulator is considered as a number of decoupled subsystems. A similar decoupling mechanism was proposed in [50] for aircraft assembly using a multi-objective posture optimization algorithm. For aerial manipulation, a multi-stage Model Predictive Control based approach was proposed in which was verified using a 3 degrees of freedom manipulator [15]. Several iterative approaches were also experimented for developing a controller for robot manipulator to follow a predefined trajectory. Iterative Learning Control was developed to identify and calculate the robot kinematic parameters and proposed an algorithm for the accurate path tracking of industrial robots [52]. An adaptive control method for trajectory tracking of robot manipulators, based on new neuro-fuzzy modelling was proposed in [43]. The proposed control 41 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX scheme used a three-layer neural fuzzy network to estimate the system uncertainties, and its performance was compared with the conventional computed torque PD control. In [48], an adaptive output feedback tracking controller was proposed to prove the uniform global stability for the robot dynamic model with unknown parameters. [8] used Differential Evolution algorithm for training the neural network to obtain the kinematic modeling of robot manipulators. [45] used neural network to track the trajectory adaptively. A radial basis function network is investigated to the joint position control of an n-link robot manipulator. Andreev and Peregudova investigated the trajectory control using uniform asymptotic stability in closed-loop system by using the dynamic position-feedback controller with feedforward [1]. A two-link planar elbow robot manipulator was used to illustrate the results. Baek et al., presented a practical Adaptive Time Delay Control scheme which was applied to robot manipulators to achieve good tracking performance with tolerant fluctuation and fast convergence speed [3]. Several control algorithms were proposed for trajectory tracking mechanisms in industrial robots. To name a few, Discrete-Time Nonlinear Optimization control [21], Sliding mode controller [30], adaptive control [49]. A sinusoidal-input describing function model along was created along with a controller for a twolink robot manipulator in [9]. A trajectory algorithm using artificial neural network and Kalman filter was proposed by [27]. The efficiency of the algorithm was verified by simulated using a 3-link Manipulator. To overcome the computational difficulties and approximations involved with the analytical methods, a machine learning based algorithm was proposed to predict the inverse kinematic solutions for parallel manipulators [44]. The computational complexity of the conventional analytical and other Jacobian based inverse kinematics led to the use of heuristic and meta-heuristic-based approach for this inverse kinematics problem [5]. Many swarm intelligent and meta-heuristics algorithms were employed to solve various optimization problems in the field of engineering. A new mutated genetic algorithm was employed to solve the problem of independent job scheduling in grid computing [51]. Study on the inner dynamics of PSO algorithm using network visualization showed the self-adaptive approaches of PSO [35]. The versatility of these algorithms had led to the development of hybrid algorithms. To name a few, genetic algorithm with fuzzy for multi objective optimization [40], Particle Swarm Optimization with clustering algorithm for system modeling and identification [23], facility location-network design model using Firefly and Invasive Weed Optimization based fuzzy system [39]. For the inverse kinematics issue considered in this paper, various meta-heuristic algorithms such as Genetic Algorithm [34,6], Particle Swarm Optimization Algorithm [10,15], Cuckoo Optimization Algorithm [4], Genetic Algorithm, Gravitational Search Algorithm [2], Artificial Bee Colony [14], Fire Fly Algorithm [38], Modified Firefly [22], Bat Algorithm [28], were employed. Heuristic and meta-heuristic algorithm-based motion planning and tracking were also employed. [25, 29,42,46] proposed an adaptive genetic algorithm for tracking the trajectory of a robot manipulator. Online optimizations method was proposed by [16] for trajectory tracking to overcome the limitation of the traditional methods. Literature survey shows the usage of meta-heuristicbased approaches for inverse kinematic problem. Since, the meta-heuristic approach does not include a Jacobian matrix and there is always a solution for forward kinematics, there are no singular configurations. At the same time, literature review shows that the usage of meta-heuristic approach for the robot manipulator trajectory following is very limited. The performance of these meta-heuristic algorithms is measured based on the faster convergence rate. From the survey, it was found that the usage of meta-heuristic-based approaches for inverse kinematics problem is limited and it was not extended to actual robot manipulators. This motivated us in using proposing a meta-heuristic-based approach to solve the inverse kinematic problem for an industrial robot manipulator. This paper’s contributions are summarized as follows. Four algorithms, viz.: BAT, GSA, PSO and WOA are experimented in this research work for solving the inverse kinematics problem for robot manipulators. Heuristic approach is proposed as it does not include a Jacobian matrix and there is always a solution for forward kinematics, there are no singular configurations as a consequence of inverse kinematics. The four proposed algorithms are compared based on their capability to generate solution that helps the robot manipulator to closely follow the pre-defined trajectory. Four predefined trajectories, viz.: linear, curvilinear, saw tooth and rose curve, are used to test the trajectory following capability of the heuristic algorithm. Wilcoxon test is employed to measure the performance of the trajectory following capability of the proposed algorithms based on three parameters viz., Minimum Average Error, Fast Convergence and Minimal variation in joint angles. The results are simulated and discussed. The remainder of the paper is structured as follows. A general overview of the problem is described in Section 2. Section 3 explains the proposed methodology and how it is used to carry out the simulation for trajectory following. Section 4 discusses about the experimental setup and the results obtained from the experiments. The paper is concluded with the key findings in section 5. 2 Problem Description The position of the end effector of a robot manipulator can be controlled or varied by changing the link angles of the robot manipulator. Calculating these link 42 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX Kanagaraj et al.: SettingsMeta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking ... Table 1: Input and output of the algorithms. Input 1. jth coordinate point (xd, yd, zd)jin the trajectory 2. link angles (θ1, . . . θl, . . . θnL)j−1generated for reaching the to (j−1)th point Output link angles (θ1, . . . θl, . . . θnL)jfor reaching the jth coordinate point (xa, ya, za)jin the generated trajectory angles for a particular end effector position in space is done generally by inverse kinematics. Inverse kinematics, in general, is computationally complex. The increased number of links in a redundant manipulator further increases this complexity [11]. The computational complexity of this link angle calculation can be reduced by employing heuristic algorithm instead of inverse kinematics. So, heuristic-based approach can be used to solve the inverse kinematics of the robot manipulator end-effector, which guarantees that the desired paths are followed precisely. 2.1 Problem Statement The efficiency of any robot manipulator in closely following the predefined trajectory can be measured in terms of its deviation from the predefined path or trajectory. This difference or deviation of the end effector from the predefined path at any point in the trajectory can be mathematically stated as: Error(ϵ) = q(xd−xa)2+ (yd−ya)2+ (zd−za)2 (1) Here, (xd, yd, zd),(xa, ya, za) denotes the Cartesian coordinates of desired and actual position of the end effector respectively. The problem is to find the link angles of each joint such that the error is minimum. Therefore, the problem can be transformed into an equivalent optimization problem as: Minimize Error(ϵ) subjected to the condition that (2) θl,min ≤θl≤θl,max Here, θlrepresents the lth link angle, θl,min and θl,max are the minimum and maximum value of the lth link angle, l= 1,2, . . . , nL with nL represents the number of links. So, based on the above problem discussion, the objective of this research work is to experiment the proposed heuristic algorithms (BAT, GSA, PSO & WOA) for calculating the link angles of a robot manipulator to accurately follow a predefined trajectory. 3 Proposed Methodology Four algorithms viz.: BAT, GSA, PSO and WOA, are employed in this research work to achieve the proposed objective. The efficiency of these algorithm is tested based on its closeness in following the generated trajectories. Table 1 shows the input and output of the algorithms. Here, j= 1,2, . . . , nPT and nP T is the number of points in the generated trajectory. The algorithm iterates repeatedly and gives the link angles (θ1, . . . θl, . . . θnL)jof the robot manipulator to the reach the point that is given as input. Here l= 1,2, . . . , nL and nL is the number of links in the robot manipulator. When these link angles (θ1, . . . θl, . . . θnL)jare applied to the robot manipulator, it will orient in a particular configuration there by reaching a point (xa, ya, za)jThe difference between (xd, yd, zd)jand (xa, ya, za)jis the error for jth point in the trajectory. For smooth movement of the robot manipulator, there should be minimal changes in the link angles and there should not be abrupt variations in the link angles. To ensure this, the link angles generated for the previous point in the trajectory is also given as input to the algorithm. 3.1 Particle Swarm Optimization (PSO) PSO algorithm was introduced by Kennedy and Eberhart [13]. It uses swarm intelligences for solving problems. It was inspired by the social behavior of the birds and fishes in finding their prey. A comprehensive study with the applications of PSO were presented in [24]. In PSO, each particle searches for the optimal solution thereby moving with a certain velocity. Each particle also remembers their best result as local best and the overall best of the entire population as global best. At each step, a particle has to move to a new position by adjusting its velocity so that its moves towards its local best (Pi,pbest) and the overall global best (Pi,gbest) record by the particle in the population. This is done iteratively until the optimal solution is achieved. Each and every particle in the population represents the following set of parameter, ⟨Position P t i,l |V elocity V t i,l⟩. Here, the position vector Pt i,l represents the lth link angles (θ1, . . . θl, . . . θnL) of ith particle in tth iteration. This is the required solution (Pi,l). The velocity Vt i,l represents an incremental change in the lth link angle of ith particle in tth iteration. With nP as total number of particles and nL as total number of links, i= 1,2, . . . , nP and l= 1,2, . . . , nL. Vt i,l =KV F Vt−1 i,l +Cpr1Pi,pbest −Pt−1 i,l  +Cgr2Pgbest −Pt−1 i,l (3) Pt i,l =Pt−1 i,l +Vt i,l (4) 43 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX The above equations represent the position and velocity of each particle in the population in which KV F represents the current velocity factor, r1&r2are uniformly distributed random number within the range [0,1], Cp&Cgrepresent the learning rates of local best and global best particles respectively. The values for these parameters are listed in are listed in Table 2. 3.2 Bat Algorithm (BAT) developed by Xin-She Yang and Amir Hossein Gandomi [47]. It was inspired from echolocation behavior of bats with the varying pulse rate of emission and loudness which assists them in seeking for prey and/or avoids obstacles in complete darkness. All bats use echolocation to sense the distance. They fly randomly with a certain velocity and with a fixed frequency. During their flight, bats emit a sound pulse with particular loudness and listens to its echo that bounces back from the surrounding environment. Based on the difference in time between the emission of sound pulse and echo, they detect the distance of the prey, its orientation and even its motion. Based on their location the bats adjust its velocity and position thereby reaching the prey. The same principle is used by them to detect an obstacle and avoid them. Each and every bat in the population represents the following set of parameters: ⟨Position P t i,l, V elocityV t i.l, frequency fi, Loudness At i, Pulse rate rt i⟩. Here, the position vector Pt i,l represents the lth link angles (θ1, . . . θl, . . . , θnL) of ith bat in tth iteration. The velocity Vt i,l represents an incremental change in the lth link angle of ith bat in tth iteration. Bats emits wave with frequency in the range [fmin fmax]. At i&rt irepresents the loudness and pulse rate of the wave emitted by the ith bat in tth iteration. Of all these parameters, the position vector Pt i,l i.e., the link angles of the robot manipulator (θ1, . . . θl, . . . , θnL) represents the solution. Here, i= 1,2, . . . , nB and l= 1,2, . . . , nL the velocity of each bat and their position are calculated based on the following equations. Vt i,l =Vt−1 i.l +Pt−1 i.l −Pgbestfi(5) Pt i,l =Pt−1 i,l +Vt i,l (6) To prevent the algorithm from getting stuck in a local minima / maximum, and to increase the exploration capability, a random walk is performed. Based on the pulse rate rt i, few bats are selected randomly to perform a random walk. A random position for the bat Pnew is generated using equation (7). This newly generated position is selected based on the fitness (RMS Error) & the loudness At iof the corresponding bat. The pulse rate rt iand loudness At iof each bat is updated using equations (8,9) only if this new bat is accepted. Pnew =Pgbest +εAt(7) At+1 i=ωAt i(8) rt+1 i=r0 i[1 −exp(−γt)] (9) Here, εis a random number in the range [−1,1], Atis the average loudness of all the bats at time t,r0 iis the initial pulse rate. Once a bat has found its prey, the pulse rate increases and the loudness decreases. This is achieved using the parameters ω, the pulse frequency increasing coefficient and γ, the pulse amplitude attenuation coefficient. The values of these parameters are listed in are listed in Table 2. 3.3 Gravitational Search Algorithm (GSA) GSA is an optimization method based on Newtonian gravity and the laws of motion. This algorithm was developed by Rashedi et al. [37]. In the GSA, each particle in the search space is considered as a mass. Therefore, the GSA may be expressed as an artificial mass system. All masses in the search space attract each other according to Newton’s gravity law and interact to exert force on each other with the force of gravity. Acting in the search space, masses exposed to these forces achieve the optimal solution. Each and every particle mass in the system is represented by its position Pt i,l. The gravitational constant is initialized at the beginning and will be reduced with time as in equation (10) to control the search accuracy. Gt=G0e(−αt/Imax)(10) Here, Gtis the gravitational constant at time any tand G0is the initial gravitational constant, αis the constant exponent factor. Here, t= 1,2, . . . , Imax, where Imax is the maximum iteration count of the algorithm. The fitness (RMS Error) of each mass Pt i,lare calculated and from this the best fitness value, i.e., the minimum RMS error (errt min) and worst fitness value, i.e., the maximum RMS error (errt max) are selected. From these, the gravitational mass ( mt i) of the ith particle in the system at time tis calculated using equations (11) and the normalized gravitational mass Mt iis calculated using equation (12). mt i=errt i−errt max errt max −errt min (11) Mt i=mt i PnP i=1 mt i (12) Here, nP is the total number of particle mass in the system. Now the force between the ith &kth particle mass in the system Ft i,kand the total force Ft i acting on any ith particle mass at any time tare calculated using equation (13) and (14) respectively. Ft i,k =Gt Mt i·Mt k eudt i,k +eps!Pt k,l −Pt i,l(13) Ft i= nP X k∈bestK,k=i rand ·Ft i,k (14) Here, Pt i,l &Pt k,l are the position of ith &kth particle mass and eudt i,k is the Euclidian distance between 44 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX Kanagaraj et al.: SettingsMeta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking ... Table 2: Parameters used for BAT, GSA, PSO and WOA. Algorithm Parameters Value BAT Frequency Range [fmin, fmax] [−0.05,0.05] Pulse frequency increasing coefficient (ω) 0.97 Pulse amplitude attenuation coefficient (γ) 0.1 Loudness (At i) Random number in the range [1 2] Initial Pulse Rate (rt i) Random number in the range [0 1] GSA Initial gravitational constant G0100 Constant exponent factor (α) 1 Random best agents (bestK) (100 to 2)% PSO Current Velocity Factor (KV F ) 0.3033 Learning rates (Cp&Cg) 2.4 & 3.2 Uniformly distributed random number (r1&r2) Random number in the range [0 1] WOA Shape of spiral (b) 5 Solution selection parameter (p) Random number in the range [0 1] them which is calculated using equation (15). eps is a constant value added to avoid divide by zero condition, rand is a random number between 0 and 1, bestK is the set of for kagents with best fitness values, i.e., minimum RMS errors. eudt i,k =||Pt k,l −Pi, lt|| =v u u t nL X l=1 (θt k,l −θt i,l)2(15) Now the new solution or the update in position of the masses in the system are updated using the following equations. Pt i,l =Pt−1 i,l +Vt i,l (16) Vt i,l =rand ·Vt i,l +At i,l (17) At i,l =Ft i Mt i (18) Here, Vt i,l &At i,l are the velocity and acceleration of the ith agent at time t. 3.4 Whale Optimization Algorithm (WOA) Mirjalili and Lewis [32] developed WOA which mimics the intelligent hunting behavior of humpback. This foraging behavior is called bubble-net feeding method which is observed only in humpback whales. The humpback whales dive down approximation 12 m and then create the bubble in a spiral shape around the prey. They then swim upward the surface following the bubbles and hunt the prey. Humpback whales can find the place of prey and encircle them. The WOA algorithm considers current best search agent position be the target prey or close to the optimum point, and other search agents will try to update their position towards the best search agent. Applications of WOA were presented in this work [20]. Each and every agent in the population is represented by its position vector Pt i,l. Here, the position vector represents the lth link angles (θ1, . . . θl, . . . θnL) of ith agent in tth iteration. This is the solution vector. Here, i= 1,2, . . . , nA and l= 1,2, . . . , nL.nA is the total number of agents in the population and nL is the dimension of the solution, i.e., it represents the number of links of the robot manipulator in this work. The position of each agent is calculated and updated based on the following equations. Pt i,l =nPt gbest −(A∗dist) if p < 0.5 dist′∗eb∗lrand ∗cos (2π∗lrand) + Pt gbest if p≥0.5 (19) dist = 2 ·rand ·Pt gbest −Pt i,l (20) dist′=Pt gbest −Pt i,l(21) Here, Pt i,l is the position of the ith agent with ldimension at time t,Pt gbest is the best agent so far at time t,rand &pare random numbers in the range [0 1]. A and lrand are calculated as follows. A= 2 ·rand ·a1−a1 (22) lrand = 1 + (a2−1) ·rand (23) Here, a1 is a linearly decreasing vector from 2 to 0 over the course of iterations and a2 is a linearly decreasing vector from -1 to -2 over the course of iterations, b represents the shape of the spiral. 3.5 Algorithm Parameters The values of the optimization parameters common to the algorithms are set as follows; Population size n= 60; the maximum iteration count used for trajectory following testing Imax = 100; the minimum acceptable error used for convergence testing is 0.1×10−4. The parameters related to BAT, GSA, PSO and WOA used in this work are listed in Table 2. 3.6 Framework The steps involved in carrying out the simulation are as follows. The 3D model of YASHKAWA MH5 robot is imported into MATLAB simulation environment [31]. 1. Generate the trajectory for the selected curve (Linear / Curvilinear / Saw Tooth / Rose Curve) 45 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX Figure 1: Flowchart of the Simulation Framework. 2. Now each and every point (xd, yd, zd)jin the generated trajectory, including waypoints, are given as input to the algorithms (BAT, GSA, PSO & WOA). Here j= 1,2, . . . , nPT and nP T is the number of points in the generated trajectory. The output from the algorithm are the link angles (θ1, . . . , θl, . . . , θnL)jto reach the jth point (xd, yd, zd)j.Here nL represents the number of links in the robot manipulator. 3. Now, using the link angles (θ1, . . . , θl, . . . , θnL)j, generated by the algorithm, the actual point reached by the robot manipulator (xa, ya, za)jis calculated. 4. The trajectory obtained using the points (xa, ya, za)j(obtained from BAT / GSA / PSO / WOA) is plotted and compared with the expected trajectory using the points (xd, yd, zd)j.and the error between these trajectories are generated. 5. The results are analyzed, interpreted and final conclusion is arrived. The entire framework of the simulation carried out in this work presented in the flowchart in Figure 1. 3.7 Solution Representation and Fitness Evaluation This section explains how the solution, i.e., link angles (θ1, . . . , θl, . . . θnL) generated by the algorithm are transformed into a cartesian point (xa, ya, za) in space, and based this how the fitness is evaluated. If we apply the link angles (θ1, . . . , θl, . . . θnL) to a robot manipulator, it will take a configuration and its end effector will reach a particular point in space. This end effector position calculation based on the link angles (θ1, . . . θl, . . . , θnL) is called forward kinematics. This can be calculated using the transformation matrix [12]. The general form of a standard transformation matrix from link l−1 to link lis given by link l−1 linkT=     cosθl−sinθl0al−1 sinθlcosαl−1cosθlcosαl−1−sinαl−1−sinαl−1dl sinθlsinαl−1cosθlsinαl−1cosαl−1cosαl−1dl 0 0 0 1     (24) Where, al−1is the link length, αl−1is the twist angle, dlis the joint offset, θlis the joint angle or the link angle of the corresponding link. Here l= 1,2, . . . , nL and nL represents the number of links in the robot manipulator. The transformation matrix in equation (24) is calculated for base to link 1 base link 1T, link 1 to link 2 link 1 link 2T, till link nL −1 to the end effector link nL−1 end effectorT. Multiplying all these matrices will give the transformation matrix from base to the end effector of the robot manipulator. base end effectorT=base link 1T·link 1 link 2T· · · link nL−1 end effectorT(25) In resulting 4 ×4 matrix base end effectorTfrom equation (24), the values in the location (1,4) ,(2,4) ,(3,4) give the (x, y, z) coordinate of the end effector of robot manipulator. Thus, the link angles (θ1, . . . θl, . . . , θnL) are transformed into coordinates (xa, ya, za). The input to this algorithm is any jth point (xd, yd, zd)jfrom the trajectory. This is the desired coordinate that the robot manipulator has to reach. From the randomly generated link angles (θ1, . . . θl, . . . , θnL), the coordinate point actually reached by the robot manipulator (xa, ya, za) is calculated. The fitness of this solution is calculated based on its deviation from the desired position (xd, yd, zd)j. This deviation is calculated based on the Root Mean Square (RMS) error, i.e., the fitness of the solution for any jth point in the trajectory is calculated as RMSerrorj=q(xa−xd)j 2+ (ya−yd)j 2+ (za−zd)j 2 (26) 4 Results and Discussions The simulation environment and the test cases used for measuring the efficiency of the proposed methodology are explained in this section. MATLAB (R2020b) environment is used to simulate the proposed work. YASKAWA’s MH5 robot manipulator is considered to carryout the stated objective. The actual robot and the simulated robot are shown in Figure 2. 46 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX Kanagaraj et al.: SettingsMeta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking ... (a) (b) Figure 2: (a) YASKAWA’s MH5 robot manipulator Physical robot (b) Robot simulated in MATLAB environment. Table 3: YASKAWA’s MH5 Robot Manipulator Parameters [33]. Name of Link Angle Maximum the Line Range (◦)Speed (◦/s) S-Axis: Swivel Base −170 ≤θ1≤+170 376 L-Axis: Lower Arm −065 ≤θ2≤+150 350 U-Axis: Upper Arm −136 ≤θ3≤+255 400 R-Axis: Arm Roll −190 ≤θ4≤+190 450 B-Axis: Wrist Bend −135 ≤θ5≤+135 450 T-Axis: Tool Flange −360 ≤θ6≤+360 720 This robot has 6 links (nL). Therefore the six link angles (θ1, θ2. . . , θ6) corresponding to the 6 links are adjusted simultaneously to control the position and orientation of the robot’s end effector. The other parameters related to this robot manipulator are shown in Table 3. The performance of the proposed algorithms for trajectory following are assessed using four different trajectories, viz., linear, curvilinear, sawtooth and rose trajectories. The trajectories simulated in MATLAB environment are shown in Figure 3. Each trajectory consists of j points, where, j= 1,2, . . . , nPT, and nP T is the number of points in the trajectory. In this paper, nP T is considered as 17, i.e., we generate 15 intermediate points between the start and endpoint of the trajectory, plus one start point and one end point. The number of intermediate points in the trajectory is selected as 15 for the sake of simplicity and ease of comparison. The number of points in the trajectory can be decided based on the application. The error between the expected trajectory and the trajectory generated by BAT, GSA, PSO and WOA are calculated for these 17 points. The objective of this work is to find the joint angles of the robot manipulator to reach the points in the trajectory with minimum error. So, the proposed work emphasis on following the trajectory with minimum error and does not depend on the number of points in the trajectory. The algorithm is tested for its closeness in following the trajectory. All experiments are deployed on the PC with Intel (R) Core (TM) i3, M370 @2.40GHz, 8GB RAM using MATLAB 2020b. (a) (b) (c) (d) Figure 3: Expected reference trajectory (a) linear (b) curvilinear (c) sawtooth (d) rose curve. 47 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX (a) (b) (c) (d) Figure 4: Average error plots (a) linear trajectory (b) curvilinear trajectory (c) sawtooth trajectory (d) rose curve trajectory. Table 4: Average Error ×10−4between expected and generated trajectory. Algorithm Trajectory Intermediate points in the Trajectory Type Start 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 End Linear BAT 0.02 0.02 0.03 0.03 0.02 0.05 0.02 0.03 0.03 0.07 0.03 0.02 0.06 0.02 0.03 0.04 0.06 GSA 0.20 0.23 0.22 0.29 0.23 0.22 0.21 0.20 0.16 0.20 0.24 0.24 0.22 0.31 0.34 0.33 0.20 PSO 3.00 2.95 2.48 1.88 2.43 1.03 1.79 3.43 3.45 1.82 1.79 2.25 2.99 3.03 1.11 2.97 3.79 WOA 6.50 6.30 5.21 6.09 3.68 2.76 4.17 2.85 3.65 0.60 0.62 1.23 0.79 2.60 2.35 2.62 4.11 Curvilinear BAT 0.03 0.03 0.02 0.03 0.02 0.03 0.02 0.03 0.02 0.03 0.03 0.03 0.02 0.02 0.02 0.03 0.04 GSA 0.18 0.18 0.19 0.21 0.25 0.29 0.24 0.28 0.21 0.33 0.28 0.35 0.31 0.36 0.33 0.25 0.31 PSO 1.50 1.56 0.78 0.77 1.00 0.78 0.72 0.65 0.66 0.91 0.59 0.58 0.56 0.54 0.44 0.70 1.12 WOA 5.21 4.36 2.15 1.35 1.22 1.63 0.89 0.58 0.77 0.31 0.50 0.55 0.39 0.92 1.92 1.67 2.83 Sawtooth BAT 0.06 0.04 0.04 0.04 0.09 0.02 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.02 GSA 0.19 0.20 0.21 0.22 0.21 0.17 0.18 0.26 0.21 0.26 0.35 0.28 0.24 0.21 0.18 0.20 0.22 PSO 1.66 1.60 1.14 1.57 1.21 1.16 0.75 1.22 0.96 0.65 0.89 1.99 1.49 1.64 1.07 1.28 1.48 WOA 0.83 0.53 1.07 4.49 4.87 4.62 3.49 2.35 2.16 3.52 4.89 5.06 3.90 4.15 4.53 4.66 4.64 Rose curve BAT 0.04 0.03 0.02 0.03 0.03 0.02 0.03 0.02 0.03 0.03 0.02 0.03 0.03 0.03 0.02 0.03 0.03 GSA 0.18 0.21 0.21 0.27 0.19 0.22 0.19 0.14 0.15 0.14 0.18 0.17 0.16 0.17 0.20 0.17 0.17 PSO 2.36 1.38 2.01 1.01 1.67 2.33 1.24 1.91 1.50 0.74 1.40 1.38 1.17 1.35 1.01 1.08 1.47 WOA 6.02 6.02 2.26 0.43 0.70 3.16 2.63 4.86 2.51 4.11 6.70 7.09 4.95 1.92 1.97 6.81 5.34 Table 5: Statistical analysis of average error ×10−4, bold face represents best values. Algorithm Trajectory Type Mean Median Standard Deviation Minimum Maximum Linear BAT0.03 0.03 0.01 0.02 0.07 GSA 0.23 0.22 0.05 0.16 0.34 PSO 2.33 2.48 0.79 1.03 3.79 WOA 2.63 2.85 1.88 0.60 6.50 Curvilinear BAT0.03 0.03 0.00 0.02 0.04 GSA 0.26 0.28 0.06 0.18 0.36 PSO 0.77 0.72 0.31 0.44 1.56 WOA 1.17 1.22 1.35 0.31 5.21 Sawtooth BAT 0.03 0.03 0.02 0.02 0.09 GSA 0.22 0.21 0.04 0.17 0.35 PSO 1.23 1.22 0.35 0.65 1.99 WOA 2.98 4.15 1.49 0.53 5.06 Rose curve BAT 0.03 0.03 0.00 0.02 0.04 GSA 0.18 0.18 0.03 0.14 0.27 PSO 1.41 1.38 0.44 0.74 2.36 WOA 3.13 4.11 2.02 0.43 6.81 48 MENDEL — Soft Computing Journal, Volume 28, No. 1, June 2022, Brno, Czech RepublicX Kanagaraj et al.: SettingsMeta-Heuristics Based Inverse Kinematics of Robot Manipulator’s Path Tracking ... (a) (b) (c) (d) Figure 5: Statistical plot of average error for different trajectories (a) Linear (b) Curvilinear (c) Sawtooth (d) Rose curve. 4.1 Trajectory Following Testing Each of the 4 different trajectories are tested with four different algorithms. So, a total of 16 different test cases are carried out for trajectory following experiment. To investigate the robustness, each of these 16 different cases are run for 20 times and the average objective values are recorded. The performance of the algorithm is measured in terms of its closeness in following the trajectory. For each and every point in the reference trajectory as input, the algorithms are executed for a maximum iteration count Imax = 100. At the end of 100th iteration, the best solution is taken as the output. The trajectory is generated for all the four curves considered in this paper using the proposed algorithms. Each of this case is run for 20 times and the average value of error between the expected and the generated trajectories are recorded. The plot of average error obtained for these algorithms are shown in Figure 4 and the values are listed in Table 4. Average error values in Table 4 shows that BAT and GSA closely follows the reference trajectory with minimum error when compared to other algorithms considered in this paper. Both BAT and GSA produced consistent result at all times. This is evident from the plot as there are very less variation in the average error. This is further investigated using the statistical analysis which are listed in Table 5. Table 5 shows that GSA’s performance is next best. Both BAT and GSA were able to provide solution, with consistent error values, for all points in the trajectory. Where as PSO and WOA have a higher standard deviation value. This is evident from the box whisker plot shown in Figure 5. 4.2 Convergence Testing The objective of every optimization algorithm is to converge to a global optimal solution. This convergence capability of BAT, GSA, PSO and WOA for the inverse kinematics problem is tested as follows. In trajectory following testing, the stopping condition was the maximum iteration count Imax = 100, i.e., the algorithm is executed for 100 iterations and at the end of 100th iteration, the best solution is taken as the output. In convergence testing, a minimum acceptable error is set as the stopping condition. The algorithm is executed till the minimum acceptable error is achieved. So, a random point in the trajectory is selected and given as input to the algorithm. The minimum acceptable error is set as 0.1×10−4. To prevent the algorithm from running continuously, the iteration limit is set as 1000. Each algorithm is executed 20 times and the average value of loop count and execution time required to achieve this acceptable error are recorded and listed in Table 6. 49