On the Use of a Genetic Algorithm for Determining Ho–Cook Coefficients in Continuous Path Planning of Industrial Robotic Manipulators

Abstract

Path planning is one of the key steps in the application of industrial robotic manipulators. The process of determining trajectories can be time-intensive and mathematically complex, which raises the complexity and error proneness of this task. For these reasons, the authors tested the application of a genetic algorithm (GA) on the problem of continuous path planning based on the Ho–Cook method. The generation of trajectories was optimized with regard to the distance between individual segments. A boundary condition was set regarding the minimal values that the trajectory parameters can be set in order to avoid stationary solutions. Any distances between segments introduced by this condition were addressed with Bezier spline interpolation applied between evolved segments. The developed algorithm was shown to generate trajectories and can easily be applied for the further path planning of various robotic manipulators, which indicates great promise for the use of such algorithms.

1. Introduction

The most significant element in the application of industrial robotic manipulators on realistic tasks is the path planning process [1]. This process determines the trajectory of the joint movements—their positions, speeds, and accelerations—allowing the robotic manipulator to perform operations within its environment [2]. The goal of path planning is to calculate paths that satisfy the preset conditions that have to be fulfilled in order to accomplish the task. There are two main paradigms of trajectory determination: point-to-point and continuous path planning [3,4]. Point-to-point planning concerns itself with generating paths between two points in space and is commonly used for operations such as the pick-and-place transfer of objects [5]. Continuous path planning, on the other hand, takes into account the movement not just between the initial and final points in space but also the positions and speeds between them. Such an approach is commonly used for tasks such as welding or painting objects in space [6]. Multiple deterministic methods can be used to perform the path planning of industrial robotic manipulators continuously: Ho–Cook [7], Taylor’s polynomial approach [8], or interpolation between trajectory points [9]. While these algorithms perform well, they can be computationally complex and, depending on the mode of application, error-prone. Evolutionary computing is a branch of artificial intelligence that deals with the study of algorithms that imitate natural processes [10]. The basic algorithm in this area is the so-called genetic algorithm, which, by its design, imitates the natural process of evolution [11]. Evolutionary computing algorithms have been shown to have many uses in robotics [12].

Shukla et al. (2021) [13] demonstrated the use of evolutionary computing for robotic grasp manipulation. The authors applied EC algorithms to assist in training deep learning models, in an approach known as hybrid models. They managed to achieve high-precision models with this approach. Ferigo et al. (2021) [14] applied evolutionary computing algorithms for evolving sensory apparatus in soft computing applications. Kim et al. (2021) [15] demonstrated the application of evolutionary computing for the issue of quadruped robot gait optimization. The authors utilized GA to create paths for mobile robots’ legs to achieve a more controlled gait. Liu et al. (2022) [16] created a digital twin, which is a virtualized copy of a real robot. The authors then demonstrated the ability to apply GA for the path planning of such a robot, which they managed to successfully transfer to a real robot. Li et al. (2021) [17] addressed another important issue in robotics: task allocation through the application of a differential evolution algorithm. The authors managed to achieve state-of-the-art results on a multitask optimization problem. Task allocation was also addressed by Martin et al. (2021) [18]. In this paper, the authors specifically used GA to distribute the tasks amongst robots based on the nonlinear branching criteria. Path planning for robots was also addressed by Hao et al. (2021) [19]. The authors utilized GA to optimize a path concerning the possible collision risks, achieving paths that are capable of avoiding obstructions in the space. A similar approach was demonstrated by Rahmaniar and Rakhmania (2022) [20] for mobile robots. The authors demonstrated that GA-optimized paths are capable of achieving smoother paths when compared to classical methods. Tuning the paths based on Bezier splines can also be achieved using evolutionary computing algorithms, as demonstrated by Song et al. (2021) [21], who applied particle swarm optimization to determine the parameters of the splines. Li et al. (2022) [22] demonstrated the subgoal hybrid planning for the path smoothing of paths obtained with forward search optimization in contrast to Djikstra, A*, D*, and D*-lite algorithms. The authors compared the performance of forward search optimization and the newly developed subgoal-based hybrid path planning algorithm, demonstrating that smooth global paths can be achieved with such an approach. A similar approach to the one presented in this paper, where a deterministic algorithm A* is combined with an evolutionary approach, namely the coevolutionary algorithm, was adapted by Garcia et al. (2023) [23]. The authors demonstrated that such an approach has a high enough performance for applications on edge nodes, and performs well in conditions where alternatives such as M* or WHCA would fail to generate valid paths. Yu et al. (2023) [24] applied the artificial bee colony to optimize a multi-objective path planning issue. The authors demonstrated that such an algorithm can be successfully applied to multi-objective problems, namely path efficiency and path security. Another nature-inspired algorithm was shown by Wu et al. [25], who applied the ant colony algorithm to the path planning of the mobile robot. The modified version of the ant colony algorithm that the authors developed shows a significant improvement in comparison to the state-of-the-art methods. Lou et al. (2023) [26] applied a graphical computing method for the problem of continuous path planning for welding. The simulations performed by the authors on the generated paths demonstrate the possibility of applying the investigated method in real-word applications. Another industrial production application, namely surface grinding, was discussed by Li et al. (2023) [27]. The authors applied a revised Levenberg–Marquardt and differential evolution hybrid algorithm, with real-world validation on the problem of grinding titanium blades.

