Time-Optimal Trajectory Planning for the Manipulator Based on Improved Non-Dominated Sorting Genetic Algorithm II

Abstract

To address the issues of low efficiency and lengthy running time associated with trajectory planning for 6-degree-of-freedom manipulators, this paper introduces a novel solution that generates a time-optimal path for a manipulator while adhering to its kinematic limitations. The proposed method comprises several stages. Firstly, the kinematics of the manipulator are analyzed. Secondly, the manipulator’s joint-space path points are interpolated via the quintic B-spline curve. Subsequently, the non-dominated sorting genetic algorithm II (NSGA-II) is improved by applying reinforcement learning to optimize its crossover and mutation probabilities, and the variable neighborhood search (VNS) algorithm is integrated to enhance its local search capability. Finally, the joint increments and running time of the manipulator are optimized using the improved NSGA-II, and the time-optimal trajectory is then determined by simulation on MATLAB. Furthermore, compared with other conventional optimization methods, the proposed approach has reduced the total running time by 19.26%, effectively improving the working efficiency of the manipulator.

1. Introduction

Recent years have seen the rise of manipulators as critical components in numerous industries, as advancements in industrial technology have led to their growing importance. Consequently, the demand for high-performance manipulators has also increased significantly. Research on manipulators primarily includes aspects such as structural design [1], trajectory planning, and trajectory tracking [2,3]. Among these, trajectory planning is crucial in the design process of manipulators.

Many experts and scholars have conducted extensive research on various optimization problems related to manipulators, such as time-optimal [4], energy-optimal [5], impact-optimal [6], and multi-objective optimization [7], and have proposed a variety of solutions. In recent years, research and discussion have increasingly focused on developing efficient trajectory planning methods for manipulators, with a particular emphasis on minimizing the time required to complete a task.

When planning a time-optimal trajectory for manipulators, the goal is to generate a path that satisfies the manipulator’s kinematic constraints while also minimizing the time required to complete the task. One approach that was proposed in paper [8] is to use multiple quadratic radial basis functions to generate smooth motion trajectories for the manipulator. This method uses two objective functions, one based on the total running time and the other on the square of the jerk, which are minimized in order to achieve a manipulator’s optimal trajectory in time. Experimental results show that this method can generate highly smooth trajectories within a relatively short running time.

Paper [9] proposed a smooth and efficient trajectory planning method by combining an adaptive elitist genetic algorithm with singularity avoidance mechanism and quintic polynomials. The algorithm includes an elite population and an adaptive adjustment mechanism, which improves the recognition ability of the optimal solution. Furthermore, the singularity avoidance mechanism is introduced to prevent singular points in the generated trajectory. The combination of the adaptive elitist genetic algorithm with quintic polynomials generates smooth and efficient motion trajectories.

In paper [10], an enhanced teaching-based optimization algorithm is proposed, which incorporates variable neighborhood search to keep the algorithm out of local optima. The algorithm is further combined with the quintic B-spline function to generate time-optimal motion trajectories.

The particle swarm optimization algorithm was first improved by introducing population exchange operations and a tail elimination principle to enhance the algorithm’s search ability [11]. Then, fuzzy reward–penalty theory was introduced to plan manipulator trajectory accurately while effectively reducing computational complexity for the manipulator, which leads to improved optimization efficiency by minimizing redundant decoupling operations.

An improved adaptive cuckoo search algorithm was presented in paper [12] based on the continuous quintic B-spline curve algorithm for motion trajectory planning of manipulators, which can effectively generate time-optimal joint space smooth trajectories.

Although the papers mentioned above have made positive contributions to time-optimal trajectory planning for manipulators, they did not address the issue of the motion path of manipulators. Specifically, they only achieved time-optimal trajectory planning on a fixed path without considering whether that path was the optimal solution.

In contrast to the previous approaches, paper [13] introduced a method that not only plans the time-optimal trajectory for a manipulator but also searches for the optimal path, considering both aspects simultaneously. Given the start and final points, the intermediate path points and the running time are treated as variables, and the time-optimal trajectory is planned by means of an improved genetic algorithm. Although this method achieved shorter running times compared to other methods, it only utilized a single-objective genetic algorithm that treats running time as the objective function, without considering whether the path itself is the shortest.

This paper presents an improved method for time-optimal trajectory planning of the manipulator. Firstly, the quintic non-uniform rational B-spline (NURBS) curve is used to fit the manipulator trajectory passing through the path points. Then, to address the shortcomings of the traditional fast non-dominated sorting genetic algorithm II (NSGA-II), such as slow convergence speed and easy entrapment in local optimal values, the traditional NSGA-II is improved by using the reinforcement-learning algorithm to optimize the crossover and mutation probabilities and by introducing the variable neighborhood search (VNS) algorithm to enhance the capabilities of local search. The improved algorithm and other multi-objective optimization algorithms are tested on the test function to verify the reliability of the improved algorithm. Finally, the proposed approach utilizes an improved version of the NSGA-II to plan the time-optimal trajectory of the manipulator, taking into account both the path and time optimization simultaneously. The variables are designed as path points and their running time. The first objective function is designed as the total joint increment of the six joints of the manipulator along the whole trajectory, the second objective function is the total time of the current path, and a penalty function is introduced to limit the manipulator’s angular velocity, angular acceleration, and angular jerk. By utilizing the iterative improvement of NSGA-II, the proposed method generates the optimized trajectory for the manipulator, incorporating both path planning and time optimization simultaneously. This study conducts simulation verification on MATLAB using PUMA560 and demonstrates that the proposed method yields positive effects on the trajectory optimization of the manipulator compared with other trajectory planning methods.

