Inverse Kinematic Solver Based on Bat Algorithm for Robotic Arm Path Planning

Abstract

The bat algorithm (BA) is a nature inspired algorithm which is mimicking the bio-sensing characteristics of bats, known as echolocation. This paper suggests a Bat-based meta-heuristic for the inverse kinematics problem of a robotic arm. An intrinsically modified BA is proposed to find an inverse kinematics (IK) solution, respecting a minimum variation of the joints’ elongation from the initial configuration of the robot manipulator to the proposed new pause position. The proposed method is called IK-BA, it stands for a specific bat algorithm dedicated to robotic-arms’ inverse geometric solution, and where the elongation control mechanism is embedded in bat agents update equations. Performances analysis and comparatives to related state of art meta-heuristics solvers showed the effectiveness of the proposed IK bat solver for single point IK planning as well as for geometric path planning, which may have several industrial applications. IK-BA was also applied to a real robotic arm with a spherical wrist as a proof of concept and pertinence of the proposed approach.

1. Introduction

Robotic manipulators are characterized by their versatility due to means of their mechanical structure that consists of serial links joined together to be like a human arm, where the motion of their end-effector is dependent on joint movements [1]. This grants them the adaptability to perform in multiple processes, from simple applications such as pick and place operations, to more complex processes that involve applications such as welding, assembly, and painting [2]. These tasks typically require adequate paths to be tracked by the robot’s end-effector [3,4]. These paths are generally described in the Cartesian space, while the end-effector position is set according to the robot joint angles. Therefore, the key issue in path tracking is the mapping between the Cartesian space of the tracked path and the joint space of the robot manipulator. Such mapping is achieved through a kinematic modeling of the robotic arm, where a forward kinematics (FK) model allows determining of the end-effector location from a specific configuration of the robot joints. The FK problem possesses a unique solution, considering that each position of the robot joints drives to a one single location of the end-effector. A common way of solving the FK problem is to use the Denavit-Hartenberg notation [5]; whereas, the inverse kinematics (IK) intends the opposite, i.e., determine the corresponding joint variables that enables the end-effector to reach a specific target position [6]. Figure 1 interprets the relation between the forward and inverse kinematics.

Unlike FK, the IK problem often has many different solutions, since one target point could be reached through multiple positions of the manipulator. Hence, solving the IK problem for robotic manipulators is a challenging task. The geometry of the robot and the nonlinear trigonometric equations that define the mapping between Cartesian space and joint space contribute to the problem’s complexity [7,8]. This problem has been approached by several conventional methods with deterministic characteristics, where their final solution is determined by the initial conditions of the problem [9]. For robots with high degrees of freedom (DOF), the IK problem is more sophisticated, hence, conventional methods require more computational processes, which slow down controlling the manipulator in real time [6,10].

Conventional methods are generally divided into algebraic, geometric, and numerical methods. Algebraic and geometric methods aim to find closed-form solutions to the problem. While they are indeed faster than numerical methods, their provided solutions are not generic but rather robot-dependent. Meanwhile the numerical methods are applicable to any kinematic structure, therefore they are not robot-dependent; however, they are generally time consuming and unable to calculate all possible solutions [11].

With the expansive usage of robot manipulators in production lines, and the prerequisite for higher productivity with high accuracy, many researchers have directed their efforts in searching for alternative approaches that tackle the IK problem while overcoming the disadvantages of conventional methods. The application of artificial neural networks (ANN) is an example of advanced investigations [12,13]. Cursi et al. [14] proposed the use of a type of ANN known as a probabilistic neural network to enhance the task accuracy of a surgical macro–micro robot manipulator. However, using ANNs suffers from long training time and lack of accuracy [15], which has motivated further investigations using other soft-computing alternatives, such as meta-heuristic (MH) algorithms. MH algorithms are characterized by non-deterministic behavior, which eases the obtaining of approximate satisfying solutions of an optimization problem with a plausible computational cost [16,17]. They are exploited to tackle the IK problem using only the FK model. Applying MH algorithms requires modeling the IK problem as an optimization problem [18], where an objective function is defined according to the manipulator’s FK model and the desired pose of the robot’s end-effector. The MH algorithm is dedicated to find the corresponding joints values that enable the robot manipulator to achieve a desired pose. In a case where the objective function is intractable or impossible to be formulated analytically, it is preferable to employ a black-box optimization approach. In this approach, a close form of the objective function is not required [19]. For instance, Cursi et al. [20] developed an open source global optimization framework for robot design. Considering that the robot design is regarded as an optimization problem with a black-box cost function; therefore, their developed framework implements a Bayesian optimization method to approximate an appropriate formulation of the objective function.

This paper suggests a sensible adjustment of a MH algorithm, known as the bat algorithm (BA) for the IK problem, IK-BA, of robot manipulators with spherical wrists, which is typical for industrial robotic arms. The main motivation of the use of the MH approach for open chain robotic arm is to avoid the singularities of the classical invers solvers. Singularities consist in configurations where the inverse kinematic equations of the robotic arm do not have a unique solution [21]. Classical IK solvers may not provide a meaningful solution at singular configurations, since the Jacobian matrix becomes singular, which means that it cannot be inverted to calculate the required joint velocities [22]. By using MH optimization approaches as IK solvers, the robotic arm can avoid these singularities by exploring alternative configurations and selecting the optimal solution based on a predefined fitness function. This approach can help to improve the overall performance and reliability of the robotic arm, allowing it to carry out complex tasks with greater accuracy and efficiency [23].

The proposed IK-BA is an intrinsically hybridization of BA, which aims to accord an IK solution, while assuring a minimal variation from the initial configuration of the robot manipulator until the solution’s position. The elongation control mechanism is embedded within BA update equations, making the proposal specific to IK. This contraction between the initial configuration and the generated solution, guarantees a shorter path compared to another IK solution provided without considering the initial configuration. The adjustment of the BA algorithm lies in modifying its update equations in a way to generate solutions affected by the initial configuration. The remaining of the paper is organized as follows: paragraph 2 details the key related works and similar proposals; paragraph 3 stands for the IK problem statements with elongation control. In paragraph 4, the BA algorithm is presented before introducing the IK-BA proposal, it ends with a flowchart describing how a path planning can be conducted using the suggested method. An experimental investigation is presented in paragraph 5, including a comparative analysis and a time assessment test, based on an industrial robotic arm. Conclusion and discussions are detailed in paragraph 6.

