Path planning optimization of six-degree-of-freedom robotic manipulators using evolutionary algorithms

Lowering joint torques of a robotic manipulator enables lowering the energy it uses as well as increase in the longevity of the robotic manipulator. This article proposes the use of evolutionary computation algorithms for optimizing the paths of the robotic manipulator with the goal of lowering the joint torques. The robotic manipulator used for optimization is modelled after a realistic six-degree-of-freedom robotic manipulator. Two cases are observed and these are a single robotic manipulator carrying a weight in a point-to-point trajectory and two robotic manipulators cooperating and moving the same weight along a calculated point-to-point trajectory. The article describes the process used for determining the kinematic properties using Denavit–Hartenberg method and the dynamic equations of the robotic manipulator using Lagrange–Euler and Newton–Euler algorithms. Then, the description of used artificial intelligence optimization algorithms is given – genetic algorithm using random and average recombination, simulated annealing using linear and geometric cooling strategy and differential evolution. The methods are compared and the results show that the genetic algorithm provides best results in regard to torque minimization, with differential evolution also providing comparatively good results and simulated annealing giving the comparatively weakest results while providing smoother torque curves.


Introduction
To achieve movement, the torque of robot manipulators actuator (e.g. an electrical motor) needs to be higher than the inertial torsion of the joint the actuator is moving. 1,2 The movement requires use of energy directly proportional to the actuator torque. Successfully minimizing the actuator torque, which is dependent on the joint trajectory, will mean that the movement of the robotic manipulator requires less energy, as well as prolong the work life of the robotic manipulator. 3 It is possible to define the problem of joint torque optimization in a way that makes it possible to use artificial intelligence (AI) methods, namely the evolutionary computing optimization methods, in an attempt to achieve the lower amounts of actuator torque. [4][5][6] Related work Garg and Kumar 7 show that optimization of robot manipulator paths in regard to joint torque is possible using a genetic algorithm (GA) in the case of simple robotic manipulators with two degrees of freedom (DOFs), and the authors propose that similar method can be used for a more complex robotic manipulators. Kazem et al. 4 show the application of a GA on a three-link (redundant) robot arm for point-to-point planning. The paper demonstrated that the GA is an effective method of optimization in that case, capable of achieving multi-objective optimization. Albert et al. 5 have applied GA to optimize joint angles of a threelink planar manipulator system in regard to path control and accuracy in a point-to-point pick-and-place situation. The paper concludes that GA is applicable when dealing with complex path control -giving improved search speed and approximate solution. Sharma et al. 3 attempt to optimize a point-to-point path, interpolated with a fourth-order polynomial, in regard to energy consumed by the robotic manipulator and travelling time with upper limit on the torque. Furthermore, the paper compares the results of a GA with a geometric approach to path planning. The results show a drop in energy consumed when using the path planning with a GA. Wei et al. 8 show the usage of evolutionary algorithms for the development of collective behaviours, specifically task allocation. Authors propose a two-step scheme for task partitioning and allocation to overcome problems presented by the bootstrapping and deception problems. The results obtained by authors show that the proposed approach leads to better performance in comparison with conventional evolutionary robotics approach. Abu-Dakka et al. 9 evaluate the efficiency of the time and trajectory smoothness optimization by a multiple population GA using analysis of variance test. Authors conclude that the approach of using a multiple population GA is valid and produces results quickly. Larsen et al. 10 show the application of evolutionary algorithms and rapidly exploring random trees (RRT) with various parameters on a robot system with KUKA KR210 and KUKA R3100 robots on a common axis to calculate collision-free paths. The results show superior results when using evolutionary algorithms in comparison with sampling methods, as well as concluding that the results gained using GA are superior to ones gained with RRT. Wei et al. 11 use deep Q-learning AI algorithm to develop end-to-end control policies for robotic swarms. It is shown that use of such techniques enables development of control policies using fewer computation resources, especially in cases with large solution spaces. El Haiek et al. 12 show the trajectory planning of a three-DOF spherical robot with evolutionary algorithms in regard to the minimization of travel time and joint travelling distance. The path is point-to-point and interpolated using polynomial interpolation. Based on the results, by using GA the robot can successfully find an optimal, collisionfree, trajectory while minimizing the time and the length of the trajectory.
This article observes the usage of three different algorithms GA (in two different configurations) which is an algorithm based on natural evolution, 13,14 simulated annealing (SA) (in two different configurations) which is an algorithm based on the metal annealing process 15,16 and differential evolution (DE), an algorithm inspired by the differential equations. 17,18 The above algorithms are heuristic space search algorithm, designed to search the solution space in an attempt to find the optimal solution. Unlike the deterministic methods, like an extensive search of the solution space, the algorithms used do not guarantee the solution to be the optimal solution, but they often provide good solutions significantly faster in comparison with deterministic methods. 19,20 This article is based on research done in 'Optimization techniques applied to multiple manipulators for path planning and torque minimization' by Garg and Kumar,7 and attempts to apply the modified techniques on to a more complex case -using a robotic manipulator with six DOFs, as opposed to a robotic manipulator with two DOFs.