While the state-of-the-art research shows many applications of evolutionary computing and GA in the area of robotics and path planning, none address the combination of algorithms such as Ho–Cook with GA or similar evolutionary algorithms.

The Ho–Cook algorithm was selected in particular as the topic of the research due to several advantages that it possesses compared to similar algorithms: mainly the ability to manually select as many points as desired (allowing for the granular control of the trajectory precision) and the fact that it includes the orientation of the tool, by design—as will be shown in Section 2. Due to the lack of previous research focusing on the Ho–Cook algorithm parameter tuning with evolutionary algorithms, the authors selected GA as the second focus of the research. GA is the basic evolutionary algorithm, which means that it should serve as a good indicator of performance for more advanced algorithms [28,29]. The gene setup that was developed and is presented in Section 2.2 of this manuscript is also novel, and customized to the Ho–Cook problem. This chromosome encoding describes the Ho–Cook parametrization and may be used as the basis for the further research of additional, more advanced evolutionary algorithms, as most algorithms of this kind will require this encoding to be performed in the same manner [12,30].

The goal of this paper was to test whether the process of path planning in a continuous environment, based on the Ho–Cook algorithm, can be simplified through the application of the GA. In addition, the parameters of GA that provide the best performance were also determined. In this approach, the common issues in continuous path planning are addressed through the Ho–Cook algorithm, whose shortcoming is the analytical complexity of the coefficient determination, as will be shown in Section 2.

2. Materials and Methods

This section will serve to present the basic idea of the Ho–Cook path planning process to point out which part of it will be tuned using the GA-based approach. Then, the process of GA development will be described.

2.1. Ho–Cook Path Planning

The Ho–Cook path planning algorithm is a continuous path planning algorithm. It is based on determining n points that will be the elements of the trajectory. In the case of the obstacles being present in the tool space inside of which the path planning is being performed, the aforementioned points should be placed in such a way that obstacles are not present, as obstacle avoidance is not a built-in feature of the Ho–Cook method [31]. Due to the Ho–Cook method having to pass through the points defined initially, as they are starting/ending points of the segments, the movement from the source to the destination is guaranteed. The Ho–Cook method works by interpolating between the aforementioned points of the desired trajectory to achieve the shortest possible path through the calculation of in-between segments. Higher-order polynomials (4th and 3rd) are used to assure the smoothness of the final trajectory. Then, the $n-1$ segments between the points are interpolated using the polynomials given as [32]:

$$q_k\left(t\right)=B_{0k}+B_{1k}t+B_{2k}{t}^2+B_{3k}{t}^3+B_{4k}{t}^4$$

(1)

for the first and the last segments, whereas the other segments are interpolated using:

$$q_k\left(t\right)=B_{0k}+B_{1k}t+B_{2k}{t}^2++B_{3k}{t}^3$$

(2)

The speeds and accelerations can be derived from the above. For the first and the last segment, they are given as:

$$\dot{q_k\left(t\right)}=B_{1k}+2B_{2k}t+3B_{3k}{t}^2+4B_{4k}{t}^3,$$

(3)

and

$$\ddot{q_{k}t}=2B_{2k}+6B_{1k}t+12B_{4k}{t}^2$$

(4)

respectively. The speed and acceleration for the other segments are given similarly with:

$$\dot{q_k\left(t\right)}=B_{1k}+2B_{2k}t+3B_{3k}{t}^2,$$

(5)

and

$$\ddot{q_{k}t}=2B_{2k}+6B_{1k}t+12B_{4k}{t}^2.$$

(6)

In the above equations, the symbols used are as follows:

• k—the trajectory point;
• m—the total number of trajectory points, and
• B—coefficients of the interpolation polynomials.

For each of the segments, a number of coefficients need to be determined. They can be defined within the matrix of the shape $k\times n$, where n is the polynomial degree (third or fourth) and k is the number of the joints of the robotic manipulator. In the presented research, the number of joints was assumed to be six, since this is a common number of degrees of freedom for industry-standard articulated robots [33].

To determine the above, the Dennavit–Hartenberg (D-H) algorithm was performed. First, the simplified kinematic schematic with the noted rotation axes was designed for the robot that is being modeled. An example of the ABB IRB 120 robot can be seen in Figure 1. On the schematic, the joints of the robot have orthonormal coordinate systems adjoined, with the orientation dependent on the rotation axes of the robot in question [34].

This will allow us to determine the kinematic properties of the robot, which can be read from the schematic: joint distance d, joint angle $\theta$, link length a, and link rotation angle $\alpha$ [35]. The obtained kinematics properties for the robot in question are given in

Table 1. The kinematic parameters of the analyzed robotic manipulator.

$\mathbf{\Theta }$ [rad]	d [mm]	a [mm]	$\mathit{\alpha }$ [rad]
${\mathrm{\Theta }}_1=q_1$	$d_1=290$	$a_1=0$	${\alpha }_1=-\pi /2$
${\mathrm{\Theta }}_2=q_2$	$d_2=0$	$a_2=270$	${\alpha }_2=0$
${\mathrm{\Theta }}_3=q_3$	$d_3=0$	$a_3=70$	${\alpha }_3=-\pi /2$
${\mathrm{\Theta }}_4=q_4$	$d_4=302$	$a_4=0$	${\alpha }_4=\pi /2$
${\mathrm{\Theta }}_5=q_5$	$d_5=0$	$a_5=0$	${\alpha }_5=-\pi /2$
${\mathrm{\Theta }}_6=q_6$	$d_6=72$	$a_6=0$	${\alpha }_6=0$