2. Related Works

Early applications of meta-heuristic algorithms for solving the IK problem for redundant robotic manipulators have emanated from using one of the initial paradigms of MH algorithms, which is the genetic algorithm (GA). The GA is an MH developed by Holland [24]. It is inspired by the process of natural selection and evolution, which is used to find approximate solutions to optimization and search problems. Genetic algorithms are commonly used in a variety of fields, including machine learning, engineering, to solve complex optimization problems that are difficult to solve using traditional methods [25].

Parker et al. [26] tested the GA in solving the IK problem of a 4 DOF serial manipulator, where the solution provides a minimized largest joint displacement in a point-to-point positioning task. Furthermore, Nearchou [27] suggested using a modified version of GA to solve the IK problem as a way to minimize both the end-effector’s positional error and the robot’s joint movements, while keeping the solution free from collisions with obstacles in the robot’s workspace. Both contributions used the weighted sum approach to combine the minimization of the positional error and the minimization of joint motions objectives, into one objective function. Tarokh and Zhang [28] used the same approach with an additional objective function dedicated to the orientation employing an adaptive GA to tackle the trajectory tracking problem of robot manipulator. In the weighted sum approach, a multi-objective optimization problem is transformed into a single objective, where each objective function is multiplied by a user-supplied weight. The weight value assigned to an objective function is generally proportional to the objective’s relative importance. However, the importance of an objective is a vague concept, hence, the choice of these coefficients might be inaccurate, whereas the final solution is strongly affected by the chosen weight coefficients, which eventually might lead to an inaccurate solution [29].

Another paradigm of MH algorithms is the particle swarm optimization (PSO) algorithm. It was developed in 1995 by Kennedy and Eberhart [30]. The PSO is inspired by the social behavior of birds’ flocks and fish schools. It has been applied as an IK solver, as in [31,32,33].

Adly and Abd-El-Hafiz [34] have intended to solve the IK problem with minimizing the range difference between the initial joints position and the IK solution values that represent the final position. They suggested a multi-objective optimization approach, where the first objective function is dedicated for the IK solution; meanwhile, the second function formulates the variation from the initial joint values and the solution values. They applied an extension of PSO for multi-objective optimization problems called the Time Variant Multi-Objective Particle Swarm Optimization (TV-MOPSO) algorithm [35]. This approach provides an IK solution that guaranties an optimal trajectory of the end-effector to move from one point to another. However, compared to the single objective approach, the multi-objective method is generally more time-consuming as it requires additional processes for generating a set of trade-off solutions that satisfies two or more conflicting objectives.

Rokbani et al. [36] proposed a multi-objective version of a new modified particle swarm optimization MO-m-PSO to solve the inverse kinematics of a 5-DOF robotic arm. Its IK problem is formulated as a multi-objective problem that takes into consideration the end-effector pose (position and orientation).

Various other swarm-based meta-heuristics were employed for the IK issue, such as Cuckoo Search Algorithm [37], Artificial Bee Colony [38,39], Firefly Algorithm [40,41], and Salp Swarm Algorithm [42]. Most investigations about the application of MH algorithms to the IK problem proved that MH algorithms are suitable and reliable IK solvers. Accordingly, some researchers have made further analyses attempting to employ MH algorithms for robot manipulator trajectory tracking. Lopez-Franco et al. [43] presented a comparative study on six soft-computing algorithms in solving the IK and path tracking problems, where the best algorithm according to the comparative simulations is applied on a real robot manipulator. Kanagaraj et al. [44] tested the performance of four MHs in path tracking simulations, where the algorithms are required to solve the IK problem for every point of a generated trajectory. The tested algorithms are particle swarm optimization (PSO), whale optimization algorithm (WOA), gravitational search algorithm (GSA), and the bat algorithm (BA). The simulation results confirmed the outperformance of the BA compared to the others.

When robots need to follow a specific orientation with a minimal angles variation during point-to-point tasks, the two most common approaches used are the classical weighted sum technique or multi-objective optimization methods. The weighted sum technique simplifies the problem, but may not be so accurate, while the multi-objective optimization requires more computing time. In this paper, the IK-BA is suggested to solve the IK problem assuring a minimal variation of joint angles (elongation). The proposed method takes into account the initial joint configuration that restricts the exploration of the algorithm, hence guaranteeing a limited variation of angles from the initial configuration towards the final solution. It is tested on solving the IK problem of KUKA LBR iiwa 14 R820 robotic arm; it is compared with ordinary BA, and other MH-IK solvers from the literature. The comparatives are in terms of accuracy and angles variation of their provided IK solutions.

3. Inverse Kinematics Problem

The positional accuracy of a robotic manipulator is assessed by the distance between the desired position to reach and the actual position of the robot’s end-effector, which is mathematically calculated using the following equation:

$${\epsilon }_p=\sqrt{{\left(x-x_d\right)}^2+{\left(y-y_d\right)}^2+{\left(z-z_d\right)}^2}$$

(1)

where $\left(x,y,z\right), \left(x_d, y_d, z_d\right)$ denotes the end-effector’s coordinates and the desired coordinates, respectively. The lower the distance error $\epsilon$ is, the higher the accuracy.

The orientation of the end-effector could be described by a $3\times 3$ rotation matrix as represented in Equation (2) below:

$$R=\left[\begin{array}{ccc}{n}_x& o_x& a_x\\ n_y& o_y& a_y\\ n_z& o_z& a_z\end{array}\right]$$

(2)

where $\mathit{n}={\left[n_x{n}_y{n}_z\right]}^T$, $\mathit{o}={\left[o_x{o}_y{o}_z\right]}^T$ and $\mathit{a}={\left[a_x{a}_y{a}_z\right]}^T$ are the expression of $\overrightarrow{x}$, $\overrightarrow{y}$ and $\overrightarrow{z}$ axes, associated to the end-effector frame, relative to the base frame.

Let $R_d$ be required orientation of the end-effector. The orientation error ${\epsilon }_r$ could be defined as in [45] using Equation (3), where $\Vert .\Vert$ indicates the Euclidean norm of a matrix.

$${\epsilon }_r=\frac{\Vert R-R_d\Vert }{\Vert R_d\Vert }$$

(3)

The IK problem aims to determine the appropriate joint angles that minimize the position error, which could be formulated as an optimization problem as denoted in Equation (4):