Hypothesis
To summarize the novelty of this article, the idea is to investigate the possibility of evolutionary computing algorithms implementation for optimization of robot manipulator trajectory with respect to minimization of joint torques in cases of complex robotic manipulators and performance comparison of multiple AI algorithms on solving trajectory planning of complex robot manipulators with minimization of joint torques.
This article compares the results and execution times of aforementioned algorithms and their configurations in two cases -a single robotic manipulator and two cooperating robotic manipulators, both for point-to-point path planning.

Research methodology
In this section, the methodology of the research article will be presented. First, the cases studied in this article will be described, followed with the description of algorithms used to calculate the kinematic and dynamic robotic manipulator equations, followed by the description of algorithms used for optimization.

Cases
This article observes two cases: a single robot manipulator with six DOFs ( Figure 1) and two such manipulators working cooperatively by following the same path ( Figure 2). Cooperating robotic manipulators are often used in industrial applications. 10,21 In both cases, the method is validated using ABB IRB 120 robotic manipulator. In observed cases, robotic manipulators are transporting a prismatic shaped object with weight of 2 kg along the calculated path.
For path planning, all joints of robotic manipulators start the movement in a position of 0 radian, and finish the movement in 1 radian. The robotic manipulators are stationary at the beginning and the end of the movement making their starting linear and angular velocities as well as the linear and angular accelerations equal zero. 7,22 Robot manipulator dynamics To determine the dynamic equations, two different algorithms are used -Lagrange-Euler (L-E) and Newton-Euler (N-E). Both algorithms should produce the same results and can be used to check the validity of the calculated results. [23][24][25] First step in determining the dynamics of the robotic manipulators is calculating the kinematic matrix. This is also needed for determining the inverse kinematic equations which are necessary for determining the path of the second robotic manipulator in the case with two cooperating robotic manipulators. 26 Kinematics. Kinematics are determined using the Denavit-Hartenberg (D-H) method. The D-H method sets the orthonormal coordinate systems in each joint of the robotic manipulator, in such way that axis z k (where k is a counter determining the joint) matches the axis of movement of the joint k. 27 Once all the coordinate systems are positioned, the parameters q k , d k , a k and a k are determined and they are placed in the homogeneous transformation matrix for the joint 28-30 Because of the size of the equations, the shortened trigonometric format is used. When using this format, the trigonometric functions are written using only the first letter of their name, so sine is written as S and cosine as C. In addition, the arguments of the functions are written as indexes -so sin a 1 ð Þ becomes S a 1 . If all the arguments are marked with the same letter, for example q 1 and q 2 , then only the indexes of the variables can be used, for example sin q1 ð Þ becomes S 1 and cos q2 ð Þ becomes C 2 . If the argument consists of multiple variables, then it can be written as sin q1 þ q2 ð Þ!S 12 or cos q3 À q4 ð Þ!C 3À4 . The homogeneous transformation matrix of the robot manipulator is calculated as a product of matrices of each joint using The resultant matrix is formatted as where R q ð Þ ¼ r 1 r 2 r 3 ½ is tool orientation matrix (3 Â 3), with r 1 being the perpendicular vector, r 2 being movement vector and r 3 being the approach vector; p(q) is the tool end position, v T 1 being perspective vector (usually 0 0 0 ½ ) 31 and s being scaling coefficient (usually 1). 32,33 The inverse kinematics equations can be extracted from the homogeneous transformation matrix, where the  position of the tool is defined with the vector w, with w defined with 34 where x, y and z are tool positions in the tool configuration space and F, Y and C are Euler angles (spin, nutation and precession). 35,36 Vectors w 1 and w 2 are defined from homogeneous transformation matrix given in equation (3) with equations shown in 37 and L-E algorithm. L-E algorithm is an iterative method of determining the differential equations defining the joint torque of the robotic manipulator. 22,38 The L-E algorithm defines the differential equations per 26,39 where i represents the robot manipulator link, P n j¼1 D ij q ð Þq j h i represents the inertial forces and moments, P n k¼1 P n j¼1 C i kj q ð Þq k q j h i represents Coriolis forces, h i q ð Þ is the effect of the gravity on the robotic manipulator, b i q ð Þ is the friction within the robotic manipulator and t i is the moment of the actuator. 40 Before the start of L-E algorithm, additional starting values need to be set. Those values are and i ¼ 1; where i is the counter that determines which joint of the robotic manipulator are the values being calculated for in the given step, where i is 1 for the joint next to the base of the robotic manipulator. The first step in L-E algorithm is the determination of homogeneous coordinates of the link the calculation is being done for Dc i , as well as the tensor of inertia D i q ð Þ using The next step is calculating the vector z iÀ1 q ð Þ as Then the matrix of a complex homogeneous transformation is calculated as a product of transformation matrices for all the joints between the base and the current link as The matrix T i 0 allows for calculating the position of the centre of mass coordinates for the link i in relation to the coordinate system of the base of the robotic manipulator per where H 1 is defined as The tensor of inertia in relation to the base coordinate system is also determined as Final step in the first part of the L-E algorithm is calculating the Jacobian matrix. The Jacobian matrix connects the infinitesimal movements of the manipulator joints to the infinitesimal movements of the tool and is defined using where matrices A and B are respectively defined with and With these matrices calculated, the manipulator tensor of inertia can be determined as When the manipulator's tensor of inertia has been calculated, the value of the counter i is increased by one. If i is smaller or equal to n, where n is the number of joints of the robotic manipulator, the algorithm returns to the beginning of the process. If i is greater than n, the L-E algorithm proceeds with the second iterative part.
The beginning of the second iterative part is resetting the value of i to 1. Then, for the joint given by i, the speed connectivity matrix is calculated as The gravity influence vector is calculated as Furthermore, the friction is determined as With this, all the necessary values are calculated and the L-E equation for the given link can be determined using The counter i is increased. If the value is smaller or equal to n, the algorithm returns to the first step of second iterative part and the calculation is repeated for the given link. If the value of i is greater than n, the L-E algorithm is finished.
N-E algorithm. N-E algorithm is a recursive method used to calculate the dynamic differential equations. It has two parts: the calculation 'forward' in which the linear and angular speeds and accelerations are calculated, and the calculation 'backward' in which the forces and momenta acting upon the link are calculated. 41 In the calculation 'forward', the values are calculated for the links in the direction from the base of the robotic manipulator towards the tool, while in the calculation 'backward' the values are calculated from the direction from the tool to the base. 22,42,43 This is illustrated in Figure 3.
Before starting the calculation, the starting values need to be set to T 0 When these values are set, the algorithm progresses into its first part -calculation 'forward'.
The first part of the calculation 'forward' is the same as the start of the first iterative part of L-E algorithm in which the vector z iÀ1 is determined as The next step is the calculation of the angular speed as where x i is defined as Then the value of angular acceleration is calculated per As in the L-E algorithm, the complex homogeneous transformation matrix is determined followed by the vector Ds as Final part of calculation forward is calculating the linear acceleration using After the linear acceleration is calculated, the counter i is increased by one. If the value of i is lesser than n, the calculation repeats for the next link. If the value is greater than n, the algorithm proceeds to the calculation 'backward'. Calculation 'backward' is started by calculating the vector r k and the tensor of inertia of the robotic manipulator link given by i in relation to the base of robotic manipulator Then the forces and the momentum acting upon the link given by i are calculated as and The final calculation is determining the value of the joint actuator momentum using the equation The value of i is then lowered by one. If the value of i is larger than zero, the calculation 'backward' is repeated for the next joint. If the value of i is zero, the base of the robotic manipulator has been reached and the N-E algorithm is finished.