.

The next step is to obtain the kinematic transformation matrix, $T_0^6$, which is the product of each individual joint transformation matrix [35]:

$$T_0^6={\mathrm{\Pi }}_{k=1}^6\left[\begin{array}{cccc}cos{\theta }_k& -cos{\alpha }_{k}sin{\theta }_k& sin{\alpha }_{k}sin{\theta }_k& a_{k}cos{\theta }_k\\ sin{\theta }_k& cos{\alpha }_{k}cos{\theta }_k& -sin{\alpha }_{k}cos{\theta }_k& a_{k}sin{\theta }_k\\ 0& sin{\alpha }_k& cos{\alpha }_k& d_k\\ 0& 0& 0& 1\end{array}\right]$$

(7)

Inserting the values from Table 1 into the above equation yields the following equation that describes the transformation matrix:

$$\begin{array}{c}\hfill T_0^6=\left[\right[1.0\ast ((sin\left(q_1\right)\ast sin\left(q_4\right)+cos\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3))\ast cos\left(q_5\right)-sin\left(q_5\right)\\ \hfill \ast sin(q_2+q_3)\ast cos\left(q_1\right))\ast cos\left(q_6\right)+1.0\ast (sin\left(q_1\right)\ast cos\left(q_4\right)-1.0\ast sin\left(q_4\right)\ast cos\left(q_1\right)\\ \hfill cos(q_2+q_3))\ast sin\left(q_6\right),1.0\ast (-(sin\left(q_1\right)\ast sin\left(q_4\right)+cos\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3))\\ \hfill cos\left(q_5\right)+sin\left(q_5\right)\ast sin(q_2+q_3)\ast cos\left(q_1\right))\ast sin\left(q_6\right)+1.0\ast (sin\left(q_1\right)\ast cos\left(q_4\right)\\ \hfill -1.0\ast sin\left(q_4\right)\ast cos\left(q_1\right)\ast cos(q_2+q_3))\ast cos\left(q_6\right),-1.0\ast (sin\left(q_1\right)\ast sin\left(q_4\right)\\ \hfill +cos\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3))\ast sin\left(q_5\right)-1.0\ast sin(q_2+q_3)\ast cos\left(q_1\right)\ast cos\left(q_5\right),\\ \hfill -0.072\ast sin\left(q_1\right)\ast sin\left(q_4\right)\ast sin\left(q_5\right)-0.072\ast sin\left(q_5\right)\ast cos\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3)\\ \hfill -0.072\ast sin(q_2+q_3)\ast cos\left(q_1\right)\ast cos\left(q_5\right)-0.302\ast sin(q_2+q_3)\ast cos\left(q_1\right)\\ \hfill +0.29\ast cos\left(q_1\right)\ast cos\left(q_2\right)+0.07\ast cos\left(q_1\right)\ast cos(q_2+q_3)],[1.0\ast \left(\right(sin\left(q_1\right)\ast cos\left(q_4\right)\\ \hfill cos(q_2+q_3)-sin\left(q_4\right)\ast cos\left(q_1\right))\ast cos\left(q_5\right)-sin\left(q_1\right)\ast sin\left(q_5\right)\ast sin(q_2+q_3))\\ \hfill cos\left(q_6\right)-1.0\ast (sin\left(q_1\right)\ast sin\left(q_4\right)\ast cos(q_2+q_3)+cos\left(q_1\right)\ast cos\left(q_4\right))\ast sin\left(q_6\right),\\ \hfill 1.0\ast ((-sin\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3)+sin\left(q_4\right)\ast cos\left(q_1\right))\ast cos\left(q_5\right)+sin\left(q_1\right)\\ \hfill sin\left(q_5\right)\ast sin(q_2+q_3))\ast sin\left(q_6\right)-1.0\ast (sin\left(q_1\right)\ast sin\left(q_4\right)\ast cos(q_2+q_3)+cos\left(q_1\right)\ast cos\left(q_4\right))\\ \hfill cos\left(q_6\right),1.0\ast (-sin\left(q_1\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3)+sin\left(q_4\right)\ast cos\left(q_1\right))\ast sin\left(q_5\right)\\ \hfill -1.0\ast sin\left(q_1\right)\ast sin(q_2+q_3)\ast cos\left(q_5\right),-0.072\ast sin\left(q_1\right)\ast sin\left(q_5\right)\ast cos\left(q_4\right)\ast cos(q_2+q_3)\\ \hfill -0.072\ast sin\left(q_1\right)\ast sin(q_2+q_3)\ast cos\left(q_5\right)-0.302\ast sin\left(q_1\right)\ast sin(q_2+q_3)\\ \hfill +0.29\ast sin\left(q_1\right)\ast cos\left(q_2\right)+0.07\ast sin\left(q_1\right)\ast cos(q_2+q_3)\\ \hfill +0.072\ast sin\left(q_4\right)\ast sin\left(q_5\right)\ast cos\left(q_1\right)],[-1.0\ast (sin\left(q_5\right)\ast cos(q_2+q_3)+sin(q_2+q_3)\ast cos\left(q_4\right)\\ \hfill cos\left(q_5\right))\ast cos\left(q_6\right)+1.0\ast sin\left(q_4\right)\ast sin\left(q_6\right)\ast sin(q_2+q_3),1.0\ast (sin\left(q_5\right)\ast cos(q_2+q_3)\\ \hfill +sin(q_2+q_3)\ast cos\left(q_4\right)\ast cos\left(q_5\right))\ast sin\left(q_6\right)+1.0\ast sin\left(q_4\right)\ast sin(q_2+q_3)\ast cos\left(q_6\right),\\ \hfill 1.0\ast sin\left(q_5\right)\ast sin(q_2+q_3)\ast cos\left(q_4\right)-1.0\ast cos\left(q_5\right)\ast cos(q_2+q_3),-0.29\ast sin\left(q_2\right)\\ \hfill +0.072\ast sin\left(q_5\right)\ast sin(q_2+q_3)\ast cos\left(q_4\right)-0.07\ast sin(q_2+q_3)\\ \hfill -0.072\ast cos\left(q_5\right)\ast cos(q_2+q_3)-0.302\ast cos(q_2+q_3)+0.29],[0,0,0,1]]\end{array}$$

