2.1. Formulation of the optimization problem

This article researches the problem of trajectory optimization for the seven-degree-of-freedom(DOF) robot manipulator. Given path points in Cartesian space, these points are transformed into via-points in joint space by inverse kinematics, based on the via-points, trajectory planning is carried out. Kinematic constraints include constraints on velocity, acceleration, and jerk of each joint for manipulator, the optimization objective is minimizing functions of time-energy-jerk of manipulator. It should be stated that the article supposes that the obstacle avoidance problem of the manipulator during motion has been solved in the path planning process. The 7-DOF Sawyer robot manipulator is employed as the experimental platform, as shown in Fig. 1.

Figure 1. 7-DOF Sawyer robot manipulator.

The robot manipulator’s optimization objectives are as follows:

1. 1) Minimize the traveling time of manipulator;

2. 2) Minimize the energy consumption of manipulator;

3. 3) Minimize the absolute value of mean jerk of manipulator.

The mathematical definition of robot manipulator’s optimization objectives is: Minimize:

(1) \begin{equation} \begin{aligned} & S_1=\sum _{i=1}^{n}T_i=\sum _{i=1}^{n}(t_{i}-t_{i-1})\\ & S_2=\sum _{j=1}^N \sqrt{\int _0^T\left (\ddot{q}_j(t)\right )^2 d t/ T} \\ & S_3=\sum _{j=1}^N \sqrt{\int _0^T\left (\dddot{q}_j(t)\right )^2 d t/ T} \end{aligned} \end{equation}

Subject to:

(2) \begin{equation} \quad \begin{aligned} & \left |\dot{q}_j(t)\right | \leq V C_j,\quad j=1,2,. ., N \\[3pt] & \left |\ddot{q}_j(t)\right | \leq A C_j,\quad j=1,2,. ., N \\[3pt] & \left |\dddot{q}_j(t)\right | \leq J C_j,\quad j=1,2,. ., N \end{aligned} \end{equation}

The definition of the symbols is: $S_{1}$ is traveling time objective function, $S_{2}$ is energy consumption index objective function, $S_{3}$ is jerk index objective function, $\dot{q}_j$ is velocity of the j-th joint, $\ddot{q}_{j}$ is acceleration of the j-th joint, $\dddot{q}_{j}$ is jerk of the j-th joint, $T$ is total traveling time of the trajectory, $T_i$ is time interval between the i-th via-point and (i + 1)-th via-point, $t_i$ is i-th time node, n + 1 is the number of the via-points, N is the number of manipulator joints, $Vc_{j}$ is velocity constraint for the j-th joint, $Ac_{j}$ is acceleration constraint for the j-th joint, $Jc_{j}$ is jerk constraint for the j-th joint. $S_2$ reflects the average acceleration of joints, which can be used as an indicator to measure the energy consumption of the manipulator during the motion [Reference He, Zhang, Sun and Shi45–Reference Elnagar and Hussein47].

There are contradictory relationships between $S_1$ and $S_2$ , between $S_1$ and $S_3$ . When $S_1$ decreases, the values of $S_2$ and $S_3$ will increase and vice versa. Thereby, multi-objective optimization techniques need to be adopted to solve the complex trajectory optimization problem.

2.2. Constructing the trajectory based on quintic B-spline curves

In the process of constructing the trajectory, the quintic B-spline interpolation method is utilized to achieve trajectory planning, based on the comprehensive consideration of factors such as smoothness, reduction of oscillations, and computational resource consumption. The quintic B-spline has the characteristic of local support, and its jerk is smooth with continuous. In the trajectory optimization studied in this article, it is required that the end-effector of the manipulator passes through the given path points. According to inverse kinematics, joint angles of the manipulator at each path point are obtained, and they are used as via-points of the manipulator. The B-spline interpolation method is used in joint space to perform trajectory planning, so the manipulator can accurately pass through via-points in joint space, and the end-effector can also pass through path points in Cartesian space accurately. The control points of the B-spline curve will be calculated firstly, and then, the B-spline curve is generated by the obtained control points.

The mathematical expression of the B-spline curve is:

(3) \begin{equation} \quad q\bigl (u\bigr )=\sum _{i=0}^{m}p_{i}N_{i,k}\bigl (u\bigr ),u_{0}\le u\le{u_{e}} \end{equation}

$p_{i} (i=0,1,\ldots,m)$ are the control points of the B-spline, $N_{i,k}(u)$ is the ith kth-order B-spline basis function, u is the knot vector, for the basis function $N_{i,k}(u)$ , it can be calculated by the Cox-DeBoor recursive formula:

(4) \begin{equation} \quad \left \{\begin{array}{l} N_{i,0}\left (u\right )=\begin{cases}1,&u_{i}\leq u\lt u_{i+1} \\0,&otherwise\end{cases},\quad suppose:\frac{0}{0}=0 \\[15pt] N_{i,k}\left (u\right )=\dfrac{u-u_{i}}{u_{i+k}-u_{i}}N_{i,k-1}\left (u\right )+\dfrac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}}N_{i+1,k-1}(u). \end{array}\right. \end{equation}

The node $u=[\underbrace{u_0, \ldots, u_0}_{k+1}, u_{k+1}, \ldots, u_{m}, \underbrace{u_{e} \ldots, u_{e}}_{k+1}]$ is obtained by normalizing the time node t according to the cumulative chord length parameterization method. The number of control points is $n+k$ , and the number of knots is $n+2k+1$ . The expression of each knot vector is:

(5) \begin{equation} \quad \left \{\begin{array}{l} u_0=0 \\[5pt] u_i=u_{i-1}+\dfrac{\left |\Delta t_{i-k-1}\right |}{\sum _{j=0}^{n-1}\left |\Delta t_j\right |}, \quad i=k+1, \ldots, m \\[5pt] u_e=1 \end{array}\right. \end{equation}

$\Delta t_{k} = t_{k+1} - t_{k}$ ( $t=0,1,\ldots,n-1$ ) is the chord length vector, i represents the node sequence number, $m=n+k-1$ . The B-spline function in the knot interval $\left [u_{i},u_{i+1}\right ]$ can be expressed by Eq.(6).

(6) \begin{equation} \quad q\bigl (u\bigr )=\sum _{j=i-k}^{i}p_{j}N_{j,k}\bigl (u\bigr ),u\in \bigl [u_{i},u_{i+1}\bigr ] \end{equation}

The r-th derivative $q^r(u)$ at the point q(u) on the B-spline curve can be calculated by Eq.(7):

(7) \begin{equation} \begin{aligned} q^{r}\left (u\right )=\sum _{j=i-k+r}^{i}p_{j}^{r}N_{j,k-r}\left (u\right ),\quad u_{i}\leq u\leq u_{i+1},j=i-k+r,\ldots,i \end{aligned} \end{equation}

(8) \begin{equation} \begin{aligned} \quad p_j^s=\begin{cases}p_j&s=0\\[5pt]\big (k+1-s\big )\dfrac{p_j^{s-1}-p_{j-1}^{s-1}}{u_{j+k+1-s}-u_j}&s=1,2,\ldots,r\end{cases} \end{aligned} \end{equation}

Accordingly, the joint position, velocity, acceleration, and jerk of the B-spline trajectory at any time can be calculated. The number of known path points is n + 1, and since the number of control points to be calculated is n + k, k-1 additional constraints are required. To interpolate the path points with the quintic B-spline curve, four more constraints are needed: set that the initial point velocity and acceleration are $v_0$ and $a_0$ , respectively, and the final point velocity and acceleration are $v_f$ and $a_f$ , respectively, therefore:

(9) \begin{equation} \begin{aligned} \quad v_0=q'\bigl (u_k\bigr )=\sum _{j=k-k+1}^{k}p_j^1N_{j,k-1}\bigl (u_k\bigr )=p_1^1 \end{aligned} \end{equation}

(10) \begin{equation} \begin{aligned} \quad v_f=q'\big (u_{n+k}\big )=\sum _{j=n+k-1-k+1}^{n+k-1}p_j^1N_{j,k-1}\big (u_{n+k}\big )=p_{n+k-1}^1 \end{aligned} \end{equation}

(11) \begin{equation} \quad a_0=q''\big (u_{k}\big )=\sum _{j=k-k+2}^{k}p_j^2N_{j,k-2}\big (u_k)=p_2^2 \end{equation}

(12) \begin{equation} \quad a_{f}=q''\left (u_{n+k}\right )=\sum _{j=n+k-1-k+2}^{n+k-1}p_{j}^{2}N_{j,k-2}\left (u_{n+k}\right )=p_{n+k-1}^{2} \end{equation}

The control points of the B-spline curve are computed by $B_{N}P = Q$ , $P$ is the control points, $P=\begin{bmatrix}p_0,p_1,\ldots,p_{m-1},p_m\end{bmatrix}^T$ , and $Q$ is the vector composed of the shape value points and the derivatives of the shape value points, $Q=\begin{bmatrix}q_0,q_1,\ldots,q_{n-1},q_n,\nu _0,\nu _f,a_0,a_f\end{bmatrix}^T$ .

Te control points P can be calculated according to the following formula:

(13) \begin{equation} \quad P={B_{N}}^{-1}Q \end{equation}

The intervals between adjacent via-points determine the coefficient matrix $B_N$ , and the matrix elements are obtained by Eq.(4). After obtaining the control points, the B-spline curve passing through given via-points in joint space can be obtained, and the B-spline curve is used to construct the joint’s trajectory for the manipulator.

Figure 2. The distribution of the search agent in DBO algorithm.

3.1. Improved non-dominated sorting dung beetle optimizer (INSDBO) algorithm

Given via-points of each joint of the robot manipulator, taking time-energy-jerk as the optimization objectives, introducing the non-dominated sorting idea to improve the DBO algorithm into INSDBO algorithm. In order to increase the optimized performance of algorithm, this research utilizes strategies of piecewise linear chaotic map, Lévy flight, and dynamic weight coefficient, employing INSDBO to search the optimal time intervals between all adjacent via-points, obtaining the Pareto front of optimization objectives. For the decision vector $x^{*}\in S$ , if it satisfies $\not \exists x\in S, x \prec x^{*}$ in the decision space S, the decision vector $x^{*}$ is called the Pareto optimal solution, also known as the non-dominated solution [Reference Tian, Lu, Zhang, Cheng and Jin48], the Pareto front is $\left \{F\left (x\right )|x\in X^{*}\right \}$ .

The operational process of INSDBO algorithm is as follows:

1. a. Generating the initial parent population $P_0$ randomly, with the population size of N, and generating the subpopulation $Q_0$ with the population size of N through four dung beetle behaviors.

2. b. The parent population and the subpopulation form the population $R_i$ together, use the non-dominated sorting strategy to rank all the individuals in the population $R_i$ into layers, each layer being a non-dominated set $Z_i$ , and calculate the total local crowding degree of each non-dominated set and sort them.

3. c. Filling the individuals in the non-dominated sets into the new parent population $P_{i+1}$ according to the order of sorting, when filling a certain level of non-dominated set, if all the individuals in the set will cause the population size to exceed N, fill some individuals with higher crowding degree into the population $P_{i+1}$ according to the individual crowding degree, until the size of $P_{i+1}$ is N.

4. d. The population $P_{i+1}$ generates the subpopulation $Q_{i+1}$ through four dung beetle behaviors.

5. e. Before reaching the maximum number of iterations, loop through b $\sim$ d.

6. f. When reaching the maximum number of iterations, export the Pareto front and the corresponding optimization variables. Let M be the number of optimization objectives and N be the population size, the computational complexity of the INSDBO algorithm is $O(MN^{2}$ ). The flow of INSDBO algorithm is shown in Fig. 3.

Figure 3. The flow of INSDBO algorithm.

After transforming the DBO algorithm into multi-objective algorithm, this research uses the three strategies of piecewise linear chaotic map, Lévy flight, and dynamic weight coefficient to realize further improvements of algorithm.

3.1.1. Piecewise linear chaotic map(PWLCM) strategy

To improve the search efficiency of the algorithm, the initial population should be uniformly distributed to increase the diversity of the population; however, the original DBO does not take this into account and uses pseudorandom numbers to generate the initial population. PWLCM is a chaotic mapping that can generate a uniformly distributed initial population and avoid individuals falling into local optima prematurely effectively, the PWLCM is defined in Eq.(14):

(14) \begin{equation} \left .X_{n+1}=g(X_n,p)=\left \{\begin{array}{c@{\quad}c}X_n/p&\quad 0\lt X_n\lt p\\[4pt](X_n-p)/(0.5-p)&\quad p\leq X_n\lt 0.5\\[4pt]g(1-X_n,p)&\quad 0.5\leq X_n\lt 1\end{array}\right .\right \}, \end{equation}

1000 random numbers in the interval (0,1) are generated using two methods: pseudorandom numbers and PWLCM sequence; then, the interval is divided into 10 equal sub-intervals, and the number of random numbers in each sub-interval is counted. Fig. 4 shows the position and number distribution of the random numbers generated by the two methods. In Fig. 4, it can be seen that the random numbers of PWLCM sequence are more uniform in each sub-interval.

Figure 4. Comparison of two kinds of random number generators.(a) Distribution of pseudorandom sequences; (b) Histogram of pseudorandom sequence; (c) Distribution diagram of chaotic sequences of PWLCM; (d) Histogram of the chaotic sequence of PWLCM.

Then, PWLCM is used to improve the mode of generating the initial population in DBO for generating a uniformly distributed initial population.

3.1.2. Lévy flight strategy

The flight trajectory of certain organisms exhibits characteristics of Lévy distribution. Lévy flight belongs to a class of random walk strategies, which shows a staggered movement of relatively short and long distances [Reference Viswanathan, Buldyrev, Havlin, da Luz, Raposo and Stanley49, Reference Heidari and Pahlavani50]. Introducing the Lévy flight search mechanism can effectively improve the search diversity and escape from local optima [Reference Houssein, Saad, Hashim, Shaban and Hassaballah51, Reference Iacca, dos Santos, Celso, de and Veloso52]. The Lévy flight strategy, utilized in the search process of the stealing dung beetle, will benefit the algorithm by increasing its search diversity.

For the mathematical definition of Lévy flight, the calculation formula is as follows:

(15) \begin{equation} \quad levy=\frac{\lambda U}{\left |V\right |^{1/\beta }} \end{equation}

(16) \begin{equation} \quad \delta _u=\left \{\frac{\Gamma \bigl (1+\beta \bigr )\sin \bigl (\pi \beta/2\bigr )}{\Gamma \Bigl [\bigl (1+\beta \bigr )/2\Bigr ]\beta \cdot 2^{(\beta -1)/2}}\right \}^{1/\beta } \end{equation}

$\beta =1.5$ , U, and V belong to the normal distribution: $U\sim N\bigl (0,\mathcal{\delta }_{u}^{2}\bigr )$ , $V\sim N\left (0,\mathcal{\delta }_{_v}^2\right )$ , $\delta _\nu$ =1.

As shown in Fig. 5(a), the step length generated by Lévy flight has traversability and randomness, then, the Lévy flight strategy is incorporated into the search mode of the stealing dung beetle to increase the search diversity.

3.1.3. Dynamic weight coefficient strategy

To improve the iterative computation efficiency and convergence speed of algorithm, the dynamic weight coefficient $\omega$ is added to position update formula of the stealing dung beetle. The calculation formula of $\omega$ is:

(17) \begin{equation} \quad w=\arctan \left (1.5\times \left (1-t/T_{max}\right )\right ) \end{equation}

As shown in Fig. 5(b), the dynamic weight coefficient $\omega$ has a larger value in the early stage of iteration and will promote global search, it becomes smaller in the later stage of iteration and will promote local search.

Figure 5. The improvement strategies for stealing dung beetles. (a) Lévy flight. (b) Dynamic weight coefficient.

Combining Lévy flight strategy and dynamic weight coefficient strategy, the improved stealing dung beetle is:

(18) \begin{equation} \quad \begin{aligned} x_i(t+1) &= levy \cdot X^b + S \times g \times (|x_i(t) - X^*| + |x_i(t) - w \cdot X^b|) \end{aligned} \end{equation}

$X^b$ is the global best position, $X^{*}$ is the current local best position, $x_i(t)$ is the position information of the i-th stealing dung beetle in the t-th iteration, g is a random vector of size 1 $\times$ D that follows normally distributed, D is the dimension of the optimization problem, and S is a constant. According to above strategies, the dung beetle optimizer is improved as INSDBO. Since stealing dung beetles act as the important exploiter’s role in the dung beetle optimization algorithm, thus their enhancement will improve the performance of the algorithm.

3.2. Trajectory optimization capabilities of INSDBO for robot manipulator

When planning the trajectory in joint space, finding a optimal time-energy-jerk trajectory that passes through given via-points set is desirable. The joint via-points for trajectory planning are given in Table I, the kinematic constraints of Sawyer robot manipulator’s joints are set, the constraints are given in Table II, and the penalty function method is applied to impose the kinematic constraints throughout the optimization process [Reference Ni, Chen, Ju and Chen53]. To illustrate the performance of the INSDBO algorithm in trajectory optimization, the NSGA-II algorithm is employed to compare with the INSDBO algorithm.

First of all, trajectory optimization is refrained from adopting, the time, energy consumption, and jerk indicators of the manipulator obtained are 24.7 s, 0.1196, 0.0975 through the calculation of Eq.(1), respectively. Then, the NSGA-II and INSDBO algorithms are utilized to perform trajectory optimization. The parameter settings of the NSGA-II algorithm are as follows: crossover probability: C = 0.9, mutation probability: Mu = 0.1. Set the population size to 100 and the number of iterations to 80 for both the NSGA-II and INSDBO algorithms, performing trajectory optimization experiments.

Table I. Joint via-points for trajectory planning.

Table II. Kinematic constraints of Sawyer robot manipulator’s joints.

In terms of averages of optimal time, optimal energy consumption, and optimal jerk obtained by optimization, the three indicators of NSGA-II are as follows: 18.3449 s, 0.1154, 0.0938, respectively, and the three indicators of INSDBO are as follows: 17.0446 s, 0.0941, 0.0692, respectively. Compared with the unoptimized trajectory, the trajectory obtained from NSGA-II reduces the values of these indicators (time, energy consumption, jerk) by 25.73%, 3.51%, and 3.79%, respectively. Compared with the unoptimized trajectory, the trajectory obtained from INSDBO reduces the values of these indicators by 30.99%, 21.32%, and 29.03%, respectively. The performance of the unoptimized and the optimized trajectory is shown in Table III.

Table III. The average trajectory optimization results of the NSGA-II and INSDBO algorithms.

As shown in Table III, it can be observed that the INSDBO algorithm can generate a better trajectory with a smaller optimal time, optimal energy consumption, and optimal jerk than the NSGA-II algorithm. Fig. 6 illustrates a comparison of Pareto fronts obtained by the INSDBO and NSGA-II algorithms, which demonstrates that the INSDBO algorithm can achieve better performance on the three optimization objectives of time, energy consumption, and jerk for the manipulator.

Figure 6. The comparison of Pareto fronts obtained by the INSDBO and NSGA-II algorithms.

3.3. Selecting the trajectory with the best comprehensive performance from the Pareto solution set based on INSDBO

The Pareto solution set obtained from INSDBO consists of many solutions, each of which contains the values of the three objective functions: traveling time  $S_1$ , energy consumption  $S_2$ , and jerk  $S_3$ . After obtaining the optimal solution set, according to different usage requirements, it can be chosen as the solution from the solution set with the best performance in one aspect or the solution with the best comprehensive performance. Since the units of each objective function are different, and to avoid the domination of the function with the largest range, the objective functions need to be normalized:

(19) \begin{equation} \quad f_i=\frac{F_{i \max }-F_i}{F_{i \max }-F_{i \min }},\quad i=1,2,3 \end{equation}

$F_i$  is a set of the corresponding optimization index $S_i$ in the solution set; $F_{imax}$ is the maximum value of the corresponding optimization index $S_i$ in the solution set; $F_{imin}$ is the minimum value of the corresponding optimization index $S_i$ in the solution set; $f_i$ is the normalized set of the corresponding optimization index  $S_i$ . The normalized solutions are calculated as follows:

(20) \begin{equation} \quad f_e=N_1\cdot f_1+N_2\cdot f_2+N_3\cdot f_3 \end{equation}

Set the weight coefficients  $N_1=1/2$ , $N_2=1/4, N_3=1/4, f_e$ is a one-dimensional vector that merges three indicators of performance after normalization. The objective function corresponding to the maximum value in $f_e$ reflects the optimal selection of the manipulator in terms of comprehensive performance. This research takes one solution set of 30 experiments as the research object. The Pareto front obtained by INSDBO is shown in Fig. 7, and point B is the Pareto optimal solution selected according to Eq. (20).

Point B is a tradeoff Pareto solution, which is better than point C in terms of traveling time and better than point A in terms of energy consumption and jerk index. The time intervals between adjacent via-points of the manipulator corresponding to point B are used to plan the trajectory of Sawyer robot manipulator. Table IV shows the time intervals for the manipulator to pass through all adjacent via-points.

Table IV. Time intervals between adjacent via-points of the selected solution.

Figure 7. Pareto front of optimal time-energy-jerk trajectory based on INSDBO.

The velocity and acceleration of manipulator at initial and final points are set to zero. The trajectory of seven joints with the optimal comprehensive performance of time-energy-jerk is shown in Fig. 8, and the position, velocity, acceleration, and jerk of the manipulator are shown.

The trajectory in Fig. 8 shows that the velocity and acceleration at the start and end points are zero, which conforms to the specified requirements; meanwhile, the velocity, acceleration, and jerk meet the kinematic constraints. Moreover, the joint motion is continuous and smooth; there is no sudden change in the trajectory, ensuring the smoothness of the optimized trajectory.