Used algorithms
In this article, five different algorithms were utilized and these are: GA with average recombination, GA with random recombination, SA with linear cooling strategy, SA with geometric cooling strategy and DE.
Used algorithms will be discussed in this chapter. All of the above are evolutional computing algorithms, which are based on performing the iterative method on a certain population of agents to achieve an optimal result from the solution space. Solution space is defined as the ndimensional mathematical space containing all the possible solutions to a given problem, 45,46 with each solution having a value calculated with a fitness function defining how well does the given solution satisfy the conditions given for solving the problem. 47,48 Each agent is a point in the solution space -an instance of a solution containing genotype data. Genotype is the simplified representation of the solution to the given problem which is the phenotype of the problem. 49 In this article, phenotype is the movement of the robotic manipulator through the space, while the genotype is the parameters of equations that describe the movement of the robotic manipulator. A group of agents upon which the optimization is performed is called the population 50,51 and each iteration of the evolutionary algorithm is called a generation. 52 All algorithms end when an end condition 53,54 is fulfilled. End condition is usually fulfilled after a certain amount of generations have passed or when there was not a change in the best achieved value of fitness function for a certain amount of generations. 55 All algorithms are run with population of 20 agents for the period of 50 generations.
Genetic algorithm. GA is an optimization algorithm that imitates the process of biological evolution using the mechanisms of crossover and mutation. 55,56 Crossover is a mechanism in which the new agent (child) is created by randomly choosing two existing agents from the population (parents) and generating the values of parameters of the child agent from the values of parameters of parent agents. The idea of crossover is that by generating new agents, and consistently keeping the better fitting agents in the population, through the recombination of agents, the optimal solution will be found. 56 The wish is for the crossover to result in the better fitting gene which will replace a worse fitting agent currently in the population. If the new agent is not of a higher quality, then the current gene is kept. 57,58 There are multiple ways to cross the agent's parameters to generate a new agent. In this article, two mechanisms are used: random recombination and average recombination. 59,60 When using the average recombination, the values of child agent's (g c ) parameters (a i ) are calculated as the average value of the parameter of the parent agents (g A and g B ). This is given as : : a n ¼ g A Á a n þ g B Á a n 2 In the case of random recombination, the parameter values of the child agent are randomly chosen between the appropriate parameter values of the two parent agents, as : : a n \ g A Á a n ; g B Á a n f g Mutation is the process in which an agent's parameters are randomly chosen from the solution space and the agent created in that way replaces an agent in the population, with no attention paid to the fitness function values -so a worse mutated agent can replace the better one. 45 The mutation mechanism is necessary to avoid the stagnation that can happen when only the crossover mechanism is implemented. It allows the GA to move away from the situation where all the agents are in the local optimum of solution space and give it ability to find the global optimum. 61 Mutation is used in a low amount of cases (1% of the iterations), as opposed to the crossover which is used comparatively often (80% of the iterations). 62,63 The described procedure of GA is shown in Figure 4.
Simulated annealing. SA is an optimization algorithm inspired by the process of metal annealing. The SA algorithm has a system temperature which lowers, according to the cooling strategy, through the algorithm iterations. 64,65 The algorithm selects a solution from the neighbouring solution space to the current solution and, if that solution is better, replaces the current solution with it. The system temperature determines the likelihood of the algorithm accepting the worse solution. The likelihood of the worse solution being chosen is proportional to the system temperature. This is usually implemented by generating a random value in the possible temperature range and accepting a worse solution if that value is lower than the current temperature. This mechanism serves a similar purpose as the mutation mechanism does in the GA. This gives the SA algorithm the capability to select a different solution if the current solution is located in the local minimum of the solution space. 66,67 The likelihood of choosing a worse solution at some point of the algorithm execution can be configured by choosing a cooling strategy. 68 Two strategies are used in this article -linear as shown and geometric strategy as where k is the current algorithm iteration, t k is the current system temperature, t 0 is the starting system temperature and a and b are cooling factors which determine the speed at which the system cools. 69 The process which SA algorithm uses to find the optimal solution is shown in Figure 5.
In this article, the starting temperature of the SA algorithm is set to 100, with the coefficient a for linear cooling strategy being set to 1, and coefficient b for geometric cooling strategy being set to 0.5.
Differential evolution. DE is an AI optimization algorithm, which uses differential equations as the inspiration. This algorithm is particularly fitting for the solution spaces that are not smooth and have interruptions. 70,71 DE algorithm creates a new agent by combining three agents using : : : a n ¼ g a Á a n þ FÃ g b Á a n À g c Á a n ð Þ where F is a constant that determines how distant the new solution will be in the solution space from a randomly selected agent g a , and has a value in range F 2 h0; 2. 72,73 Three agents are randomly picked from the population and the values of new agent are calculated. If the value of the new agent is better than the agent in the current iteration, it replaces the current agent. 70 Unlike SA algorithm, DE does not have a mechanism which allows it to select a worse solution, nor does it have a mechanism similar to mutation such as can be seen in GA. DE avoids stagnation by having a large search area (configured by parameter F). 74 A larger value of the parameter F will give a larger search area, but will increase the time it takes for the algorithm to converge to a solution, while the smaller value will cause a faster convergence to a solution, but it will be more likely the found solution will be a local instead of a global optimum. In this research, parameter F is set to 1.2, which is the starting value used in most research. 71 The process of DE algorithm is shown in Figure 6.

