Time-Jerk optimal Trajectory Planning of Industrial Robots Based on a Hybrid WOA-GA Algorithm

: An optimal and smooth trajectory for industrial robots has a positive impact on reducing the execution time in an operation and the vibration in their joints. In this paper, a methodology for the time-optimal and jerk-continuous trajectory planning of industrial robots is proposed. The entire trajectory is interpolated in the joint space utilizing ﬁfth-order B-splines and then optimized by a hybrid whale optimization algorithm and genetic algorithm (WOA-GA). Two objective functions, including the integral of the squared jerk along the entire trajectory and the total execution time, are minimized to obtain the optimal entire trajectory. A ﬁfth-order B-spline interpolation technique enables the achievement of a jerk-continuous trajectory, while respecting the kinematic limits of jerk, acceleration and velocity. WOA-GA is utilized to solve the time-jerk optimal trajectory planning problem with nonlinear constraints. The proposed hybrid optimization algorithm yielded good results and achieved the time-jerk optimal trajectory better under kinematic constraints compared to the genetic algorithm, whale optimization algorithm, improved whale optimization algorithm with particle swarm optimization and adaptive cuckoo search algorithm. The numerical results show the competent performances of the proposed methodology to generate trajectories with high smooth curves and short total execution time.


Introduction
Due to their outstanding performance in accuracy, repeatability and efficiency, industrial robots have become one of the most important equipment to replace human resources in the machinery, automobile, electronics, metal, and other industries [1][2][3].Generally speaking, an industrial robot is a nonlinear dynamic mechanism, which needs to follow a specific trajectory to avoid any possible collision.Therefore, trajectory planning is an essential technique in automated industrial activities using industrial robots.A reasonable trajectory can reduce the vibration and noise of robots and ensure a high trajectory tracking accuracy and motion stability of robots [4][5][6].
In the relevant literature, trajectories are mainly constructed by polynomial, spline, Bezier, NURBS and other interpolation functions.In order to obtain a smooth trajectory with jerk boundary, Jond et al. used polynomials of the orders 3, 4 and 5 to generate an acceleration-continuous trajectory, in which the acceleration values at the start and end motions are non-zero [7].Sonja Macfarlane et al. employed a series of fifth-order polynomials to plan a smooth trajectory between two directional points [8].Similarly, Bureerat et al. used quintic polynomials to connect the initial and intermediate points, to finally the end point [9].Lin et al. proposed a method of a 7-segment acceleration profile motion model based on a time-optimal trajectory, which divided all waypoints into continuous batches, and each batch had two overlaps of waypoints [10].Wu et al. proposed a novel point-to-point trajectory planning algorithm (PTPA) to improve the motion efficiency of industrial robots [11].Moreover, spline functions have been used as trajectory primitives to ensure acceleration continuity or jerk continuity.Since each control point on the B-spline curve will only affect the parameters within a range of the curve, various B-splines were used to construct joint trajectory.For example, the minimum time trajectory is generated by cubic splines for the robotic manipulator [12][13][14].Muhammad Zeeshan et al. interpolated the entire fifth-order trajectory with B-spline interpolation method to obtain the minimum time trajectory [15,16].
Trajectory optimization is a highly complicated problem due to the nonlinear characteristics of trajectories and the need of motion efficiency, smooth operation and low noise of industrial robots.Many studies take the minimum time [17][18][19][20][21], the minimum acceleration [9,16] and the minimum energy [22,23] as the trajectory planning optimization objectives.In addition, the hybrid optimization methods have become the focus of robot research in recent years.For example, Zhang et al. proposed a multi-objective optimization method for robot trajectory planning under obstacles, including time optimization and acceleration optimization.This method uses quintic B-splines to generate trajectory configurations in the joint space with motion constraints and obstacles constraints [24].Wang et al. proposed a new algorithm to define the objective function from two aspects of time and speed [14].Zhang et al. proposed an adaptive cuckoo search (ACS) algorithm to optimize the trajectory of quintic B-spline programming [25].Ameer adopted an RNN based metaheuristic approach (RNN-MA) to achieve real-time implementation [26].Cheng proposed beetle antennae search algorithm (BAS) to tackle the formulated optimization function [27].Hybrid algorithms, such as an improved whale optimization algorithm with differential evolution algorithm or particle swarm optimization (IWOA-DE and IWOA-PSO), and new genetic algorithm and adaptive elite genetic algorithm (AEGA-SA) have been proposed to improve the exploration and exploitation capacities of the original algorithm [21,28,29].Several optimization methods with their characteristics are shown in Table 1.
Table 1.Several optimization methods with their characteristics.