The remaining sections of this paper are as follows. Section 2 introduces the commonly used PUMA560 manipulator and established its model. Section 3 discusses the NURBS interpolation and establishes the interpolation matrix. Section 4 presents the improved NSGA-II and compares it with other algorithms. Section 5 performs trajectory planning for the manipulator and compares it with other methods. The final section presents the conclusion of the paper.

2. The Manipulator Modelling

The PUMA560 is a six-jointed robot with three joints used for determining the end effector’s position, while the remaining three joints determine the wrist orientation. The axis of the first joint is vertical, while the axes of the second and third joints are horizontal and parallel, with a distance of a2 between them. The joint 1 and joint 2 axes intersect perpendicularly, while the axes of joints 3 and 4 intersect alternately at a distance of a3. The manipulator’s schematic diagram is presented in Figure 1, while Figure 2 shows the corresponding link coordinate system, and

Table 1. Modified PUMA560 D-H Parameters.

Joint	θi/(°)	αi−1/(°)	αi−1/mm	di/mm	Variable Range (°)
1	90	0	0	0	−160~160
2	0	−90	0	149.09	−225~45
3	−90	0	431.8	0	−45~225
4	0	−90	20.32	433.07	−110~170
5	0	90	0	0	−100~100
6	0	−90	0	0	−266~266

displays the corresponding modified D-H data [14].

When solving for the manipulator with joint angles as variables, setting up a coordinate system for each joint of the manipulator is crucial. To construct the coordinate system for each joint, the various parameters of the manipulator listed in Table 1 are needed. The range of rotation for each joint is primarily used to eliminate irrelevant solutions when the manipulator has multiple solutions during the solution process [15]. Each link can be transformed by its transformation matrix using the coordinate system shown in Figure 2 and the data provided in Table 1, as demonstrated in Equation (1):

$${}_i{}^{i-1}T=\left[\begin{array}{cccc}c{\theta }_i& -s{\theta }_i& 0& {\alpha }_{i-1}\\ s{\theta }_{i}c{\alpha }_{i-1}& c{\theta }_{i}c{\alpha }_{i-1}& -s{\alpha }_{i-1}& -d_{i}s{\alpha }_{i-1}\\ s{\theta }_{i}s{\alpha }_{i-1}& c{\theta }_{i}s{\alpha }_{i-1}& c{\alpha }_{i-1}& d_{i}c{\alpha }_{i-1}\\ 0& 0& 0& 1\end{array}\right]$$

(1)

In Equation (1), $c\theta$ represents $\cos \theta$, and $s\theta$ represents $\sin \theta$. According to Equation (1), by substituting the data of PUMA560, the forward kinematics equation of PUMA560 can be obtained:

$${}_6{}^{0}T={}_1{}^{0}T({\theta }_1){}_2{}^{1}T({\theta }_2){}_3{}^{2}T({\theta }_3){}_4{}^{3}T({\theta }_4){}_5{}^{4}T({\theta }_5){}_6{}^{5}T({\theta }_6)$$

(2)

Simplifying Equation (2) yields the forward kinematics equation for PUMA560:

$${}_6{}^{0}T=\left[\begin{array}{cccc}{a}_x& b_x& c_x& d_x\\ a_y& b_y& c_y& d_y\\ a_z& b_z& c_z& d_z\\ 0& 0& 0& 1\end{array}\right]$$

(3)

The variables in the equation are defined as Equation (4):

$$\left\{\begin{array}{l}{a}_x=c_1\left[c_{23}(c_4{c}_5{c}_6-s_4{s}_6)-s_{23}{s}_5{c}_6\right]+s_1(s_4{c}_5{c}_6+c_4{s}_6)\\ a_y=s_1\left[c_{23}(c_4{c}_5{c}_6-s_4{s}_6)-s_{23}{s}_5{c}_6\right]-c_1(s_4{c}_5{c}_6+c_4{s}_6)\\ a_z=-s_{23}(c_4{c}_5{c}_6-s_4{s}_6)-c_{23}{s}_5{c}_6\\ b_x=c_1[c_{23}(-c_4{c}_5{s}_6-s_4{c}_6)+s_{23}{s}_5{s}_6]+s_1(c_4{c}_6-s_4{c}_5{s}_6)\\ b_y=s_1[c_{23}(-c_4{c}_5{s}_6-s_4{c}_6)+s_{23}{s}_5{s}_6]-c_1(c_4{c}_6-s_4{c}_5{c}_6)\\ b_z=-s_{23}(-c_4{c}_5{s}_6-s_4{c}_6)+c_{23}{s}_5{s}_6\\ c_x=-c_1(c_{23}{c}_4{s}_5+s_{23}{c}_5)-s_1{s}_4{s}_5\\ c_y=-s_1(c_{23}{c}_4{s}_5+s_{23}{c}_5)+c_1{s}_4{s}_5\\ c_z=s_{23}{c}_4{s}_5-c_{23}{c}_5\\ d_x=c_1[a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}]-d_2{s}_1\\ d_y=s_1[a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}]+d_2{c}_1\\ d_z=-a_3{s}_{23}-a_2{s}_2-d_4{c}_{23}\end{array}\right.$$

(4)

Equation (4) describes the position and orientation of coordinate system {6} relative to the base coordinate system {0}, where $c_{23}=\cos \left({\theta }_2+{\theta }_3\right)$ and $s_{23}=\sin \left({\theta }_2+{\theta }_3\right)$.

3. The Quintic NURBS Interpolation

3.1. The Fundamental Formulas of NURBS