Agent construction
In the cases observed in this article, the phenotype of the robotic manipulator movement is represented with a parameter of the equation describing that movement as [75][76][77] Parameters b, c, d, and e can be removed from the equation if the following conditions are observed 22,77 : the starting time is zero (t s ¼ 0s) and the movement lasts 2 s (t f ¼ 2s) the robotic manipulator joints move from 0 radian to 1 radian (q 0 ð Þ ¼ 0 rad, q 2 ð Þ ¼ 1 rad), it starts movement as stationary -the starting speeds and accelerations are zero (v t s ð Þ ¼ 0 m/s, dv dt t s ð Þ ¼ 0 m/s 2 , ! t s ð Þ ¼ 0 rad/s, d! dt t s ð Þ ¼ 0 rad/s 2 ) and it ends the movement stationarythe ending speeds and accelerations are equal to zero The parameter e is eliminated using the condition that the starting joint angle equals 0 Parameter d can be eliminated by the derivation of equation (41) to get the equation describing the angular joint velocity and using the starting condition of initial angular joint velocity being 0 Second derivation of equation (41) is the expression for the joint's angular acceleration is derived. If the condition that the initial joint's angular acceleration is zero is used, the parameter c is eliminated per Parameter b can be defined using By inserting the starting and ending conditions, parameter b can be eliminated per making it equal to a constant. With these parameters eliminated, the only remaining parameter is parameter a which will describe the movement of each joint of the robotic manipulator.
In the case of the single robotic manipulator, the joint movements are described by six parameters a i , where i is the joint number. Since this is the parameter used for optimization, it can be referred to as optimization parameter. In the case of two robotic manipulators, the joint movements need to be calculated using the inverse kinematic equations. The values of vector w are calculated according to equations (4) to (6). Vector w contains the equations containing values of all six joints. The joint variables (q 1 , q 2 , q 3 , q 4 , q 5 and q 6 ) can be extracted from these equations which gives a set of six inverse kinematic equations that define the joint values as function of the end-effector position in space (vector w). [78][79][80] By entering the values of vector w into these inverse kinematic equations, it is possible to get the values of joint variables needed to achieve the position described by vector w. Hence, if the joint variable values of the first robotic manipulator (q 1 1 , q 1 2 , q 1 3 , q 1 4 , q 1 5 and q 1 6 ) are determined from the optimization parameters a i , the direct kinematic equations (5) and (6) can be used to determine the position of endeffector in the space, defined by vector w. 21,38 Using the inverse kinematic equations to determine the joint variables of second robotic manipulator (q 2 1 w ð Þ, q 2 2 w ð Þ, q 2 3 w ð Þ, q 2 4 w ð Þ, q 2 5 w ð Þ and q 2 6 w ð Þ) will result in the endeffector position of the second robotic manipulator following the positioning of the first robotic manipulator, enabling cooperative behaviour. 39,76,78 This process is shown in Figure 7.
That means that in both observed cases generating only a single set of parameters a i (one for each joint) is necessary, as those will be used for determining the movement of the first robot and the second robotic manipulator's movement will be determined by the movement of the first, which gives the genotype shape defined by where a i is part of the interval a i 2 À10; 10 ½ ð 48Þ