IWOA-DE [21]
Good global optimization ability, strong robustness and slow convergence speed.
ACS [25] Faster convergence degree, higher accuracy and better global search ability, but the final accuracy needs to be improved.RNN-MA [26] Fast convergence speed and suitable for real-time implementation.
BAS [27] Fast convergence speed with low searching precision and convergence speed, and sometimes it is prone to local optimum.
AEGA-SA [28] Superior to the original genetic algorithm in terms of trajectory smoothness and trajectory efficiency, but it is highly dependent on population size, and is prone to low convergence speed or fall into local optimum.IWOA-PSO [29] Good global exploration ability and slow convergence speed.
This paper proposes a time-optimal and jerk-continuous trajectory planning technique for industrial robots.This technique is based on fifth-order B-spline and hybrid WOA-GA.Fifth-order B-spline is utilized to interpolate the multi-point trajectory.The WOA-GA combining the advantage of GA and WOA is proposed to solve the time-jerk trajectory optimization problem with nonlinear constraints.The rest of this paper is organized as follows.Section 2 defines the problems of trajectory optimization.In Section 3, the fifth-order B-spline is chosen to interpolate the joint trajectory, and the hybrid WOA-GA is presented.In Section 4, we compare the simulation results from different optimization algorithms, and Section 5 is the conclusion.

Problem Formulation
Trajectory planning is mostly carried out in the joint space, which can ensure the stability of joint motion and avoid motion jump of the manipulator, while planning trajectories in the operation space.Thus, this paper focus on planning trajectories in the joint space.
We define the optimal planning problem similarly to Ref. [12].Two objective functions, including the integral of the squared jerk along the whole trajectory and the total execution time, are minimized to obtain the optimal whole trajectory.These two terms have opposite effects.This implies that reducing the total execution time will lead to trajectories with large values of the kinematic quantities, while reducing the integral of the squared jerk along the whole trajectory will lead to a smoother trajectory with low efficiency.The trade-off between these two objectives can be achieved by legitimately adjusting the weights of the two objectives.Therefore, the optimization problem can be formulated as: The meaning of the symbols in Equation ( 1) is described in Table 2. Velocity limit for the jth joint (symmetrical) AL j Acceleration limit for the jth joint (symmetrical) JL j Jerk limit for the jth joint (symmetrical)