(8)

The above matrix defines the equations that may be used to calculate the position $\left[xyz\right]$ and orientation $\left[\varphi \theta \psi \right]$. These equations are given as:

$$\begin{array}{c}\hfill x=-0.072·sin\left(q_1\right)·sin\left(q_4\right)·sin\left(q_5\right)-0.072·sin\left(q_5\right)·cos\left(q_1\right)·cos\left(q_4\right)·cos(q_2+q_3)\\ \hfill -0.072·sin(q_2+q_3)·cos\left(q_1\right)·cos\left(q_5\right)-0.302·sin(q_2+q_3)·cos\left(q_1\right)\\ \hfill +0.29·cos\left(q_1\right)·cos\left(q_2\right)+0.07·cos\left(q_1\right)·cos(q_2+q_3),\end{array}$$

(9)

$$\begin{array}{c}\hfill y=-0.072·sin\left(q_1\right)·sin\left(q_5\right)·cos\left(q_4\right)·cos(q_2+q_3)-0.072·sin\left(q_1\right)·sin(q_2+q_3)·cos\left(q_5\right)\\ \hfill -0.302·sin\left(q_1\right)·sin(q_2+q_3)+0.29·sin\left(q_1\right)·cos\left(q_2\right)+0.07·sin\left(q_1\right)·cos(q_2+q_3)\\ \hfill +0.072·sin\left(q_4\right)·sin\left(q_5\right)·cos\left(q_1\right),\end{array}$$

(10)

$$\begin{array}{c}\hfill z=-0.29·sin\left(q_2\right)+0.072·sin\left(q_5\right)·sin(q_2+q_3)·cos\left(q_4\right)-0.07·sin(q_2+q_3)\\ \hfill -0.072·cos\left(q_5\right)·cos(q_2+q_3)-0.302·cos(q_2+q_3)+0.29,\end{array}$$

(11)

for the linear coordinates in the tool space, whereas the orientation is defined by:

$$\begin{array}{c}\hfill \varphi =-1.0·(sin\left(q_1\right)·sin\left(q_4\right)+cos\left(q_1\right)·cos\left(q_4\right)·cos(q_2+q_3))·sin\left(q_5\right)\\ \hfill -1.0·sin(q_2+q_3)·cos\left(q_1\right)·cos\left(q_5\right),\end{array}$$

(12)

$$\begin{array}{c}\hfill \theta =1.0·(-sin\left(q_1\right)·cos\left(q_4\right)·cos(q_2+q_3)+sin\left(q_4\right)·cos\left(q_1\right))·sin\left(q_5\right)\\ \hfill -1.0·sin\left(q_1\right)·sin(q_2+q_3)·cos\left(q_5\right),\end{array}$$

(13)

and

$$\begin{array}{c}\hfill \psi =-0.072·sin\left(q_1\right)·sin\left(q_5\right)·cos\left(q_4\right)·cos(q_2+q_3)\\ \hfill -0.072·sin\left(q_1\right)·sin(q_2+q_3)·cos\left(q_5\right)\\ \hfill -0.302·sin\left(q_1\right)·sin(q_2+q_3)+0.29·sin\left(q_1\right)·cos\left(q_2\right)\\ \hfill +0.07·sin\left(q_1\right)·cos(q_2+q_3)+0.072·sin\left(q_4\right)·sin\left(q_5\right)·cos\left(q_1\right)\end{array}$$

(14)

The above equations can be transformed in order to obtain the joint values for the given position and orientation. This is how the Ho–Cook trajectory planning process is capable of taking into account the orientation and the position, as the points that are defined as the points of the trajectory involve orientation as well as the position $\left[xyz\varphi \theta \psi \right]$ [31]. The process of defining the Ho–Cook trajectory can then be followed by defining the following relation:

$$\begin{array}{c}\hfill \left[\begin{array}{c}\dot{q_2}\\ \dot{q_3}\\ ⋮\\ \dot{q_{m-2}}\\ \dot{q_{m-1}}\end{array}\right]\left[\begin{array}{cccccc}\frac{3}{t_2}+\frac{2}{t_3}& t_4& 0& \cdots & 0& 0\\ \frac{1}{t_3}& 2·(t_3+t_4)& t_5& \cdots & 0& 0\\ 0& t_3& 2·(t_4+t_5)& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & t_{m-1}& 0\\ 0& 0& 0& \cdots & 2·(t_{m-2}+t_{m-1})& \frac{t}{t_{m-1}}\\ 0& 0& 0& \cdots & t_{m-2}& \frac{2}{t_{m-1}}+\frac{3}{t_m}\end{array}\right]\\ \hfill =\left[\begin{array}{c}\frac{6}{t^2}(q_2-q_1)+\frac{3}{t_3^2}(q_3-q_2)\\ \frac{3}{t_3{t}_4}[t_3^2(q_4-q_3)+t_4^2(q_3-q_2)]\\ \frac{3}{t_4{t}_5}[t_4^2(q_5-q_4)+t_5^2(q_4-q_3)]\\ \frac{6}{t_{m-2}{t}_{m-1}}[t_{m-2}^2(q_{m-1}-q_{m-2})+t_{m-1}^2(q_{m-2}-q_{m-3})]\\ \frac{6}{t_m^2}(q_m-q_{m-1})+\frac{3}{t_{m-1}^2}(q_{m-1}-q_{m-2})\end{array}\right]\end{array}$$

(15)

Finally, the segment coefficients of the Ho–Cook trajectories are then calculated using three different equation sets. The first segment coefficients are calculated as:

$$\left[B_1^0{B}_1^1{B}_1^2{B}_1^3{B}_1^4\right]=\left[q_1{q}_2\dot{q_1}\dot{q_2}\right]·\left[\begin{array}{ccccc}1& 0& 0& -\frac{4}{t_2^3}& \frac{3}{t_2^4}\\ 0& 0& 0& \frac{4}{t_2^3}& -\frac{3}{t_2^4}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& -\frac{1}{t_2^2}& \frac{1}{t_2^3}\end{array}\right].$$

(16)

The segments between the first and the last one can be calculated using:

$$\left[B_k^0{B}_k^1{B}_k^2{B}_k^3\right]=\left[q_k{q}_{k+1}\dot{q_k}\dot{q_{k+1}}\right]·\left[\begin{array}{cccc}1& 0& -\frac{3}{t_{k+1}^2}& \frac{2}{t_{k+1}^3}\\ 0& 0& \frac{3}{t_{k+1}^2}& -\frac{2}{t_{k+1}^3}\\ 0& 1& -\frac{2}{t_{k+1}}& \frac{1}{t_{k+1}^2}\\ 0& 0& -\frac{1}{t_{k+1}}& \frac{1}{t_{k+1}^2}\end{array}\right].$$

(17)

The final segment of the interpolated trajectory can then be defined as:

$$\left[B_{m-1}^0{B}_{m-1}^1{B}_{m-1}^2{B}_{m-1}^3{B}_{m-1}^4\right]=\left[q_{m-1}{q}_m\dot{q_{m-1}}\dot{q_m}\right]·\left[\begin{array}{ccccc}1& 0& -\frac{6}{t_m^2}& \frac{8}{t_m^3}& -\frac{3}{t_m^4}\\ 0& 0& \frac{6}{t_m^2}& -\frac{8}{t_m^3}& \frac{3}{t_m^4}\\ 0& 0& -\frac{3}{t^2}& \frac{3}{t_2^2}& \frac{1}{t_2^3}\\ 0& 0& 0& 0& 0\end{array}\right]$$

(18)

These equations also allow for the calculation of the joint positions and speeds if the coefficients are known, which will be the focus of the paper going forward. The motion of the robot is only limited by the natural limitations of the robotic manipulator at hand (for ABB IRB 120, first joint in range <−2.88, 2.88> [rad], second joint <−1.92, 1.92> [rad] third joint <−1.22, 1.92> [rad], fourth joint <−2.79, 2.79> [rad], fifth joint <−2.09, 2.09> [rad], and <0, 6.28> [rad] for the final joint) [36], or other robot positional issues, such as singularities [37]. These issues can be addressed within the GA by artificially lowering the fitness of the solutions that are outside of the bonds.

2.2. Genetic Algorithm

A genetic algorithm (GA) is an optimization algorithm based on the natural evolution process [38]. The algorithm works by generating a population of randomly selected potential solutions [28]. Then, each of these solutions is evaluated individually according to a pre-defined function: the so-called fitness function [39]. The iterating process of the GA then starts by randomly selecting the potential solutions and performing the evolutionary operations on them, which generates a new solution set [40]. This process will repeat until a satisfactory solution is achieved [16]. The process is based on the evolutionary operations, of which, there are three in GA: crossover, mutation, and reproduction. Crossover was performed on two of the randomly selected potential solutions, and they were combined into a new solution [41]. This process allows for the newly generated solutions, and when the initial random selection process is performed by weighting it towards the better-performing solutions, previous research shows that newly generated solutions tend toward the optimal solution [42]. Still, there are two issues that the crossover operation can introduce: the convergence into local optima, and the loss of quality solutions [28]. The first one was addressed by introducing the mutation operation, which will randomly modify a single randomly selected solution. This allows us to check a wider area within the possible solutions [43]. The loss of quality solutions refers to the phenomena in which a good solution is, through the application of crossover and mutation, replaced with a worse one. To address this, the reproduction operation simply transferred a quality solution into the next solution set [44]. From the above, it can be concluded that multiple values need to be defined or tested to define how the GA will be applied:

• Shape of the potential solutions;
• The way in which the crossover and mutation will be applied;
• The probabilities with which the evolutionary operations will occur;
• The fitness function that will evaluate the solutions;
• The number of iterations (generations) of the algorithm;
• The number of candidate solutions in the algorithm;
• The manner of the solution selection for the operations.

Each of these elements, in the context of the paper’s goal, will be further discussed.

2.2.1. Solution Construction

The shape of the individual solutions first needs to be determined, as it is the basis for defining the remaining elements of GA. In the discussed problem, it was decided that a six-point trajectory would be used; in other words, five segments need to be generated in the shape of the vectors $B_1$, $B_2$, $B_3$, $B_4$, $B_5$. According to the previous subsection, $B_1$ and $B_5$ consist of five elements and the remaining ones consist of four elements. For simplicity, all of these elements were joined in a single vector that can be used for the evolutionary computation operations as:

$$B=\left[B_1^0{B}_1^1{B}_1^2{B}_1^3{B}_1^4{B}_2^0{B}_2^1{B}_2^2{B}_2^3{B}_3^0{B}_3^1{B}_3^2{B}_3^3{B}_4^0{B}_4^1{B}_4^2{B}_4^3{B}_5^0{B}_5^1{B}_5^2{B}_5^3{B}_5^4\right].$$

(19)

Each of the initial population candidate solutions were uniformly randomly filled with values in the range of $[-5,5]$. Trajectory points were set randomly.

2.2.2. Application of Evolutionary Computing Operations

The three aforementioned operations—crossover, mutation, and reproduction—need to be defined, as their manner of implementation is an important element. Each of the operations also has a probability of occurring, with the previous research in the area indicating that the probability of the crossover occurring should be very high (>80%), whereas the mutation should be low (<5%), with the remaining iterations being fulfilled with the reproduction [45]. Reproduction is the simplest, as it only copies the potential solution between the generations. Its occurrence was set to 7%. The manner in which the crossover was implemented was the so-called random crossover. As shown in Equation (19), each of the segments has an individual set. The algorithm of crossover was then performed on two candidate solutions. For each of the five segments, one of the two candidate solutions was selected, and its segment was inserted into the new solution. This process is shown in Figure 2. This yields a solution that is the combination of segments from previous solutions.

The crossover operation was performed with the $90\%$ probability. Finally, the mutation needs to be performed. The mutation will randomly replace an entire segment of the randomly selected solution, as shown in Figure 3. Without this mechanic, due to how the crossover operation is performed, only the initially generated segments are present in all of the generations of the algorithm. The occurrence rate of this operation was set to $3\%$.

2.2.3. The Fitness Function

To determine the quality of the solution, a fitness function needs to be determined. Such a function should be simple to calculate to allow for a fast calculation and execution of the algorithm while providing a realistic metric of the agent performance [46]. As can be seen in the previous Figure 2 and Figure 3, the segments generated by the GA are not continuous with each other. This is one of the main points that need to be addressed when trajectory planning is performed, and, as such, this was selected to represent the measure of quality for the candidate solution. An illustration of these distances is given in Figure 4.

To calculate this, the sum of all of the vertical distances between the segments was considered and calculated. In other words, if the last value of segment $B_k$ is given as $B_k\left[i\right]$ and the first value of the following segment $B_{k+1}$ is given as $B_{k+1}\left[0\right]$, the fitness function $\mathcal{F}$ can be defined as:

$$\mathcal{F}=\sum _{k=0}^4|B_k\left[i\right]-B_{k+1}\left[0\right]|.$$

(20)

Due to this being a minimization problem, the lower value of the fitness function indicates a higher-quality solution. It is important to note that boundary conditions need to be set. The initial testing showed that, when allowing the GA to minimize the distance completely, solutions will tend to a stationary trajectory, without any movement amongst the segments. In other words, the elements of the matrix given in Equation (19) will converge to zero. As this is not the desired outcome, the boundary condition was introduced. The boundary condition stated defines that all candidate solutions must have a total sum of the trajectory parameters equal to or higher than a certain value, which will be determined through testing. All of the candidate solutions that do not satisfy this condition will be removed from the population, and replaced with another randomly generated solution that does satisfy the aforementioned condition.

2.2.4. Candidate Solution Selection

To achieve an improvement across the generations of the GA, the candidate solutions selected for the application of evolutionary operations need to be selected wisely. Only if the solutions selected are of high quality can the improvement in the overall population fitness be expected. To achieve this, a fitness proportionate selection was used. This refers to the type of randomized selection in which fitness is of a higher quality. The type of fitness proportional selection used in the presented work was the so-called roulette wheel selection. In this type of selection, the probability of selecting a certain candidate solution is equal to its ratio to the overall fitness. If the total fitness is ${\mathcal{F}}_{total}$ and the individual fitness of candidate solution A is ${\mathcal{F}}_A$, then the probability of selecting it would be equal to:

$$p_A=\frac{{\mathcal{F}}_A}{{\mathcal{F}}_{total}}.$$

(21)

Since, in the presented case, the smaller fitness indicates the better solution, the above method can be applied to calculating a vector that determines the probability of each candidate solution, sorted in descending order. The candidate solution’s probabilities are then sorted in reverse order and the individual probability is assigned to each. This process is illustrated in Figure 5.

With all of these elements defined, it can be stated that the GA used in this research is the GA based on crossover and mutation with the roulette-wheel type selection and fixed fitness function [47].

2.3. Interpolation

Previously, when discussing the fitness value used during the presented research, it was mentioned that the lower bound needed to be set on the fitness that the solutions can achieve in order to prevent the stationary solutions. This approach yields solutions that will, even when optimized fully to the extent of the GA possibility, have vertical gaps between segments, meaning that a trajectory generated in such a manner cannot be considered as continuous. To address this, an interpolation technique was introduced. As the length of each segment is defined as 1 [s], the starting and ending 0.1 [s] of each segment were deleted. The first and final segments only had their final or initial 0.1 [s] element removed, respectively. This yields four missing segments, with a length of 0.2 s that need to be interpolated.

Interpolation was performed using second-degree Bezier splines [48], defined as [49]:

$$B_{i,k}\left(x\right)=\frac{x-t_i}{t_{i+k}-t_i}{B}_{i,k-1}\left(x\right)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}}{B}_{i+1,k-1}\left(x\right),$$

(22)

and