Fitness function
To be able to determine how well does each agent solve a problem, it is necessary to define a fitness function which will provide a value describing it. In this instance, the fitness function is defined as the sum of the total torsion on each point of the trajectory, where the total joint torque is defined as the sum of joint torques on each joint of the robotic manipulator(s). 7 This is defined by the equation where M is the number of points in the trajectory, n is the total number of joints of the robotic manipulator and t i is torque of the ith joint. Squaring each torque and taking the square root of their sum in each step of the trajectory provides the absolute value of all joint torques acting on the robotic manipulator.
Since the observed cases in this article have 20 points in their trajectory and n is either 6 (in the case of a single robotic manipulator) or 12 (in the case of the two cooperating robotic manipulators), the fitness function for each of the two cases can be defined for the case with the single robotic manipulator as and for the case with two cooperating robotic manipulators as

Results
The result of the D-H method can be written in the following form where T 31 ¼ 0 ð60Þ T 32 ¼ 0 ð64Þ In the results, the algorithms used are given shortened names: GA-A represents the GA with average recombination, GA-R represents the GA with random recombination, SA-L represents the SA algorithm with linear cooling strategy and SA-G represents SA with geometric cooling strategy.
The comparison of results between algorithms is given in Table 1. The best results and shortest execution times are shown in bold in the table for emphasis. In Table 1, D represents the average over multiple runs of the difference between the best fitness value from the initial randomly picked population and the best fitness value from the final optimized population, D ¼ f g p 1 À f g p 50 .
In Figure 8, a drop in the fitness function value is shown, which was calculated using equation (49), as is expectedconsidering this is an optimization problem which has the goal of lowering the value of fitness function which represents the total joint torques of robotic manipulators. The fitness shown is the fitness function value of the best agent in the population of 20 agents in the given generation. The initial drop to a better solution happens in the first few generations, with a further drop down before the 30th generation. The last drop of the fitness function value happens around the 35th generation and then it stays at the lowest optimized value. Figure 8 shows a significant drop in the fitness function value, most probably due to a rather suboptimal initial randomly selected solution.
In Figure 9, the joint trajectory after the optimization process is shown. It can be seen that the joints start at the position of 0 radian and end in the position of 1 radian. The calculated joint trajectories are smooth, without any sudden jumps and changes that would cause large differences in joint torques or generate sudden movements of the robotic manipulator. Figures 10 and 11 show the resultant trajectories of joints of robotic manipulators in the case with two cooperating robotic manipulators. The joint trajectories for the case of two cooperating robotic manipulators have been shown on two separate figures for easier viewing. Trajectories of both robotic manipulators show smooth curves. Figure 12 shows the change of joint torque on each joint of the robotic manipulator. It can be seen that in the optimal case found in this particular optimization run, the joint torque shows a drop into negative values followed by a smooth raise. The smooth changes in the joint torques like these have shown to be common when using SA algorithms. The joint torque of the second joint exhibits the highest rise, caused by the robotic manipulator configuration where the second joint is placed under most stress. Figures 13 and 14 show the change of joint torques on the first and second robotic manipulator. As with the joint trajectories for dual manipulators, the joint torques have been split into two figures for easier viewing. Some sudden changes can be seen in the direction of the joint torques (observe 14 at 0.5 s). This is prominent with the use of GA and DE for the optimization process. The sudden changes GA-A: genetic algorithm with average recombination; GA-R: genetic algorithm with random recombination; SA-L: simulated annealing with linear cooling strategy; SA-G: simulated annealing with geometric cooling strategy; DE: differential evolution.   like these are also especially prominent on the second robotic manipulator for which the joint angles are determined using inverse kinematics as opposed to the first manipulator for which the joint angles are determined using equation (41). On both robotic manipulators, joint torque of the second robotic manipulator joint shows the highest torque, which is caused by the robotic manipulator configuration.