Multiple B-Splines
Multiple B-spline curves are widely used in the trajectory planning algorithms of robots with good local support and continuous smoothness.A B-spline curve can be expressed mathematically as follows: ( where p(u) is the joint position, d i is the ith control points of the B-spline, k is the degree of the B-spline, N i,k (u) is the ith basis function of the kth order B-spline, n is the total number of the control points and u is the normalized knot vector variable.N i,k (u) can be obtained using the Cox-de Boor recursion formula.
The trajectory of an industrial robot is discretized to some designed points on the time scale.Two virtual points are adopted at the second and the second-last positions of the designed point sequence, due to the consideration of zero jerk at the initial and final movements.Normalized knot vector u is generated from the sequence of t corresponding to the designed points [16]. where The lth order derivatives can be described as where

Whale Optimization Algorithm
Whale Optimization Algorithm (WOA) is a recently proposed population-based optimization algorithm inspired by the hunting behavior of humpback whales [30].The WOA algorithm starts with a set of random solutions, including three kinds of hunting behaviors: encircling prey, exploration and exploitation.At each iteration, the location of individuals in the population is updated according to the most promising solution [31].

Encircling Prey
Since the prey in the search space is not known, the WOA assumes that the current best-so-far solution is the target prey.Then, the other individuals will try to update their positions towards the target prey.This behavior is represented by the following equations: where A and C are coefficient vectors, X best is the current best-so-far solution, X is an individual, g is the current iteration and | | denotes the absolute value.
The vectors A and C are calculated as follows.
A = (2r − 1)a ( 9) where a is linearly decreased from 2 to 0 over the iterations and r is a random set in [0, 1].The selection of this operator depends on the coefficient vector A and a random number Po, where Po is a random number in [0, 1].For each individual, if Po < 0.5 and |A| < 1, the position is updated by encircling prey.

Bubble-Net Attacking Phase
If Po > 0.5, the bubble-net predation stage begins.The mathematical model of the search is as follows.
where b is a constant for defining the shape of the logarithmic spiral and m is a random number in [−1, 1].

Searching for Prey
If Po < 0.5 and |A| ≥ 1, the individuals need to expand their exploration areas.This phase uses a randomly generated set of solutions X rand to update the individuals' locations, which allow the WOA algorithm to perform a global search.The rules of the update solution are as follows.

Genetic Algorithm
Genetic algorithm (GA) is a random global search and optimization method developed by imitating the biological evolution mechanism of nature [32].It can automatically obtain and accumulate knowledge about the search space in the search process, and adaptively control the search process to obtain the best solution.Each individual in the population represents a possible solution to the optimization problem.The adaptation ability of individual is judged by the fitness function.The individuals with poor fitness are eliminated, and the individuals with good fitness can continue to reproduce.In the process of reproduction, it needs selection, crossover and variation to form a new population.Then, the better solution will finally be obtained.

Proposed Hybrid WOA-GA Algorithm
The flowchart of the proposed WOA-GA algorithm is shown in Figure 1, and detailed introductions are as follows: Step 1: Define the initial parameters, such as Pop_quantity, MaxGen, Pc, Pm and Er.Then, an initial population is obtained randomly.
Step 2: Calculate the fitness value of each individual using Equation (1) and find X global_best and X local_best .
Step 4: Update the binary-encoded population by GA.Roulette wheel selection method, single-point crossover and bitwise mutation are utilized in the corresponding operators.
Step 5: Calculate the fitness value of each individual and update X global_best and X local_best .
Step 6: If t < MaxGen, then t = t + 1, update a, A, C, m, b and Po, and return to step 3.If t = MaxGen, then the best fitness and X global_best are obtained.
Since finding a suitable constraint handling method for WOA-GA is out of the scope of this work, a constraint handling technique named static penalty is employed [33].
Step 6: If t < MaxGen, then t = t + 1, update a, A, C, m, b and Po, and return to step 3.If t = MaxGen, then the best fitness and X global_best are obtained.
Since finding a suitable constraint handling method for WOA-GA is out of the scope of this work, a constraint handling technique named static penalty is employed [33].

Results and Discussion
In order to verify the feasibility of the proposed method, the numerical simulations were tested by solving six unconstrainted optimization problems and a time-optimal trajectory planning problem using 5th-order B-splines.WOA-GA was programed in MATLAB R2018b, as well as WOA, GA, IWOA-PSO [29] and ACS [25], which were utilized for comparative analysis.The configurations of the personal computer were a Core i5-1035G1 CPU with 16 GB of RAM.The parameter configuration of each algorithm are presented in Table 3.All the algorithms were initialized with the same population for each run.For WOA-GA and GA, we utilized 14 bits to code each time interval.The best-so-far feasible solution in the population obtained at the end of generations was utilized to evaluate WOA-GA on effectiveness and stability.The mean ('MEAN') and standard deviations ('STD') of the best-so-far solutions are utilized as evaluation indices, formulated as follows:

Results and Discussion
In order to verify the feasibility of the proposed method, the numerical simulations were tested by solving six unconstrainted optimization problems and a time-optimal trajectory planning problem using 5th-order B-splines.WOA-GA was programed in MATLAB R2018b, as well as WOA, GA, IWOA-PSO [29] and ACS [25], which were utilized for comparative analysis.The configurations of the personal computer were a Core i5-1035G1 CPU with 16 GB of RAM.The parameter configuration of each algorithm are presented in Table 3.All the algorithms were initialized with the same population for each run.For WOA-GA and GA, we utilized 14 bits to code each time interval.The best-so-far feasible solution in the population obtained at the end of generations was utilized to evaluate WOA-GA on effectiveness and stability.The mean ('MEAN') and standard deviations ('STD') of the best-so-far solutions are utilized as evaluation indices, formulated as follows: where f i best is the best-so-far feasible solution obtained in the ith run and Nr is the number of runs.Note that the smaller the values of MEAN and STD, the more reliable and stable the solutions that the algorithm can provide [34].

Application 1: WOA-GA for Test Functions
The first six unconstrainted optimization problems are classical benchmark functions shown in Table 4.These functions include unimodal (F 1 -F 2 ), variable-dimension multimodal (F 3 -F 4 ) and fixed-dimension multimodal (F 5 -F 6 ) benchmarks.Table 4 shows the details of the test functions.Note that only one optimal solution exists in the unimodal functions to evaluate the convergence rate and exploitation capacity of the algorithm, and a global optimal solution with several local ones exist in the multimodal functions to evaluate the exploration capacity of the algorithm.For each test function, WOA-GA was performed 20 times as well as WOA, GA, IWOA-PSO and ACS.For all algorithm, the test population size and the maximum iteration were set to 30 and 500, respectively.
u(x i , 10, 100, 4) Table 5 reports the statistical results for the unimodal function and shows that WOA-GA can obtain the global optimal solution for F 1 and F 2 .Table 6 reports the statistical results for the variable-dimension multimodal and fixed-dimension multimodal functions.WOA-GA can obtain the near-optimal solution for F 3 and F 4 , and the global optimal solution for F 5 and F 6 .Tables 5 and 6 indicate that WOA-GA has advantages in global searching and can prevent falling into a local optimum, and demonstrate that WOA-GA is competitive with other algorithms.Hence, although the time consumption is high, the proposed WOA-GA can provide good exploration and exploitation capabilities.Figure 2 shows the convergence curves from WOA-GA, WOA, GA, IWOA-PSO and ACS.WOA, IWOA-PSO and ACS have an advantage on the local space search and can converge rapidly, but they are easily to trap in a local optimum in low iterations.GA has the advantage on global space search, but it converges slowly and needs more generations to obtain the best fitness.In each iteration of WOA-GA, WOA follows several generations of GA.Thus, the global search capacity of WOA is enhanced by GA, and the local search ability of GA is improved by WOA.With the positive interaction between WOA and GA, the exploration and exploitation abilities of WOA-GA are improved.The proposed WOA-GA can converge rapidly and obtain the best fitness compared with WOA, GA, IWOA-PSO and ACS.

Application 2: WOA-GA for Time-Optimal Trajectory Planning Problem
A typical multi-point trajectory task for a 6-DOF industrial robot was employed [12].The desired positions and kinematic constraints of test task are listed in Tables 7 and 8, respectively.Since the order of magnitudes of the integral of the squared jerk is much

Application 2: WOA-GA for Time-Optimal Trajectory Planning Problem
A typical multi-point trajectory task for a 6-DOF industrial robot was employed [12].The desired positions and kinematic constraints of test task are listed in Tables 7 and 8, respectively.Since the order of magnitudes of the integral of the squared jerk is much greater than that of the total execution time, w T and w J were set as 0.9999 and 0.0001 to achieve the trade-off between two objectives.For all algorithms, a test population size of 30 and the maximum iteration of 100 and 500 were utilized to solve this problem.WOA-GA was performed 10 times as well as WOA, GA, IWOA-PSO and ACS.Table 9 shows the comparisons of the optimization results obtained using WOA-GA, WOA, GA, IWOA-PSO and ACS.WOA-GA is able to converge best in the test task, and the mean and standard deviations are also the smallest, which implies that WOA-GA has a good stability.Figure 3 shows the distributions of the best-so-far solutions from WOA-GA, WOA, GA, IWOA-PSO and ACS.Additionally, the percentage of performance rank in Table 9 indicates that WOA-GA mostly achieves a better performance than other benchmark algorithms.The travelling time of the best solution from WOA-GA is 29.9930 s, where T 1 = 0.0030 s, T 2 = 10.8240 s, T 3 = 8.9690 s, T 4 = 10.1850s and T 5 = 0.0120 s.Figures 4-7 illustrate the trajectory curves of the industrial robot and their derivatives, including jerk, acceleration and velocity.As shown in Figure 7, the jerk curves of all joints were continuous, which resulted in the smoothness of the optimized trajectories.Furthermore, the values of jerk, acceleration and velocity were zero at the start and end movement.The motion curves were all bound under the kinematic constraints all the time.Therefore, the proposed method is an effective time-jerk optimal trajectory planning tool for industrial robots.The travelling time of the best solution from WOA-GA is 29.9930 s, where T1 = 0.0030 s, T2 = 10.8240 s, T3 = 8.9690 s, T4 = 10.1850s and T5 = 0.0120 s.Figures 4-7 illustrate the trajectory curves of the industrial robot and their derivatives, including jerk, acceleration and velocity.As shown in Figure 7, the jerk curves of all joints were continuous, which resulted in the smoothness of the optimized trajectories.Furthermore, the values of jerk, acceleration and velocity were zero at the start and end movement.The motion curves were all bound under the kinematic constraints all the time.Therefore, the proposed method is an effective time-jerk optimal trajectory planning tool for industrial robots.

Conclusions
In this article, we proposed a new methodology for time-optimal and jerk-continuous trajectory planning of industrial robots.Fifth-order B-spline was adopted to interpolate the multi-point trajectory in the joint space.The proposed WOA-GA combining the advantages of WOA and GA is utilized to solve the time-jerk optimal trajectory planning problem with nonlinear constraints.Simulations were carried out to verify the effectiveness of the proposed methodology.The comparison results of the best-so-far feasible solutions showed that WOA-GA provides the best results compared with WOA, GA, IWOA-PSO and ACS.Furthermore, the proposed WOA-GA algorithm has good stability.The optimized motion curves obtained by the proposed methodology were all bound under the kinematic constraints all the time.Therefore, these results demonstrate that the proposed methodology is an effective tool for industrial robots to generate time-jerk optimal trajectories.

Conclusions
In this article, we proposed a new methodology for time-optimal and jerk-continuous trajectory planning of industrial robots.Fifth-order B-spline was adopted to interpolate the multi-point trajectory in the joint space.The proposed WOA-GA combining the advantages of WOA and GA is utilized to solve the time-jerk optimal trajectory planning problem with nonlinear constraints.Simulations were carried out to verify the effectiveness of the proposed methodology.The comparison results of the best-so-far feasible solutions showed that WOA-GA provides the best results compared with WOA, GA, IWOA-PSO and ACS.Furthermore, the proposed WOA-GA algorithm has good stability.The optimized motion curves obtained by the proposed methodology were all bound under the kinematic constraints all the time.Therefore, these results demonstrate that the proposed methodology is an effective tool for industrial robots to generate time-jerk optimal trajectories.Future work will be devoted to extending WOA-GA to optimize other trajectory models based on polynomial and NURBS.The effectiveness of the proposed technique may further be investigated by experiments.

Figure 4 .
Figure 4. Position curves of all joints by the proposed technique (Tf = 29.9930s).

Figure 5 .
Figure 5. Velocity curves of all joints by the proposed technique (Tf = 29.9930s).

Figure 4 .
Figure 4. Position curves of all joints by the proposed technique (T f = 29.9930s).

Figure 4 .
Figure 4. Position curves of all joints by the proposed technique (Tf = 29.9930s).

Figure 5 .
Figure 5. Velocity curves of all joints by the proposed technique (Tf = 29.9930s).Figure 5. Velocity curves of all joints by the proposed technique (T f = 29.9930s).

Figure 5 .
Figure 5. Velocity curves of all joints by the proposed technique (Tf = 29.9930s).Figure 5. Velocity curves of all joints by the proposed technique (T f = 29.9930s).

Figure 5 .
Figure 5. Velocity curves of all joints by the proposed technique (Tf = 29.9930s).

Figure 6 .
Figure 6.Acceleration curves of all joints by the proposed technique (Tf = 29.9930s).Figure 6. Acceleration curves of all joints by the proposed technique (T f = 29.9930s).

Figure 6 . 14 Figure 7 .
Figure 6.Acceleration curves of all joints by the proposed technique (Tf = 29.9930s).Figure 6. Acceleration curves of all joints by the proposed technique (T f = 29.9930s).Processes 2022, 10, x FOR PEER REVIEW 12 of 14 Future work will be devoted to extending WOA-GA to optimize other trajectory models based on polynomial and NURBS.The effectiveness of the proposed technique may further be investigated by experiments.Author Contributions: Conceptualization, F.W. and Z.W.; methodology, F.W. and Z.W.; writingoriginal draft preparation, F.W., Z.W. and T.B.; writing-review and editing, F.W., Z.W. and T.B.; supervision, T.B.; project administration, Z.W.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by [General Scientific Research Project of Education Department of Zhejiang Province] grant number [Y202147343] and [Research on Public Welfare Technology Application Projects of Zhejiang Province] grant number [GG19E050005].Acknowledgments: The authors acknowledge the support and inspiration of the research on public welfare technology application projects of Zhejiang province and the general scientific research project of the education department of Zhejiang province.

Figure 7 .
Figure 7. Jerk curves of all joints by the proposed technique (T f = 29.9930s).

Table 2 .
Meaning of the symbols in Equation (1).

Table 3 .
Parameter configuration of each algorithm.

Table 4 .
Details of the test functions.

Table 5 .
Comparisons of the optimization results for the unimodal functions.

Table 6 .
Comparisons of the optimization results for the multimodal functions.

Table 7 .
Desired positions of the test trajectory planning task.

Table 8 .
Kinematic constraints of the test trajectory planning task.

Table 9 .
Comparisons of the optimization results for the test trajectory planning task.