$$\mathrm{Find}\mathsf{\Theta } =\left[{\theta }_1,{\theta }_2,\dots ,{\theta }_n\right]\mathrm{that}:\{\begin{array}{c}Minimize {\epsilon }_p and {\epsilon }_r\\ Subject to:\\ {\mathsf{\Theta }}_{\min }\le \mathsf{\Theta }\le {\mathsf{\Theta }}_{\max }\end{array}$$

(4)

where $\mathsf{\Theta }$ is a vector of joint angles of an n-dof robot manipulator. ${\mathsf{\Theta }}_{\min }$ and ${\mathsf{\Theta }}_{\max }$ are vectors that represent the upper bound and the lower bound of joint angles, respectively.

Formulating the IK problem as an optimization problem enables the employment of an optimization an algorithm to determine an adequate solution.

4. Proposed Method

4.1. Bat Algorithm

The bat algorithm (BA) is a nature-inspired algorithm, developed by Yang [46]. It is built around the bio-sonar characteristics of micro-bats, known as echolocation. Most bats use echolocation to navigate and catch prey in dark places. Micro-bats are able to emit ultrasonic pulses that bounce as an echo when clashing into a surrounding object. When the echo bounces to their ears, they could estimate the location of the object, based on the delay time between the emitted sound pulse and the echo. Each pulse has a constant frequency in the range of 25 KHz to 150 KHz, and they typically last 5 to 20 ms. Micro-bats emit around 10 to 20 sound pulses every second with a loudness in the range of 110 dB. While they became quieter when hovering close to a prey, the rate of pulses per second could increase to 200.

The BA is formulated based on the principle of echolocation characteristics of micro-bats, when it is associated with an objective function as a subject of optimization.

To numerically simulate the bats behaviour for an optimization problem, each bat is formulated with a set of parameters:

• Position $X$: corresponds to the position of the bat in the search space, representing a candidate solution of the optimization problem.
• Velocity $V$: is the velocity of the bat, which models an incremental variation of the position between 2 successive iterations.
• Frequency $f$: is the frequency of the emitted pulse. It is used to adjust the velocity change.
• Loudness $A$: represents the loudness of the emitted pulse. It enables a local search around the best solution.
• Pulse rate $r$: is the rate of the emitted pulses. It increases gradually along the iterations of the algorithm, assuming a consistent gradual compromise from an exploration phase to an exploitation phase.

The position $X_i$ and velocity $V_i$ of bat $i$ are updated in the search space at each iteration $t$, according to the following equations:

$$f_i=f_{min}+\left(f_{max}-f_{min}\right)\beta$$

(5)

$$V_i^t=V_i^{t-1}+\left(X_i^t-X_{best}\right)f_i$$

(6)

$$X_i^t=X_i^{t-1}+V_i^t$$

(7)

where $\beta \in \left[0,1\right]$ is a random generated vector drawn from a uniform distribution. $X_{best}$ is the current best position (solution) among the positions of all bats. Initially each bat is assigned by a random frequency $f_i$ uniformly drawn from $\left[f_{min},f_{max}\right]$ using Equation (5).

A local search is performed around the selected best solution, where a new solution is locally generated using random walk using Equation (8):

$$X_{new}=X_{gbest}+\epsilon A^t$$

(8)

where $\epsilon$ is random number in the range of [−1, 1] generated from a normal distribution. $A^t$ is the loudness at iteration $t$.

The loudness $A$ decreases progressively while the rate pulse increases along the iterations. This mimics the natural behaviour of bats when they get closer to their prey. This formulated in Equations (9) and (10):

$$A_i^{t+1}=\alpha A_i^t$$

(9)

$$r_i^{t+1}=r_i^0\left[1-exp\left(-\gamma t\right)\right]$$

(10)

where $\alpha$ and $\gamma$ are constants. $r_i^0$ is the initial rate of pulse emission that corresponds to the maximal value $r$ attained when $t\to \infty$. According to Equations (9) and (10), for any $0<\alpha <1$ and $\gamma >0$, the variation of the loudness $A_i$ and the pulse rate $r_i$, along the iterations proceeding is described as below:

$$A_i^t\to 0, r_i^t\to r_i^0, \mathit{as}t\to \infty$$

(11)

Algorithm 1 is a pseudo code that elucidates the approach of the BA.

Algorithm 1: Pseudo code of BA

Input: Fitness function: $fit(X={\left[x_1,\dots ,x_n\right]}^T)$, Number of bats: $N$, Maximum number of iterations: $t_{max}$, frequency range: $f_{min}$ and $f_{max}$, increasing coefficient: $\alpha$, Attenuation coefficient: $\gamma$

Output: Best solution $X_{best}$

Initialize the bat population $X_i\left(i=1, 2, \dots , N\right)$ and $V_i$;

Define pulse frequency $f_i$ at $x_i$;

Initialize pulse rates $r_i$ and the loudness $A_i$;

while $t<t_{max}$ do

for $i=1:N$ do

Adjust frequency using Equation (5);

Update velocity and position using Equations (6) and (7), respectively;

if $r_i<rand$ then // rand is a random number in [0, 1]

Select a solution among the best solutions $X_{best}$;

Generate a local solution around the selected best solution using Equation (8);

End

evaluate the new solution according to the objective;

if $rand<A_i \& f\left(X_i\right)<f\left(X_{best}\right)$ then

Accept the new solution;

Increase $r_i$ and decrease $A_i$ using Equations (9) and (10);

End

Rank the bats according to their fitness and update the best solution $X_{best}$;

End

End

4.2. Bat Algorithm for Inverse Kinematics

In this study, a simplified version of BA is employed, where the loudness $A$ and the pulse rate $r$ are equivalent for all bats as in [47]; in order to make it easier to implement, reduce the computational complexity for a faster convergence of the algorithm, and to provide a better balance between exploration and exploitation.

The proposed IK-BA ensures a minimum angular variation, which could be measured using Equation (12).

$$Variation\left({\mathsf{\Theta }}_s,{\mathsf{\Theta }}_f\right)=\left|\left|{\mathsf{\Theta }}_s-{\mathsf{\Theta }}_f\right|\right|=\sqrt{{\displaystyle \sum }_{i=1}^d{\left({\mathsf{\Theta }}_s-{\mathsf{\Theta }}_f\right)}^2}$$

(12)

where ${\mathsf{\Theta }}_s=\left[{\theta }_{s1}, {\theta }_{s2},\dots ,{\theta }_{sd}\right] \mathrm{and} {\mathsf{\Theta }}_f=\left[{\theta }_{f1}, {\theta }_{f2},\dots ,{\theta }_{fd}\right]$ represent the positions of the initial configuration and solution’s configuration, respectively. $d$ is the the number of DOF of a robot manipulator.

The velocity update equation is modified as in Equation (13), where the initial configuration is included in the equation to restrict the bats positions in the search space according to the initial position.

$$V_i^t=V_i^{t-1}+\left(X_i^t-\frac{\left(X_{best}+{\mathsf{\Theta }}_s\right)}{2}\right)f_i$$

(13)

Another modification concerns the local search generated by Equation (8). It is modified to take into account the initial configuration at the expense of the best solution in the last iterations as in Equation (14). Furthermore, the local search acquired by this equation is realized at each iteration; hence, the condition $(r_i<rand)$ in Algorithm 1 is removed in the IK-BA.

$$X_{new}=r^t{X}_{gbest}+\left(1-r^t\right){\mathsf{\Theta }}_s+\epsilon A^t$$

(14)

Note that it is preferable to use the initial configuration instead of random positions, as initial position for at least one bat, to assure a minimal angle variation.

For the IK-BA, the bats position vector $X=\left( {\left[x_1,\dots ,x_n\right]}^T\right)$ corresponds to the joint positions of the robot, where each bat’s position represents a potential solution to the IK problem. The fitness function $fit()$ denotes the difference between the desired end-effector position and the actual end-effector position, as in Equation (1). $fit()$ is used to evaluate the quality of the solution $X$ at each iteration of the algorithm. New potential IK solutions are generated iteratively by modifying the current solutions using the echolocation of the bat, which is numerically simulated by the algorithm’s intrinsic parameters $f, A$ and $r$. Algorithm 2 presents the methodology of the proposed IK-BA used in this study.

Algorithm 2: IK-BA pseudo code

Input: IK fitness function: $fit\left(X={\left[x_1,\dots ,x_n\right]}^T\right)$, Number of bats: $N,$ Maximum number of iterations: $t_{max}$, frequency range: $f_{min}$ and $f_{max}$, increasing coefficient: $\alpha$, Attenuation coefficient: $\gamma$, Initial configuration: ${\mathsf{\Theta }}_s=\left[{\theta }_{s1}, {\theta }_{s2},\dots ,{\theta }_{sd}\right]$, Target position: $P=\left[p_x p_y p_z\right]$

Output: Best solution $X_{best}$

Initialize the bat population $X_i\left(i=1, 2, \dots , N\right)$ and $V_i$;

Define pulse frequency $f_i$ at $x_i$;

Initialize pulse rates $r_i$ and the loudness $A_i$;

while $t<t_{max}$ do

Increase $r_i$ and decrease $A_i$ using Equations (9) and (10);

for $i=1:N$ do

Adjust frequency using Equation (5);

Update velocity depending on ${\mathsf{\Theta }}_s=\left[{\theta }_{s1}, {\theta }_{s2},\dots ,{\theta }_{sd}\right]$ using Equation (13) and the position and using Equation (7);

Select a solution among the best solutions $X_{best};$

Generate a local solution around the selected best solution depending on${\mathsf{\Theta }}_s=\left[{\theta }_{s1}, {\theta }_{s2},\dots ,{\theta }_{sd}\right]$ using Equation (14);

evaluate the new solution according to the fitness function fit(X) and target position $P=\left[p_x p_y p_z\right]$

if $rand<A_i \& f\left(X_i\right)<f\left(X_{best}\right)$ then

Accept the new solution;

end

Rank the bats according to their fitness and update the best solution $X_{best}$;

end

end

4.3. IK-BA for Decoupled Position-Orientation

For a robot manipulator of $n$-dof with spherical wrist, the axes of the last three rotational joints, dedicated to the orientation of the end-effector, intersect in one point. The manipulator can be accordingly decoupled into two parts. The first $\left(n-3\right)$ joints are responsible for the position and their values would be determined using the IK-BA, whereas the last three joints would be calculated according to the required orientation. Based on the formulation in [21], the origin of the $\left(n-2\right)$ frame $O_{\left(n-2\right)}$ (point of intersection of the three last joints, that assign the orientation) defines the wrist’s position of the arm in relation with the end-effector’s position as shown in Figure 2 and formulated in Equation (15):

$$P_w=P-d_n\left[R\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$$

(15)

where, $d_n$ is the distance between $O_{n-1}$ and $O_n$ along $Z_{n-1}.$ $\left[R\right]$ is a $3\times 3$ rotation matrix of the end-effector relative to the base (frame 0), where its third column indicates the direction of $z_n$ and $z_{n-1}$ axes with respect to the base (since $z_n$ and $z_{n-1}$ are always collinear). $P_w$ and $P$ are the vectors denoting the position of the wrist and the position of the end-effector, relative to the base frame.

To include a required orientation in addition to the position, the calculation of the three last joints value using the IK-BA is excluded. In this case, the algorithm attends to finding the values of ${\theta }_i, i=1,\dots ,n-3$, that enable the wrist to be positioned in $Pw$ position, calculated by Equation (15). The remaining angles ${\theta }_{n-2}$, ${\theta }_{n-1}$ and ${\theta }_n$ are calculated according to the orientation of the wrist based on the ${\theta }_i, i=1,\dots ,n-3$ values provided by the IK-BA, and the required orientation ${}^0{R}_n$, using Equations (16) and (17):

$${}^0{R}_{n-3}\left({\theta }_1, {\theta }_2, \dots , {\theta }_{n-3}\right)={}^{0}R_1\left({\theta }_1 \right)\times {}^1{R}_2\left({\theta }_2\right)\times \dots \times {}^{n-4}{R}_{n-3}\left({\theta }_{n-3}\right)$$

(16)

$${}^{n-3}{R}_n\left({\theta }_{n-2}, {\theta }_{n-1}, {\theta }_n\right)={}^{0}R_{n-3}^{\left(-1\right)}\times {}^{0}R_n$$

(17)

In a robotic arm with a spherical wrist, the orientation of the end-effector is set according to the positions of the last three joints, in which their values could be deduced from the ${}^{n-3}{R}_n$ orientation matrix. Thus, all the required values ${\theta }_i, i=1,\dots ,n$ to attain a specific pose ${}^0{T}_n$ are obtained in two phases. In a first phase, the values of ${\theta }_i, i=1,\dots ,n-3$ are obtained using the proposed IK-BA in order to enable the wrist of the arm to reach the position calculated using Equation (15). Furthermore, in a second phase, the required rotational matrix to reach the desired orientation is calculated using Equations (16) and (17), and the remaining three joints value ${\theta }_i, i=n-2,\dots ,n$ are deduced from the required rotational matrix.

Figure 3 is a flowchart that summarizes the entire methodology.

5. Simulation Results

The proposed IK-BA is dedicated to providing an IK solution dependent on the initial configuration. The advantage of this approach lies in solving the IK problem while minimizing the joints variation from one point to another. The proposed methodology highlights its importance in path tracking tasks for assuring a smooth variation of the joints values from point to point of the trajectory, hence a path tracking test is suggested to feature this advantage.

The performance of the IK-BA is compared with the ordinary BA, and other MH algorithms from the literature.

The kinematic structure of KUKA LBR iiwa 14 R820 robot manipulator is considered for the path tracking test.

5.1. Kinematic Modeling of the KUKA LBR Iiwa 14 R820

KUKA LBR iiwa 14 R820 is a 7-axes robot manipulator. Its kinematic structure provides it with a maximum versatility that makes it effective in collaborative tasks with human operators. Figure 4 elucidate the allocation of the KUKA LBR’s joints and their associated frames. The intrinsic parameters of the algorithms used in the simulation are indicated in

Table 1. Intrinsic parameters values for the IK-BA and the ordinary BA.

Algorithm	Parameters	Values
IK-BA	Frequency range $\left[f_{min},f_{max}\right]$	[0, 2]
	Increasing coefficient $\left(\alpha \right)$	0.97
	Attenuation coefficient ($\gamma$)	0.1
	Loudness for each bat ($A_i^t$)	1
	Initial pulse rate ($r_i^t$)	1
BA [48]	Frequency range $\left[f_{min}, f_{max}\right]$	[0, 2]
	Increasing coefficient $\left(\alpha \right)$	0.5
	Attenuation coefficient ($\gamma$)	0.5
	Loudness for each bat ($A_i^t$)	Random number in the range of [0, 1]
	Initial pulse rate ($r_i^t$)	0.001
DE [45]	Scaling factor ($F$)	0.6
	Cross-over factor ($C$)	0.9
PSO [43]	Inertia weight ($w$)	1.1312
	Individual confidence factor ($c_1$)	2.0149
	Swarm confidence factor ($c_2$)	0.53514
K-ABC [38]	No intrinsic parameters
MO-PSO [49]	Inertia weight ($w$)	0.5
	Individual confidence factor ($c_1$)	2
	Swarm confidence factor ($c_2$)	2
	Number of grids in each dimension ($N_{grid}$)	20
	Maximum velocity ($V_{max}$)	5
	Uniform mutation percentage ($u_{mut}$)	0.5

.

In the Denavit–Hartenberg notation, each link relating two neighbouring joints $j-1$ and $j$ is described by 4 parameters, whereas each parameter denotes a spatial transformation between the frames associated to each joint. These parameters are:

• ${\theta }_i$: indicates a rotational transformation from frame $j-1$ to frame $j$ around the $z$ axis, therefore ${\theta }_j$ is the angle between $x_{j-1}$ and $x_j$ through $z_{j-1}$.
• $d_j$: indicates a displacement from frame $j-1$ to frame $j$ along the $z$ axis, it is measured as the distance from $o_{j-1}$ to $x_j$ along the $z_{j-1}.$
• $a_j$: indicates a displacement from frame $j-1$ to frame $j$ along the new $x$ axis, therefore it represents the distance between $z_{j-1}$ and $z_j$ axes along $x_j$.
• ${\alpha }_j$: indicates a rotational movement from frame $j-1$ to the frame $j$ around the new $z$ axis, therefore ${\alpha }_j$ is the angle between from $z_{j-1}$ to $z_j$ axes about the $x_j$ axis.

The pose (orientation and position) of frame $j$ relative to $j-1$, could be defined by a $4\times 4$ homogeneous matrix, that corresponds to the product of the four transformations denoted by (${\theta }_i$, $d_j$, $a_j$ and ${\alpha }_j$) as implied in Equation (18):

$$\begin{array}{c}{}^{j-1}T_j=Ro{t}_{z,{\theta }_j}·Tran{s}_{z,d_j}·Tran{s}_{x,a_j}·Ro{t}_{x,{\alpha }_j}\\ =\left[\begin{array}{cccc}C\left({\theta }_j\right)& -S\left({\theta }_j\right)& 0& 0\\ S\left({\theta }_j\right)& C\left({\theta }_j\right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]·\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& d_j\\ 0& 0& 0& 1\end{array}\right]·\left[\begin{array}{cccc}1& 0& 0& a_i\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]·\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& C\left({\alpha }_j\right)& -S\left({\alpha }_j\right)& 0\\ 0& S\left({\alpha }_j\right)& C\left({\alpha }_j\right)& 0\\ 0& 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{cccc}C\left({\theta }_j\right)& -S\left({\theta }_j\right)C\left({\alpha }_j\right)& S\left({\theta }_j\right)S\left({\alpha }_j\right)& a_{i}C\left({\theta }_j\right)\\ S\left({\theta }_j\right)& C\left({\theta }_j\right)C\left({\alpha }_j\right)& -C\left({\theta }_j\right)C\left({\alpha }_j\right)& a_{i}S\left({\theta }_j\right)\\ 0& S\left({\alpha }_j\right)& C\left({\alpha }_j\right)& d_j\\ 0& 0& 0& 1\end{array}\right]\end{array}$$

(18)

Matrix ${}^{j-1}T_j$ represents the transformation matrix from frame $j-1$ to frame $j$. $C()$ and $S()$ stands for $cos()$ and $sin()$ functions respectively.

Table 2. DH parameters of KUKA LBR iiwa 18 R820 robot manipulator.

Link $\mathit{i}$	${\mathit{\theta }}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$	${\mathit{a}}_{\mathit{i}}$	${\mathit{\alpha }}_{\mathit{i}}$ [deg]	Range [deg]
1	${\theta }_1$	$d_1$	0	−90	[−170, 170]
2	${\theta }_2$	0	0	90	[−120, 120]
3	${\theta }_3$	$d_3$	0	90	[−170, 170]
4	${\theta }_4$	0	0	−90	[−120, 120]
5	${\theta }_5$	$d_5$	0	−90	[−170, 170]
6	${\theta }_6$	0	0	90	[−120, 120]
7	${\theta }_7$	$d_7$	0	0	[−175, 175]

represents the DH parameters of the KUKA LBR iiwa 14 R820 robot manipulator according to Figure 4.

Based on the generic homogeneous matrix of Equation (18) and the DH parameters in Table 2, the last column in the table indicates the limits of each joint. The transformation matrix between each two successive joints of KUKA LBR iiwa 18 R820 are denoted from Equations (19)–(25).

$${}^{0}T_1=\left[\begin{array}{cccc}C\left({\theta }_1\right)& 0& -S\left({\theta }_1\right)& 0\\ S\left({\theta }_1\right)& 0& C\left({\theta }_1\right)& 0\\ 0& -1& 0& d_1\\ 0& 0& 0& 1\end{array}\right]$$

(19)

$${}^{1}T_2=\left[\begin{array}{cccc}C\left({\theta }_2\right)& 0& S\left({\theta }_2\right)& 0\\ S\left({\theta }_2\right)& 0& -C\left({\theta }_2\right)& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(20)

$${}^{2}T_3=\left[\begin{array}{cccc}C\left({\theta }_3\right)& 0& S\left({\theta }_3\right)& 0\\ S\left({\theta }_3\right)& 0& -C\left({\theta }_3\right)& 0\\ 0& 1& 0& d_3\\ 0& 0& 0& 1\end{array}\right]$$

(21)

$${}^{3}T_4=\left[\begin{array}{cccc}C\left({\theta }_4\right)& 0& -S\left({\theta }_4\right)& 0\\ S\left({\theta }_4\right)& 0& C\left({\theta }_4\right)& 0\\ 0& -1& 0& d_4\\ 0& 0& 0& 1\end{array}\right]$$

(22)

$${}^{4}T_5=\left[\begin{array}{cccc}C\left({\theta }_5\right)& 0& -S\left({\theta }_5\right)& 0\\ S\left({\theta }_5\right)& 0& C\left({\theta }_5\right)& 0\\ 0& -1& 0& d_5\\ 0& 0& 0& 1\end{array}\right]$$

(23)

$${}^{5}T_6=\left[\begin{array}{cccc}C\left({\theta }_6\right)& 0& S\left({\theta }_6\right)& 0\\ S\left({\theta }_6\right)& 0& -C\left({\theta }_6\right)& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(24)

$${}^{6}T_7=\left[\begin{array}{cccc}C\left({\theta }_7\right)& -S\left({\theta }_7\right)& 0& 0\\ S\left({\theta }_7\right)& C\left({\theta }_7\right)& 0& 0\\ 0& 0& 1& d_7\\ 0& 0& 0& 1\end{array}\right]$$

(25)

The representation of the end-effector’s pose relative to the base frame is obtained through a successive multiplication of transformation matrices between two adjacent joint frames, as denoted in Equation (26).

$${}^{0}T_7={}^{0}T_1·{}^{1}T_2·{}^{2}T_3·{}^{3}T_4·{}^{4}T_5·{}^{5}T_6·{}^{6}T_7=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ 0& 0& 0& 1\end{array}\right]$$

(26)

where $\mathit{n}={\left[n_x n_y n_z\right]}^T$, $\mathit{o}={\left[o_x o_y o_z\right]}^T$ and $\mathit{a}={\left[a_x a_y a_z\right]}^T$ are the elements of rotational transformation matrix, that describe the orientation of the end-effector in reference to the base frame (frame 0). $p_x, p_y$ and $p_z$ that indicates the coordinates of the end-effector’s position along the $\overrightarrow{{\mathit{x}}_0}$, $\overrightarrow{{\mathit{y}}_0}$, $\overrightarrow{{\mathit{z}}_0}$ axes.

Each set of values of these angles ${\theta }_i, i=1,\dots ,7$ results in a unique pose defined by ${}^{0}T_7$. In the case of redundant robot manipulators where the number of degree of freedom is higher than the sufficient end-effector’s active joints for a certain task [50], a particular position of the end-effector might be reached through multiple different configurations of the manipulator. Unlike the forward kinematics problem, the inverse kinematics problem has infinite solutions.

The IK-BA is used in this study to determine the corresponding joint values ${\theta }_i, i=1, \dots , 7$ that enable the end-effector to reach a desired pose defined by ${}^{0}T_7$, along with granting a minimal possible joints displacement from an initial configuration to the target pose, using the proposed IK-BA.

5.2. Path Tracking Test

The path tracking test is proposed to evaluate the performance of the algorithm in solving the IK problem for each point along a trajectory. The efficiency of the IK-BA in terms of the accuracy of the solution and angles variation from point to point within the track, is compared with the ordinary BA, Differential Evolution (DE), PSO, K-ABC and Multi Objective PSO, and MO-PSO algorithm, which tackles the IK problem as a multi-objective optimization problem, where the algorithm aims to minimize the positional error and the angles variation simultaneously, as in [34].

The proposed path for this test, consists of 20 points distributed along a helical track within the workspace of the KUKA LBR iiwa 14 R820 as shown in Figure 5. The helical trajectory is defined by its center point from the bottom $C=\left[x_c=60 \mathrm{cm}, y_c=0 \mathrm{cm}, z_c=10 \mathrm{cm}\right]$ relative to the base of the robot, a radius of the helix $r=20 \mathrm{cm}$, its length $l=50$ cm, and a number of helices h = 2. In this test, the required orientation of the robot’s end-effector during the path tracking process is defined by Equation (27), which enables the end-effector to be normal to the base floor.

$${}^{0}R_7=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right]$$

(27)

The common parameters for all the compared algorithms used in this test are population size $N=50$, and the maximum number of iterations $t_{max}=200$. For reliability inception, the path tracking test is run 20 times and the average position error, angles variation and the computing taken time of the IK solution for each point of the trajectory is recorded as evaluation’s criteria of the performance of the IK-BA in comparison with the other algorithms, in path tracking task. Figure 6 shows the average error of the IK solutions for each point. Figure 7 displays the angles variation from point to point of the trajectory, noting that the angles variation of the 1st point is relative to the initial configuration of the robot. Figure 8 reveals the average calculation time taken to provide an IK solution for each point of the trajectory. Figure 9 represents the joints movement along the trajectory of the first run of the path tracking test.

Figure 6 demonstrates that the proposed IK-BA provides finer IK solutions in terms of accuracy and repeatability since it has the least average position error among the compared algorithms. Figure 7 shows that both the IK-BA and MO-PSO algorithms present IK solutions with the least angles variation from point to point of the trajectory. This is affirmed in Figure 9, which shows that IK-BA and MO-PSO generate smooth angles variation along the trajectory, in contrary to the uneven angles variation generated by the other used algorithms. In term of computing time, Figure 8 shows that PSO has the least computational time in 200 iterations, whereas the MO-PSO ranks last in term of computing time.

The performance of the compared algorithms is statistically analysed to give an accurate measurement of the effectiveness of IK-BA compared to other algorithms.

Table 3. Statistical analyses of the performance of the employed algorithms.

Algorithm	Performance	Mean	Median	Std	Maximum	Minimum
IK-BA	Position error [cm]	0.0019	0.0019	1.7824 × 10−4	0.0022	0.0016
	Angles variation [deg]	87.1916	87.5779	7.8247	101.8740	72.8021
	Computing time [s]	0.2439	0.2359	0.0229	0.3022	0.2195
BA	Position error [cm]	1.6849	1.5374	0.9069	3.9487	0.1830
	Angles variation [deg]	301.5794	301.5530	25.2466	345.7477	260.4859
	Computing time [s]	0.2800	0.2642	0.0335	0.3511	0.2462
DE	Position error [cm]	0.3653	0.2793	0.2929	0.9981	0.0387
	Angles variation [deg]	303.6915	297.8383	32.8049	377.9527	264.5489
	Computing time [s]	0.4154	0.3924	0.0617	0.5850	0.3534
PSO	Position error [cm]	0.4382	0.4019	0.3312	1.0384	0.0000
	Angles variation [deg]	274.7692	276.8618	21.9459	314.8747	223.3415
	Computing time [s]	0.2069	0.1946	0.0336	0.3029	0.1769
K-ABC	Position error [cm]	0.2318	0.2232	0.0722	0.3515	0.1270
	Angles variation [deg]	276.6456	273.5645	20.9967	311.8878	233.6387
	Computing time [s]	0.4591	0.4200	0.0926	0.7318	0.3980
MO-PSO	Position error [cm]	0.4075	0.3960	0.0891	0.7014	0.2871
	Angles variation [deg]	89.9602	86.8036	25.6630	164.1794	25.7217
	Computing time [s]	1.2853	1.2222	0.2282	1.9888	1.0135

shows the results of the statistical analyses of the employed algorithms for the path tracking test. A one-way analysis of variance (ANOVA) is also suggested to check if there is a statistical difference between the suggested algorithm and the other algorithms.

Therefore, to inspect the significant difference among the compared algorithms, the parametric statistical test ANOVA is employed, using 0.05 significance level.

Table 4. ANOVA test results for position error comparatives.

Source	SS	df	MS	F	Prob > F
Algorithms	35.0338	5	7.0068	40.7685	7.2225 × 10−24
Error	19.5929	114	0.1719
Total	54.6267	119

,

Table 5. ANOVA test for angles variation comparatives.

Source	SS	df	MS	F	Prob > F
Algorithms	1.0877 × 106	5	2.1754 × 105	388.9396	7.1885 × 10−70
Error	6.3761 × 104	114	559.3091
Total	1.1514 × 106	119

and

Table 6. ANOVA test for computing time comparatives.

Source	SS	df	MS	F	Prob > F
Algorithms	16.4700	5	3.2940	293.9766	2.1670 × 10−63
Error	1.2774	114	0.0112
Total	17.7474	119

, show the ANOVA test results for position error, angles variation and computational time, respectively.

The components of ANOVA tables are briefly explained as follows [51]:

• Source: stands for source of variation, which could be a variation between groups (algorithms’ sample results), or error variation, indicating a variation within the groups.
• SS: stands for sums of squares.
• df: degrees of freedom, which indicates the number of independent data.
• MS: mean squares, they are calculated by dividing sums of squares (SS) by their appropriate degrees of freedom.
• F: refers to the F-statistic, which used to test the null hypothesis that the means of several groups are equal. It is calculated as the ratio of the variation between the groups to the variation within the groups.
• Prop > F: indicates the p-value which is used to determine the level of significance of the F-statistic and helps to decide whether the null hypothesis of equal means is rejected or not.

The box plots for ANOVA tests are illustrated in Figure 10, Figure 11 and Figure 12. They represent the results of the performance of each algorithm in terms of position error, angles variation and computing time, respectively.

5.3. Time Assessment

Two-time assessments have been conducted to evaluate the potentials of IK-BA regarding time complexity and its real time potentials.

The first test was conducted using the following hardware configuration: Intel core i5 processor, with 8 Go RAM DDR3.

The time assessment for a possible evaluation of the IK-BA was done using the following configuration: intel i7 processor, with 12 Go RAM (4 Go DDR4 + 8 Go DIMM)

5.3.1. Time Assessment of the IK-BA for Real-Time Path Planning

The time assessment for a possible evaluation of the IK-BA was carried out using the following configuration: Intel i7 processor, with 4 Go RAM DDR4 + 8 Go DIMM memory. The system was running under MS Windows 11 and using MATLAB 2019. The test was performed using $N=12, 15 \mathrm{and}20$ bats with the same intrinsic parameters of IK-BA presented in Table 1.

All data and variables used were cleared from one test to another. Only the operating system and MATLAB were active during the tests. The termination criterion of the algorithm processing was to a positional error ${\epsilon }_p=0.1 \mathrm{cm}$. The algorithm had run 20 times, where a random target point was generated at each run. All target points were generated within the workspace of the KUKA LBR iiwa 14 R820 robot.

Statistical results of the time assessment test of the proposed IK-BA are presented in

Table 7. Statistical results of time assessment test (time is in ms).

Number of Bats	Mean Time in ms	Median	Std	Maximum	Minimum
12	48.01	45.8	10.1	67.2	32.6
15	61.09	64.10	12.85	92.35	41.81
20	83.75	73.03	28.96	80.06	47.88
30	180,94	101.29	149.67	448.44	69.52

. They show that IK-BA may be used as a real-time path planner when the number of agents is about 12 or less; meanwhile, for higher configuration results confirmed, it can be used offline only, meaning it can be used to plan the path before coding the robot.

5.3.2. Discussion on Time Processing of MH Inverse Solvers

The industrial processing consists in reproducing the path for a needed task, as part of a production plan. In general, the path planning is offline-generated, verified and tested through simulations, then applied, for a real validation on the target robot.

In Table 3, the computing time comparatives correspond to the processing time, where the termination criterion is defined by a maximum number of iteration $t_{max}=200$ iterations. The next section presents a time assessment of the IK-BA, where the termination criterion is set according to a maximum acceptable IK positional error.

The configuration used in test detailed by Table 3 is based on 20 agents for all algorithms. The termination criterion is used to ensure that the algorithm does not run indefinitely. The termination criterion can greatly affect the performance of the algorithm. Commonly, the termination criterion is defined by a maximum number of iterations, or according to an acceptable quality of the provided solution. The algorithm is stopped when the quality of the best solution found, while processing reaches a certain threshold [52].

Assuming this and in regards of the results in Table 7, it is evident that the proposed algorithm is suitable for offline path planning. It shows better time processing than its challengers, while not meeting a real-time path planning constraint (less than 50 ms). The potentials of the IK-BA as a real-time path planner are discussed in the next paragraph.

A real-time system is one that is capable to procure a deterministic response in a specified time frame. If this time frame is critical and must be strictly respected, the system is assumed to be a hard real-time system. In contrast, when the time may be subject to delay, the system is assumed to be a soft real-time system [53,54].

According to our previous experimentation, we proved that the proposal may respect a real-time constraint for each given configuration; the most effective number of bats is the one with 12 agents, where the response time is less than 50 ms.

6. Real Application

The proposed IK-BA is employed as an IK solver for the Dobot magician robotic arm. The Dobot magician robot presented in Figure 13 is a desktop 4-Dof robot that could be used for different tasks such as laser engraving, 3D printing, drawing and calligraphy.

In this application, the robot’s end-effector is intended to follow a circular trajectory discredited by 20 points distributed along the trajectory as shown in Figure 14. The circular path is defined by a center $C=\left[16 \mathrm{cm}, 16 \mathrm{cm}, 0\right].$ and a radius $r=7 \mathrm{cm}$, which is generated along the workspace of the Dobot magician robotic arm as revealed in Figure 15.

The points reached by the robot’s end-effector along with the actual trajectory points are presented in Figure 15. Figure 16 shows the corresponding generated joint angles for each point of the trajectory using IK-BA. ${\theta }_4$ is set to be equal to ${\theta }_1$, to obtain a fixed orientation of the end-effector. The positional error of the robot’s end-effector is relative to the trajectory points represented in Figure 17.

The proposed IK-BA provided accurate solutions of an error in a range of 10⁻¹⁴ cm. Figure 18 shows that the average error among the 20 points of the trajectory reaches 2.1196 × 10⁻¹⁴ cm along 1000 iterations.

Table 8. Generated joint angles and their corresponding error using IK-BA for the 20 points of the trajectory.

Trajectory Points	Joint Angles [deg]				Position Error [cm]
	${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{\theta }}_3$	${\mathit{\theta }}_4$
1	34.8	66.5	49.5	34.8	2.04 × 10−14
2	38.7	69.6	45.6	38.7	1.28 × 10−14
3	42.9	71.4	43.5	42.9	1.14 × 10−14
4	47.1	71.4	43.5	47.1	3.37 × 10−14
5	51.3	69.6	45.6	51.3	3.04 × 10−14
6	55.2	66.5	49.5	55.2	1.71 × 10−14
7	58.6	62.6	54.8	58.6	2.20 × 10−14
8	61.2	58.4	61.1	61.2	1.68 × 10−14
9	62.8	54.6	68.0	62.8	2.29 × 10−14
10	62.8	51.5	75.2	62.8	1.81 × 10−14
11	60.6	49.5	82.2	60.6	1.43 × 10−14
12	56.0	48.8	88.3	56.0	9.61 × 10−14
13	49.0	48.8	92.0	49.0	2.02 × 10−14
14	41.0	48.8	92.0	41.0	1.49 × 10−14
15	34.0	48.8	88.3	34.0	3.36 × 10−14
16	29.4	49.5	82.2	29.4	2.37 × 10−14
17	27.2	51.5	75.2	27.2	2.92 × 10−14
18	27.2	54.6	68.0	27.2	2.95 × 10−14
19	28.8	58.4	61.1	28.8	1.34 × 10−14
20	31.4	62.6	54.8	31.4	3.00 × 10−14

represents the twenty joint positions used for the test, the position error of the solver in respect to pen orientation. All returned solutions obtained a precision higher than the system needs and the measured observed positions can not exceed the mechanical frame constraints of the real robot as testified by its provider.

It is important to note that the proposed solver is retuning solutions which are much higher in precision than the system can mechanically perform, since the mechanical system precision is limited to 0.02 cm [55]. A solver should always produce solutions as precisely as possible, but usually where the precision is higher than what the mechanical system can perform.

7. Conclusions

A specially improved BA algorithm variant dedicated for the IK robotic problem is suggested in this paper. The proposed method is called IK-BA; it stands for a specific bat algorithm dedicated to robotic-arms’ inverse geometric solution, where the elongation control mechanism is embedded in its update equations. The proposed IK-BA is qualified of giving accurate solutions while guaranteeing minimal joints variation from the initial position towards the target point. The importance of this feature stands out in path tracking tasks, which enables generation of smooth angle variation from point-to-point of the trajectory. This improvement is established through adjusting the update equation in such a way that it takes into account the initial angular position of the robot. When the orientation of the end-effector is required in addition to a desired end-effector position, a manipulator decoupling method is suggested, where the three last joints of the manipulator are calculated deterministically; whereas, the remaining first joints are determined using the proposed IK-BA. This methodology is tested using the simulation of path tracking task of the KUKA LBR iiwa 14 R820, in which it is required to follow a helical trajectory. The effectiveness of the proposed approach is proved by a comparative study with the standard BA algorithm and other optimization algorithms from the literature, including DE, PSO, K-ABC and MO-PSO. The statistical analyses supported by a one-way ANOVA test show that the IK-BA has the best compromise between the accuracy of the solution, the angles variation and the computing time.