Discussion
The results show a successful decrease in the total joint torque in the robotic manipulator in comparison with the initial randomly selected parameters. The best results and the fastest execution time are shown by the GA using average recombination for the dual cooperating manipulator and random recombination for the case with the single manipulator. The DE shows good results in both cases, having better results than the average recombination GA in the single manipulator case and better results than the random recombination algorithm in the dual manipulator case. DE shows fast execution times in both cases, comparable with the GA. SA shows mixed results, with the linear cooling strategy showing the weakest results in both cases. The SA with geometric cooling strategy shows results comparable with the average recombination GA (although significantly weaker compared to the random recombination GA) in the single manipulator case, but it shows the weakest results in the dual manipulator case. SA shows the longest execution time necessary to complete optimization, with the linear cooling strategy being the longest running algorithm of all compared in this article. The results consistently show that the largest joint torque is the torque of the second joint of the robotic manipulator (t 2 ). These large torque values are caused due to the configuration of the robotic manipulator where the second joint has most stress placed on it during the movement of the robotic manipulator. This can be observed when optimizing with all used algorithms, and is especially pronounced when optimizing with DE. Large values of joint torques near the end of the trajectory that sometimes happen are caused by the fact that the end conditions guarantee the robotic manipulator ending its movement in the defined position, making it harder to optimize the trajectory near the end of    the movement. It is important to note that, despite providing the weakest improvements in regard to the amount of total joint torque lowered, the SA (in both configurations) is the only algorithm that provides results with the smooth torque and joint curves in almost all instances. The curves provided by the GA have a tendency to show sudden changes in joint torque, which can negatively impact the durability of the robotic manipulator, while the ones provided by the DE have a tendency to have one of the joints of the robotic manipulator have a high torque, with very low joint torques on the other joints. While this does result in a low total joint torque, it might not be appropriate depending on the robotic manipulator used. This points towards the SA being the best algorithm for producing a smooth joint torque curves after the optimization.