Non-uniform rational B-splines (NURBS) is a mathematical curve representation method that uses a small number of control points and weights to describe a curve. NURBS curves have good local support characteristics, as changing one control point only affects the adjacent control points rather than the entire curve [16]. Using NURBS functions with non-uniform distribution of knot vectors can effectively prevent abrupt changes. The NURBS curves are represented by piecewise rational polynomial vector functions as Equation (5):

$$P(u)=\frac{{\displaystyle \sum _{i=0}^n{\omega }_i{d}_i{N}_{i,k}(u)}}{{\displaystyle \sum _{i=0}^n{\omega }_i{N}_{i,k}(u)}}$$

(5)

where $d_i$ represents the control points, ${\omega }_i$ represents the weight factors of the control points, and $N_i$ is the basis function of the B-spline function [17], which can be expressed from the recursion formula as Equations (6) and (7):

$${\mathrm{N}}_{i,0}(u)=\left\{\begin{array}{l}1\begin{array}{cc}& u_i\le u\le u_{i+1}\end{array}\\ 0\begin{array}{cc}& \begin{array}{cc}& others\end{array}\end{array}\end{array}\right.$$

(6)

$$N_{i,k}(u)=\frac{u-u_i}{u_{i+k}-u_i}{N}_{i,k-1}(u)+\frac{u_{i,k+1}-u}{u_{i+k+1}-u_{i+1}}{N}_{i+1,k-1}(u)$$

(7)

To obtain the trajectory passing through $m$ path points. The NURBS curve requires the definition of control points $d_i\left(i=0,1,\cdots ,m\right)$ and the knot vector $U=[u_0,u_1,$$\cdots ,u_{n+k},u_{n+k+1}]$. Regulating $u_0=u_1=\cdots =u_k=0$, $u_n=u_{n+1}=\cdots u_{n+k+1}=1$, and the other knot vectors are determined by the time intervals $h_j$ between the path points using the chord-length parameterization method [18] where $n=m+k-1$.

$$u_i=\frac{{\displaystyle \sum _{j=0}^{i-k-1}{h}_j}}{{\displaystyle \sum _{j=0}^{m-1}{h}_j}}(i=k+1,k+2,\dots ,n)$$

(8)

3.2. Matrix Expression of Quintic NURBS Curve

Based on the derivative formula of k-th NURBS curves, the angular velocity, angular acceleration, and angular jerk of the manipulator at the moment $u$ can be obtained by Equation (9):

$$P^{(k)}(u)=\frac{A^{}(u)-{\displaystyle \sum _{i=1}^k{C}_k^i{\omega }^{}(u)P^{(k-i)}(u)}}{\omega (u)}$$

(9)

where $A(u)={\displaystyle \sum _{i=0}^n{N}_{i,k}^{}(u)d_i{\omega }_i}$, $\omega (u)={\displaystyle \sum _{i=0}^n{N}_{i,k}^{}(u){\omega }_i}$.

According to Equation (10):

$$Q_m={\displaystyle \sum _{i=0}^{n+2}{N}_{i,k}(u)d_i\hspace{1em}\hspace{1em}i=0,1,\cdots ,\begin{array}{cc}n& \end{array}m=2,3,\cdots ,n-1}$$

(10)

To solve the NURBS curve with n + 1 unknowns, four endpoint constraints need to be added to make the system of equations have a unique solution [19]. Therefore, four constraint equations are introduced as Equation (11):

$$\left\{\begin{array}{l}{d}_0=Q_0\\ d_{n+2}=Q_n\\ d_0=d_1\\ d_{n+1}=d_{n+2}\end{array}\right.$$

(11)

This constraint equation constrains the angular velocity and angular acceleration at the start and final positions of the manipulator. Substituting it into Equation (5) yields a matrix equation for solving all the control points.

$$\left[\begin{array}{lll}{M}_0& & \\ M_1& \cdots & 0\\ M_2& & \\ & \ddots & \\ & & M_{n-2}\\ 0& \cdots & M_{n-1}\\ & & M_n\end{array}\right]\left[\begin{array}{c}{d}_0\\ d_1\\ d_2\\ ⋮\\ d_{n-2}\\ d_{n-1}\\ d_n\end{array}\right]=\left[\begin{array}{c}{Q}_0\\ Q_1\\ Q_2\\ ⋮\\ Q_{n-2}\\ Q_{n-1}\\ Q_n\end{array}\right]$$

(12)

where $M_i=\left[N_{i,5(u)},N_{i+1,5(u)},N_{i+2,5(u)},N_{i+3,5(u)},N_{i+4,5(u)}\right]$, $i=2,3,\cdots ,n-2$, $u$ is the time corresponding to the $U_i$, $M_1=\left[\begin{array}{ccccccc}-1& 1& 0& \cdots & 0& 0& 0\end{array}\right]$, $M_{n-1}=[\begin{array}{cccc}0& 0& 0& \cdots \end{array}$$\begin{array}{ccc}0& -1& 1]\end{array}$, $M_0=\left[\begin{array}{ccccccc}1& 0& 0& \cdots & 0& 0& 0\end{array}\right]$, $M_n=\left[\begin{array}{ccccccc}0& 0& 0& \cdots & 0& 0& 1\end{array}\right]$, $Q_1$ $=Q_{n-1}=0$. The remaining points are the known numerical values.

The displacement equation of the manipulator can be obtained from Equation (5), and the angular velocity, angular acceleration, and angular jerk equations of each joint of the manipulator can be obtained from Equation (9).

4. The Improved NSGA-II

The NSGA-II optimization method, proposed by Strinivas and Deb in 2002, is an improved version of the NSGA. NSGA-II was developed based on the non-dominated sorting concept of Goldberg and is derived from the genetic algorithm [20]. To distribute the population into different non-dominated levels, NSGA adopts a slow-sorting method, and to maintain the diversity of the population, it utilizes the sharing function method [21]. Specifically, NSGA-II is a non-dominated sorting algorithm that utilizes a fast-sorting method and has a worst-case computational complexity that is relatively low. It has been proven in practice that NSGA-II has improved in terms of optimization effect and computation time compared to NSGA and is an excellent multi-objective optimization algorithm [22,23,24].

Some experts and scholars have also made improvements to NSGA-II. In paper [25], the K-means algorithm was used to divide each generation of parent populations into two categories, which were selected for genetic evolution through tournament selection, and then, the differential evolution (DE) algorithm was integrated to effectively improve the algorithm’s search capability. In paper [26], an adaptive NSGA-II(AONSGA-II) was proposed in which the crossover and mutation probabilities change with the iteration number and was combined with the 2-optimization algorithm, effectively improving the algorithm’s performance.

The genetic evolutionary behavior of NSGA-II mainly depends on the crossover and mutation operations. Due to the randomness of the probabilities of crossover and mutation, their crucial parameters cannot be adaptively optimized based on the current population state during the algorithm evolution process. Therefore, there are also drawbacks, such as in the context of solving large-scale multi-objective task scheduling problems, some common issues that arise include the risk of getting stuck in local optima and the potential for decreased computational efficiency, resulting in the efficiency and quality of the solution failing to meet the requirements. In previous algorithms, the probabilities of crossover and mutation are often determined through extensive experimentation or prior experience, which are not conducive to the preservation of excellent individuals in the population. During the population evolution iteration process, to maintain population diversity and preserve the individuals with higher fitness, the selection strategy should prioritize their preservation in the population while avoiding disrupting the population by excluding them from crossover and mutation. Individuals with relatively lower fitness should have a higher probability of participating in crossover and mutation so that the population can evolve towards a better direction.

Reinforcement learning focuses on the autonomous interaction process between the agent and the environment, which does not involve external interference. Reinforcement learning provides a new solution for the self-learning of evolutionary parameters.

4.1. Reinforcement-Learning Algorithm

Reinforcement learning is a trial-and-error learning method where an agent interacts with the environment and receives rewards that guide its behavior towards maximizing the total reward obtained over a sequence of actions. Reinforcement learning is different from supervised learning in connectionist learning in that the reinforcement signal in reinforcement learning is a scalar evaluation of the goodness of the actions taken provided by the environment rather than an instruction on how to produce the correct action. Since the external information provided by the environment is limited, the reinforcement-learning system (RLS) learns from its own experience by interacting with the environment and adjusting its actions based on the rewards received. In this way, the RLS acquires knowledge about the action-evaluation environment and adapts its strategies accordingly.

4.1.1. Design Q Table

SARSA and Q-learning are two main algorithms in reinforcement learning used to train the value function Q [27]. The Q-table is initialized as a matrix with a number of rows and columns equal to the number of states and actions, respectively, and all its entries are set to zero. As shown in Equation (13), the feedback on actions is updated based on the environment and state.

$$Q(state,action)=\begin{array}{cc}& \begin{array}{cccc}{a}_1& a_2& \cdots & a_n\end{array}\\ \begin{array}{l}{S}_1\\ S_2\\ ⋮\\ S_n\end{array}& \left[\begin{array}{llll}0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ 0\hfill & \ddots \hfill & \hfill & 0\hfill \\ ⋮\hfill & \hfill & \ddots \hfill & ⋮\hfill \\ 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]\end{array}$$

(13)

$$Q(state,action)=(1-\alpha )Q(state,action)+\alpha (r+\gamma \underset{a^{\prime }}{\max }Q({state}^{\prime },{action}^{\prime }))$$

(14)

$$Q(state,action)=(1-\alpha )Q(state,action)+\alpha (r+\gamma Q({state}^{\prime },{action}^{\prime }))$$

(15)

Q-learning is a model-free method that does not require previous state values and can be executed in a stochastic dynamic environment. It updates the Q-table based on the maximum expected reward using Equation (14). SARSA, on the other hand, is an improved version of Q-learning, which trains the action value function using Equation (15). Both algorithms have their own advantages, and the experimental results presented by Sutton and Barto [28] showed that SARSA has a faster convergence speed, while Q-learning achieves better performance in the end. Furthermore, paper [29] has shown that combining SARSA and Q-learning can improve the convergence performance of the single objective genetic algorithm. In this paper, we adopt the SARSA algorithm in the initial phase and switch to the Q-learning algorithm in the subsequent phase, as shown in Equation (16):

$$\left\{\begin{array}{ll}SARSA\hspace{1em}\hspace{1em}& gen\le ge{n}_{\max }/2\\ Q-learning\hspace{1em}& gen>ge{n}_{\max }/2\end{array}\right.$$

(16)

4.1.2. Design States

The objective of using the reinforcement-learning algorithm is to find the optimal solution for the NSGA-II, while considering the diversity and balance of the optimization process. The design of the environment state for the algorithm should consider the fitness of the population and other factors, such as the crowding distance. In the algorithm, the construction of the environment state mainly considers the following aspects:

$$f(x_i^{gen})=w_1{f}_1(x^{gen})+w_2{f}_2(x^{gen})+\cdots +w_n{f}_n(x^{gen})$$

(17)

$$g=\frac{{\displaystyle \sum _{i=1}^{N}f(x_i^{gen})}}{{\displaystyle \sum _{i=1}^{N}f(x_i^1)}}$$

(18)

$$P{[N^{gen}]}_{dis}={\displaystyle \sum _{i=2}^{N-1}{\displaystyle \sum _{k=1}^M(P[i+1]·f_k(x^{gen})-P[i-1]·f_k(x^{gen}))}}/(f_k^{\max }-f_k^{\min })$$

(19)

$$d=\frac{1}{P}$$

(20)

$$balance={\displaystyle \sum _{j=1}^F\min {\displaystyle \sum _{k=1}^M\left|f_k-f_{j,k}\right|}}$$

(21)

$$h=({\displaystyle \sum _{i=1}^{N}balance(i))/N}$$

(22)

$$S=c_{1}g+c_{2}d+c_{3}h$$

(23)

Among them, Equation (17) obtains the fitness value of an individual by weighting the objective function, where $gen$ represents the iteration number, $w_i$ is the weight parameter that determines the specific value based on the requirements of each objective function, $f(x^{gen})$ represents the $gen$-th generations linear weighted sum of the objective functions of individual $i$ in the $gen$-th population, $P{[N^{gen}]}_{dis}$ is the crowding distance of the individual in the $gen$-th generations population, $balance$ represents the balance degree of each individual, $N$ is the number of populations, $M$ is the number of objectives, and $F$ represents the non-dominated layers. Equation (18) normalizes the average fitness of the population by the average fitness of the first generation to reflect the convergence of the population. Equation (19) represents the total crowding distance of the population in each generation. Equation (20) converts the crowding distance into a minimum value, where a smaller value represents a larger crowding distance, reflecting the diversity of the population. Equation (21) represents the balance degree of an individual, and a smaller balance degree indicates a more uniform distribution of the individual on each objective and a better multi-objective balance. Equation (22) represents the balance degree of the population. Equation (23) obtains the total state value of the population by weighting. The average fitness and diversity of the population are beneficial to improve the quality of the entire population and obtain excellent individuals so the three weights are designed as 0.4, 0.4, and 0.2.

4.1.3. Designing Actions

Actions refer to the responses made by the $agent$ to the current environmental state. For each generation of the population, the $agent$ takes different actions based on the environmental state, that is $p_c$ and $p_m$. For $p_c$, the interval value between each action is set to 0.05, with the exploration rate ranging from 0.4 to 0.9 as an example; when 0.4 to 0.45 is selected, a random number is chosen within this range. The same method is also applied to $p_m$, which has a range of 0.01 to 0.21 and an interval value of 0.02 [29].

The action set for crossover probability is shown in

Table 2. p_c Action Value Table.

Actions	pc Ranges	Actions	pc Ranges
$a_1$	$\left[0.40,0.45\right)$	$a_6$	$\left[0.65,0.70\right)$
$a_2$	$\left[0.45,0.50\right)$	$a_6$	$\left[0.70,0.75\right)$
$a_3$	$\left[0.50,0.55\right)$	$a_7$	$\left[0.75,0.80\right)$
$a_4$	$\left[0.55,0.60\right)$	$a_9$	$\left[0.80,0.85\right)$
$a_5$	$\left[0.60,0.65\right)$	$a_{10}$	$\left[0.85,0.90\right)$

, and the action set for mutation probability is shown in

Table 3. p_m Action Value Table.

Actions	pm Ranges	Actions	pm Ranges
$a_1$	$\left[0.01,0.03\right)$	$a_6$	$\left[0.11,0.13\right)$
$a_2$	$\left[0.03,0.05\right)$	$a_6$	$\left[0.13,0.15\right)$
$a_3$	$\left[0.05,0.07\right)$	$a_7$	$\left[0.15,0.17\right)$
$a_4$	$\left[0.07,0.09\right)$	$a_9$	$\left[0.17,0.19\right)$
$a_5$	$\left[0.09,0.11\right)$	$a_{10}$	$\left[0.19,0.21\right)$

.

4.1.4. Designing Reward Functions

The $agent$ does not take actions by itself but rather continuously tries to determine which actions will result in higher rewards in order to steer the algorithm towards better directions. The crossover probability $p_c$ is adjusted by comparing the overall states of each generation using Equation (24), while the mutation probability $p_m$ is adjusted by comparing the crowding distances of each generation using Equation (25). As the goal of the objective function in this paper is to minimize, the reward is positive if the population in generation $gen$-th is better than that in generation $gen-1$-th, and negative otherwise.

$$R_c=S^{gen-1}-S^{gen}$$

(24)

$$R_m=d^{gen-1}-d^{gen}$$

(25)

4.1.5. Movement Selection Strategy

At the beginning of the algorithm, the $agent$ has no learning experience to use, and all the data in the Q-table are zeros. Therefore, the $agent$ can only explore and learn through trial and error. By constantly exploring the unknown environment, the $agent$ accumulates more and more experience and uses this knowledge and experience to guide the selection of actions. A formula $\epsilon -greedy$ that considers both exploration and exploitation is shown in Equation (26):

$$\pi (state,action)=\left\{\begin{array}{cc}\max Q(state,action)\hspace{1em}& \epsilon \ge k_{0-1}\hspace{1em}\\ action(Rand)& \epsilon <k_{0-1}\end{array}\right.$$

(26)

In the formula, $\epsilon$ is the greedy rate, and $k_{0-1}$ is a random number between 0 and 1. When the $agent$ selects an action, it can either choose the action that maximizes the Q value, which is called the greedy policy, or it can perform exploration and choose a random action.

4.2. VNS Algorithm

The variable neighborhood search (VNS) is a metaheuristic algorithm proposed by Mladenović and Hansen in 1997 [30]. The VNS algorithm expands the search space to obtain a local optimum by changing the neighborhood structure.

4.2.1. VNS Algorithm Concept

The basic idea of the VNS algorithm is to use multiple different neighborhood structures to systematically search the solution space, obtain a local optimal solution first, and then modify the neighborhood structure to explore a broader search space and discover another optimal solution. The algorithm first uses the smallest neighborhood structure for the search. If the obtained result cannot be further improved, the algorithm switches to a slightly larger neighborhood structure for the search. If the obtained result can still be improved, the algorithm returns to the smallest neighborhood structure for further search.

4.2.2. Neighborhood Design

Considering the optimization objective of this paper is the joint increments and total time of the manipulator; the time-optimal trajectory of the manipulator can be obtained by minimizing the joint increments while reducing the total time as much as possible. To achieve this goal, five neighborhoods are designed as follows:

• Randomly select a displacement variable and subtract a certain value from it.
• Randomly select two time variables and exchange them.
• Randomly select a time variable and subtract a certain value from it.
• Randomly select a displacement variable and multiply it by a random number.
• Randomly select a time variable and multiply it by a random number.

4.3. Flow Chart of the Improved NSGA-II

Combining reinforcement learning with the NSGA-II requires considering rewards, states, actions, and selection policies. The state of NSGA-II is viewed as the environment state, and each iteration updates the state from $S_{gen}$ to $S_{gen+1}$. The learning module consists of a reinforcement-learning $agent$ and a reward function $R$. The overall execution process involves four steps. Firstly, the $agent$ acquires the environmental state $S_{gen}$ at $gen$-th generation, and then selects an action $a_i$ based on the action selection policy. The $agent$ adjusts $p_c$ and $p_m$ by executing the action, and NSGA-II uses the adjusted $p_c$ and $p_m$ for genetic operations. After completing the evolution operation, the environmental state is updated from $S_{gen}$ to $S_{gen+1}$. Finally, the reward function is used to calculate the reward and update the Q-table values for the corresponding state and action. After $gen$ times iterations, the action selection for $p_c$ and $p_m$ are optimized based on past learning experiences. The flowchart is shown in Figure 3.

4.4. Comparative Verification

This paper selects four test functions of multi-objective algorithms, including ZDT1, ZDT2, ZDT3, and ZDT6, to assess the effectiveness of the improved NSGA-II [31], which are shown in

Table 4. Representations of four test functions.

Test Function	Representation	Variable Range
ZDT1	$\begin{array}{l}{f}_1(x)=x_1\\ f_2(x)=h(1-\sqrt{f_1/h})\\ h(x)=1+9{\displaystyle \sum _{i=2}^n{x}_i/(n-1)}\end{array}$	$\begin{array}{l}n=30\\ 0\le x_i\le 1,(i=1,2,\cdots ,n)\end{array}$
ZDT2	$\begin{array}{l}{f}_1(x)=x_1\\ f_2(x)=h(1-{(f_1/h)}^2)\\ h(x)=1+9{\displaystyle \sum _{i=2}^n{x}_i/(n-1)}\end{array}$	$\begin{array}{l}n=30\\ 0\le x_i\le 1,(i=1,2,\cdots ,n)\end{array}$
ZDT3	$\begin{array}{l}{f}_1(x)=x_1\\ f_2(x)=h(1-\sqrt{f_1/h}-(f_1/h)\sin (10\pi f_1))\\ h(x)=1+9{\displaystyle \sum _{i=2}^n{x}_i/(n-1)}\end{array}$	$\begin{array}{l}n=30\\ 0\le x_i\le 1,(i=1,2,\cdots ,n)\end{array}$
ZDT6	$\begin{array}{l}{f}_1(x)=1-\exp (-4x_1){\sin }^6(6\pi x_1)\\ f_2(x)=h(1-{(f_1/h)}^2)\\ h(x)=1+9{({\displaystyle \sum _{i=2}^n{x}_i/(n-1)})}^{0.25}\end{array}$	$\begin{array}{l}n=10\\ 0\le x_i\le 1,(i=1,2,\cdots ,n)\end{array}$

. These functions cover convex, concave, and discontinuous features, which can effectively evaluate the comprehensive performance of most multi-objective evolutionary algorithms. The performance of the proposed algorithm is compared with other multi-objective optimization algorithms, such as the multi-objective particle swarm optimization algorithm (MOPSO), the multi-objective evolutionary algorithm based on decomposition (MOEA/D), the strength pareto evolutionary algorithm 2(SPEA2), KNSGA-II [25], and AONSGA-II [26].

The standard NSGA-II is set population size $N=100$, iteration times $gen=300$, crossover probability $p_c=0.8$, and mutation probability $p_m=0.1$. The improved NSGA-II is set with population size $N=100$, iteration times $gen=300$, learning rate of $\alpha =0.06$, discount rate $\gamma =0.85$, and greediness rate $\epsilon =0.6$.

To avoid making the article too lengthy, only the Pareto optimal frontiers of the traditional NSGA-II and the improved NSGA-II are shown in Figure 4 for each test function. Other optimization algorithms were only compared in terms of the mean inverted generational distance (IGD) and hypervolume (HV) values, as shown in

Table 5. Comparison of IGD Mean Values.

Algorithm	Test Function
	ZDT1	ZDT2	ZDT3	ZDT6
Traditional NSGA-II	0.1533	0.2182	0.1653	0.0649
MOPSO	0.0188	0.5849	0.0438	0.1326
MOEA/D	0.0077	0.0649	0.0563	0.0033
SPEA2	0.0893	0.1069	0.0982	0.0231
KNSGA-II	0.0056	0.0796	0.0264	0.0470
AONSGA-II	0.1047	0.1606	0.1192	0.0482
Improved NSGA-II	0.0055	0.0055	0.0066	0.0044

and

Table 6. Comparison of HV Mean Values.

Algorithm	Test Function
	ZDT1	ZDT2	ZDT3	ZDT6
Traditional NSGA-II	0.5318	0.2209	0.5274	0.3330
MOPSO	0.6992	0.1534	0.5930	0.3062
MOEA/D	0.7126	0.4083	0.5597	0.3885
SPEA2	0.6027	0.3009	0.5564	0.3590
KNSGA-II	0.7175	0.5189	0.6375	0.4400
AONSGA-II	0.5908	0.2471	0.5744	0.3460
Improved NSGA-II	0.7181	0.4427	0.5986	0.3875

.

The IGD and HV mean values obtained from running each function ten times are shown in Table 5 and Table 6. The bold portion represents the algorithms with the best results in terms of IGD and HV. To ensure fairness in the comparison, the same population size and number of iterations were used for the other multi-objective algorithms, and the remaining parameters were set to commonly used values.

According to the above results, although the MOEA/D algorithm has a lower average IGD value on the ZDT6 function, it does not demonstrate greater advantages on other test functions. KNSGA-II has higher average HV values, but its IGD results are not very good. Through comparison, it can be found that the proposed improved NSGA-II has better comprehensive performance and better search capability than other optimization algorithms, providing a foundation for the subsequent trajectory planning and enabling the search for better trajectories.

5. Simulation Verification

5.1. Model Building

The proposed method’s feasibility was demonstrated by implementing it on a PUMA560 manipulator model in MATLAB, utilizing the data presented in Table 1. The model of the PUMA560 manipulator is shown in Figure 5:

Considering the working performance and limits of each joint of the manipulator, the kinematic constraints are given in the

Table 7. PUMA560 kinematic constraints.

Joint	1	2	3	4	5	6
Angular velocity V/(°/s)	100	95	100	150	130	110
Angular acceleration A (°/s2)	45	40	75	70	90	80
Angular jerk J (°/s3)	60	60	55	70	75	70

:

To fulfill the kinematic constraints of the manipulator, penalty functions were introduced on each objective function to ensure that the motion trajectory meets the kinematic constraints. The fitness function is designed as shown in Equation (27):

$$\begin{array}{l}f(1)=({\displaystyle \sum _{i=1}^6{\displaystyle \sum _{j=1}^u\left|P_{i,j}-P_{i,j+1}\right|}})\times E_a+E_b\\ f(2)=({\displaystyle \sum _{i=1}^9{t}_i)\ast E_a+E_c}\end{array}$$

(27)

In the Equation (27), $u$ represents the number of control points, and $f\left(1\right)$ represents the joint increment of the six joints of the manipulator on the current motion trajectory, $P_{i,j}$ represents the joint increment of the th joint based on the th control point of the quintic NURBS curve, $f\left(2\right)$ represents the total running time of the manipulator on the current trajectory, and $E_a$, $E_b$, and $E_c$ are penalty functions for each fitness function. When each joint exceeds the kinematic constraint, the penalty function is activated:

$$\begin{array}{ll}{E}_a=& {\displaystyle \sum _{i=1}^6(\min P_i/p_{i,\min }+}\max P_i/p_{i,\max }+\max \left|P{\prime }_i\right|/V_{i,\max }\\ & +\max \left|{P^{″}}_i\right|/a_{i,\max }+\max \left|P^{(3)}{}_i\right|/J_{i,\max })\end{array}$$

(28)

The data of $p_{\min }$, $p_{\max }$, $V_{\max }$, $a_{\max }$ and $J_{\max }$ from Table 1 and Table 7, $E_b$ and $E_c$, are added once to the adaptation thresholds of 1000° and 5 s, respectively, otherwise $E_a=1$, $E_b=E_c=0$.

5.2. Experimental Simulation

According to the parameters of PUMA560 in Table 1, the start and final points of the manipulator are designed as shown in Equation (29):

$$\begin{array}{l}{q}_{start}=[10,-10,-30,-25,20,0](°)\\ q_{final}=[-20,25,90,20,-60,120](°)\end{array}$$

(29)

The population size $N=100$, the number of iterations $gen=300$, the learning rate $\alpha =0.06$, the discount rate $\gamma =0.85$, and the greediness rate $\epsilon =0.06$. The improved algorithm was used to solve for the time-optimal trajectory of the PUMA560 manipulator, and the path points and time are shown in

Table 8. Time-optimal path and time.

Path Points	Joint						ti (s)
	1	2	3	4	5	6
1	10.000	−10.000	−30.000	−25.000	20.000	0.000	0.000
2	6.675	−6.139	−16.259	−19.999	11.178	13.371	1.325
3	3.389	−2.235	−3.317	−15.010	2.339	26.820	1.802
4	0.024	1.507	9.819	−9.995	−6.719	39.981	2.271
5	−3.298	5.546	23.202	−5.010	−15.416	53.244	2.720
6	−6.481	9.328	36.681	0.334	−23.942	66.927	3.171
7	−9.894	13.365	49.564	5.021	−33.380	80.028	3.645
8	−13.229	17.044	63.336	9.922	−42.085	93.406	4.127
9	−16.659	21.081	76.673	15.034	−51.122	106.819	4.598
10	−20.000	25.000	90.000	20.000	−60.000	120.000	5.899

. Based on the data in Table 8, the angular displacement, angular velocity, angular acceleration, and angular jerk of the manipulator are plotted in Figure 6.

As shown in Table 8 and Figure 6, the optimal time for the manipulator is 5.899 s, and the resulting angular displacement, angular velocity, angular acceleration, and angular jerk curves are smooth and continuous and all within the constraint range.

5.3. Comparative Experiment

Through comparative experiments, it has been verified that the proposed method in this paper can achieve shorter robot-arm-running time compared to other methods. The comparative experiments all use the quintic NURBS curves to fit the path points of the manipulator, with the difference being that a single-objective genetic algorithm is used to find the optimal time based on both uniform and random paths. Specifically, uniform path refers to a path where the points between the start and final points increase or decrease uniformly. Random path, on the other hand, consists of arbitrary points that satisfy the constraints given in Table 7. The experimental results are shown in

Table 9. Comparison of experimental results.

ti	Random Path	Uniform Path	Method In Ihis Paper
1	0.000	0.000	0.000
2	1.755	1.766	1.325
3	2.534	2.316	1.802
4	3.709	2.825	2.271
5	5.112	3.441	2.720
6	6.396	4.093	3.171
7	7.336	4.656	3.645
8	8.396	5.224	4.127
9	9.121	5.811	4.598
10	10.970	7.306	5.899

:

According to the data in Table 9, using the PUMA560 model generated from Figure 5, three trajectory planning images were plotted, as shown in Figure 7. The angular displacement, angular velocity, angular acceleration, and angular jerk curves for the random path are plotted in Figure 8. Similarly, the curves for the uniform path are plotted in Figure 9.

According to the data in Table 9, it can be observed that the new time-optimal trajectory planning method proposed in this paper reduces the total running time by 46.23% compared to the traditional method of random path (or fixed path) and by 19.26% compared to the uniform path method. This effectively improves the efficiency of the manipulator. The four images in Figure 6 indicate that the obtained curves satisfy the kinematic constraints in Table 7, and the curves are uniform and continuous, demonstrating the effectiveness of the generated trajectory. By comparing the three trajectory planning methods in Figure 7, it can be observed that the proposed method passes through the start and final points, does not encounter singular points, and has a similar running trajectory to the uniform path method (which has a shorter running trajectory), demonstrating the effectiveness and rationality of the proposed method. In conclusion, the new time-optimal trajectory planning method proposed in this paper has a positive significance compared to conventional methods.

To verify the search capability of the improved NSGA-II, trajectory planning for the manipulator was performed using NSGA-II, MOPSO, MOEA/D, SPEA2, KNSGA-II [25] and AONSGA-II [26]. The resulting running times for the manipulator optimization using different algorithms are shown in

Table 10. Comparison of different algorithms.

ti (s)	NSGA-II	MOPSO	MOEA/D	SPEA2	KNSGA-II	AONSGA-II
t1	0.000	0.000	0.000	0.000	0.000	0.000
t2	3.077	2.748	4.927	3.814	2.672	1.902
t3	4.810	6.543	8.132	6.206	3.257	2.864
t4	6.345	8.778	10.656	11.206	3.961	3.960
t5	8.642	11.825	12.241	16.026	4.623	5.212
t6	11.315	14.781	14.615	20.195	5.003	6.312
t7	13.978	19.323	16.819	22.908	5.707	7.876
t8	15.571	22.165	18.732	26.877	6.548	8.582
t9	16.761	26.079	21.035	29.912	7.145	10.257
t10	19.600	28.842	23.712	34.254	9.591	12.939

. Additionally, trajectory planning curve graphs were plotted for each algorithm, as shown in Figure 10, Figure 11, Figure 12 and Figure 13.

Based on the data in Table 10 and the trajectory curve figures in Figure 10, Figure 11, Figure 12 and Figure 13, all algorithms satisfy the kinematic constraints in Table 7, validating the constraining effect of the penalty function on the results. By comparing the results, it can be observed that the traditional NSGA-II achieves shorter time durations compared to other optimization algorithms, demonstrating that NSGA-II is more suitable for solving the problem presented in this paper. KNSGA-II and AONSGA-II exhibit better search capabilities than the traditional NSGA-II, resulting in shorter time durations. However, when comparing with the improved NSGA-II proposed in this paper, the total runtime is optimized to 5.899 s, which is 38.4% and 54.4% shorter than KNSGA-II (9.591 s) and AONSGA-II (12.939 s), respectively. This further confirms that the improved NSGA-II possesses superior search capabilities.

6. Conclusions

The main work of this paper is to use an improved NSGA-II to simultaneously optimize the path point and total time of the manipulator to obtain the time-optimal trajectory. However, it cannot guarantee that the obtained trajectory is the shortest path of the manipulator.

Firstly, a uniformly and continuously moving trajectory of the manipulator is obtained by using quintic NURBS curve interpolation. Then, the reinforcement-learning algorithm and VNS algorithm are introduced to enhance the search capability of the algorithm. Subsequently, the improved NSGA-II is used to optimize the joint increments and running time of the manipulator arm. Finally, the correctness of the results is verified using the PUMA560 manipulator in MATLAB and compared with other optimization methods. The results show that the proposed method can achieve shorter running time.

The main contributions of this paper are the improvement of NSGA-II and a new method of trajectory planning for manipulator on shorter running time, which seeks shorter total running time while minimizing the joint increments. Based on comparative experiments, the improved NSGA-II has better overall performance compared to other optimization algorithms. Moreover, the new time-optimal trajectory planning method proposed in this paper effectively reduces the total runtime and improves the working efficiency of the manipulator compared to traditional methods. While this paper proposes an effective trajectory planning method, there are still some limitations to be addressed. Specifically, the method does not ensure finding the global optimal path and may not be able to guarantee the optimal running time of the trajectory, both of which are highly dependent on the search capability of the optimization algorithm. Therefore, for future work, it is crucial to improve the algorithm to have better search capabilities and design more suitable fitness functions to ensure that the path obtained is the shortest and the time obtained is the optimal time.