$$B_{i,0}\left(x\right)=\left\{\begin{array}{cc}1,\hfill & x\in [t_i,t_{i+1}>\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.$$

(23)

In the equations, t is the number of Bezier spline nodes, c is the coefficient of the spline determined based on the values of segments being connected, and k is the spline degree, equal to 2. An example of the result concerning the spline is given in the following section.

3. Results and Discussion

In this section, the results of the applied methodology will be discussed. The results in the process determination of boundary conditions for the obtained matrices will be demonstrated and commented on, followed by the results of GA parameter testing.

3.1. Determining the Optimal Boundary Condition

As mentioned previously, the boundary conditions were set as the minimal value of the sum of all elements of the candidate solution vector as given in Equation (19). Five different values of the bound were tested: no boundary conditions, 0.2, 0.1, 0.05, and 0.005. For the testing of the bound, the population size of the candidate solutions and the number of generations were both set to 10. While these parameter values are too small to achieve any significant results, as shown by further testing, they perform well enough to indicate the performance of the boundary condition tested without requiring too much execution time. For each of the tests performed, the results will be given as a visualization of the best-performing candidate solution at the end of the optimization process (joint positions, speeds, and accelerations given), along with a graph showing the change in the fitness of the best solution through the optimization process.

The first test was performed without any boundary condition set, and the results are shown in Figure 6. As can be seen, there are minimal movements of the robotic manipulator present, with the joint movement ranging from 0.1 to −0.25.

The first introduced boundary condition has a value of 0.2 and the results are shown in Figure 7. Here, it can be seen that the range of movement is higher, but, observing the graph of fitness through generations, we can see that the value to which the algorithm converges is very high.

Similar but slightly improved results are achieved with the bound set to 0.1, seen in Figure 8. The range of the movement is kept, but the algorithm converges to a better overall solution. The improved results achieved by setting this value indicate the need to continue testing lower values.

When the boundary condition is set to the value of 0.05, the graphs given in Figure 9 indicate that the movement range is kept with a further improvement in the convergence value of optimization, with the distances between segments below 0.05 [rad] and the range of motion above 1.0 [rad].

The final value for the boundary condition tested was 0.005, as demonstrated in Figure 10. This value shows a similar behavior to the configuration in which no boundary condition is set. This indicates that values this low should not be considered.

The tests performed indicate that lowering the boundary condition improves the results of the GA, as long as the condition is not set too low. When the condition is set too low, the GA performs in the same manner as it does when no condition is set. For this reason, the boundary condition of 0.05 was selected for further testing.

3.2. Determining the GA Parameters

After the boundary condition was determined and set to 0.05, the main parameters of the GA—the population size and the number of generations that the algorithm will run for—needed to be determined. Three separate configurations were tested:

• Population size of 100 executed for 100 generations;
• Population size of 1000 executed for 50 generations;
• Population size of 10,000 executed for 20 generations.

The first configuration, with 100 candidate solutions and 100 generations, is the least memory-intensive and computationally complex due to the fact that it has the lowest population value. It is shown to achieve a fitness function of 1.001, which is not satisfactory in the context of the problem. The achieved solution shows the large distances between multiple sections of the trajectory. The results are presented in Figure 11.

An increase in the population size to 1000 yields a significant increase in performance, even with the lower generation bound of 50. As shown in Figure 12, the lower number of generations does not negatively affect the performance, as the algorithm is shown to converge significantly before the fiftieth generation. The lowest achieved fitness function value is 0.33, which may be considered satisfactory.

The final tested configuration has 10,000 candidate solutions and was optimized for 20 generations, the results of which are shown in Figure 13. The trend from the previous configuration continues, as the larger configuration achieves a significantly improved result of 0.098. The GA converges between generations 12 and 15, indicating that 20 generations selected for this algorithm are enough for it to converge, despite the larger population size.

To discuss the overall results, it is shown that an increase in population size yields significant performance increases. The number of generations seems to have less of an influence, and can be more limited in further research, as the algorithm tends to converge to a solution around the 15th generation in all of the tested cases.

Execution Time

An important consideration for the application of algorithms is the execution times of the algorithms in tested configurations. The tests were performed on all the previously tested configurations and averaged across 10 runs. The configuration used for testing was a laptop computer with CPU Intel(R) Core(TM) i7-1065G7 CPU, with the CPU clock locked at 1.30 GHz for the test. The machine was equipped with 16 GB of RAM. The code was executed in a single-threaded mode. The results of the test are given in

Table 2. The execution times of various configurations, averaged over ten runs, with an average time per generation and the total execution time given.

Population	Generations	Total Time [s]	Average Time per Generation [s]
100	100	23.4	0.234
1000	50	119.5	2.39
10,000	20	403.4	20.17

.

It can be seen that the average time per each agent in the population remains relatively uniform (around 2$·{10}^{-3}$ s). This means that the increase in population size has a significant influence on execution times. The smallest configuration, which yields the poorest results, finishes the execution in around twenty seconds. The two following configurations take almost two minutes, or slightly below seven minutes, respectively.

3.3. Result Illustration

The path as given in the Figure 14 represents the motion of the joint over the course of five seconds between the position of −0.47 [rad] to 1.0 [rad]. The generated path is smooth, without any sudden changes in the joint motion direction, owing to the interpolation accomplished with Bezier splines. Some minor changes to the torque and negative vibrations may be present due to the changes in the direction present between the points, but as these transitions are smooth, these effects should not be overly negative on the motion of the path. The motion not being monotonous may have a negative effect in the sense of using more energy than would be necessary if manually tuned coefficients were used, but this negative effect should be minor and outweighed by the complexity of manually tuning the parameters, except in situations where lowering the energy use is crucial. The motion is continuous through the path, without discrete points of the Ho–Cook method being clearly visible. Any delays in the movement should not be present with the generated path, as the entire generated trajectory is continuous.

The generated trajectories were applied within the RobotStudio software in order to illustrate a possible path obtained from the method. Twenty-five trajectory points were inserted into the simulation within the software and the simulation was run. Figure 15 shows the initial, sixth, twelfth, eighteenth, and final (twenty-fourth) trajectory points.

4. Conclusions

In this paper, the GA approach to determining Ho–Cook algorithm parameters was investigated. The results of this investigation point towards the fact that this approach is a valid alternative to analytically determining the Ho–Cook coefficients. This approach can be used in further research, with additional optimization parameters, such as energy efficiency or torque optimization for continuous path planning. Still, it is shown that GA is not capable of generating the final trajectories by itself due to the boundary conditions that need to be set in order to avoid stationary results, and Bezier spline interpolation needs to be utilized to address this. The investigation of the parameters of the algorithm shows that the algorithm has the best performance when the boundary condition is set to 0.05 and the population size is set to 10,000. The number of generations does not seem to influence the convergence point, which is shown to be between 15 and 20 generations, no matter the population size. While the results are satisfactory within the context of the research, they do not particularly improve the current state-of-the-art results in the research [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. Still, as there is a clear lack of the Ho–Cook algorithm being used for planning, the value of the achieved results lies in the fact that it can achieve results similar to the existing ones, indicating that there is room for improvement when more advanced evolutionary or swarm-based optimization algorithms are used. The execution times of the algorithm point out that it cannot be used in online real-time planning as currently presented, with the algorithm only being usable in offline planning, where the trajectories are tuned before their application in manufacturing. Still, in the current form, the algorithm as tested is executed in a single-threaded mode, and executing the algorithm on multiple threads simultaneously could significantly improve the results. Other limitations concerning the shown approach include the need to be familiar with the Ho–Cook path planning algorithm in order to apply the GA developed in this paper, as well as the generated path not being necessarily optimal, as only a near-optimal path is guaranteed by the GA. Path planning in the demonstrated manner has certain natural limitations when compared to other methods such as the dynamic position adjustment of robotic manipulator during the operation. The paths planned in this manner are inflexible to later adjustment, as the algorithm needs to be re-run to obtain the adjusted path. This can cause issues in which only the key changes are made to the paths (e.g., new paths are calculated), as opposed to the continuous tuning of the paths for a higher efficiency (production and energy-wise) [50]. Sometimes, detailed path planning is unnecessary and simply time-consuming compared to the use of simple trajectories generated by online path planning [51]. This manner of online movement training is less skill-intensive compared to detailed offline path-planning, causing an increase in the cost [52]. Finally, fine-tuned offline planned paths such as these are only applicable in predictable environments, with a lack of capability in dynamic planning [53], which would take into account the stochastic nature of the realistic environments. The rigid planning such as that presented can cause a false sense of security, as a static environment is intrinsically assumed, which is rarely the case in real production environments, where issues that may cause faults are rife [54].

Future work in the area should focus on the testing of other evolutionary algorithms on the framework developed for GA, such as differential evolution or particle swarm optimization, to determine whether those algorithms can achieve better results.