Conclusion
The results show a successful optimization in both observed cases using evolutionary algorithms. The best results among the used algorithms are provided using GA with random recombination for single robotic manipulator and with average recombination for dual cooperating manipulators. DE has also shown good results in comparison with other algorithms. For the objectives optimized in this article, SA shows the worst results. But, the generated results are not optimal in regard to the smoothness of the curves, where SA shows the best results in comparison with other algorithms used, and stress placed on second joint of robotic manipulator(s). This article shows successful application of evolutionary algorithms on the problem of path optimization of the realistic six-DOF robotic manipulator, where earlier research had only shown application on simpler manipulators. In practice, this means that by applying evolutionary algorithms to perform optimization on the planned paths, the amount of energy used for the robotic manipulators to perform work could be lowered, as well as the longevity of the robotic manipulators raised due to less damages caused by the torsion developed in the joints during the movement. The work done in the article does come with certain limitations. First is that the dynamic model used is not general. While this means the optimization results are high quality for the used robotic manipulator, it also means that to apply the algorithms described on a different robotic manipulator entire dynamics need to be recalculated to be used for optimization. Additionally, this article describes methods as used on the point-to-point planning, and do not apply to widely used continuous path planning -which would require the way joints values are calculated to be adjusted. Future work could concentrate on multi-objective optimization that would attempt to lower the joint torques as one objective, and provide a smooth curve as the other and more equal joint torque distribution. Additionally, future work could concentrate on wider parameter variations, friction implementation and continuous path planning. Instead of using the recommended algorithm parameter values, such as they were in this article, performance of algorithms could be tested with different variations of their parameters in an attempt to increase the performance further. Friction, as mentioned, was not included in the model used in this article, because of complexity of determining the appropriate friction model. Friction in the manipulator could be modelled using different approaches, and the influence of them on the optimization process measured. Finally, while point-to-point path planning has wide applications, modifying the algorithm to perform optimization over continuous paths could provide wider application for this research.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research has been (partly) supported by the CEEPUS network CIII-HR-0108,