Multi-objective optimal trajectory planning for robot manipulator attention to end-effector path limitation

Abstract

In the process of trajectory optimization for robot manipulator, the path that is generated may deviate from the intended path because of the adjustment of trajectory parameters, if there is limitation of end-effector path in Cartesian space for specific tasks, this phenomenon is dangerous. This paper proposes a methodology that is based on the Pareto front to address this issue, and the methodology takes into account both the multi-objective optimization of robotic arm and the quality of end-effector path. Based on dung beetle optimizer, this research proposes improved non-dominated sorting dung beetle optimizer. This paper interpolates manipulator trajectory with quintic B-spline curves, achieves multi-objective trajectory optimization that simultaneously optimizes traveling time, energy consumption, and mean jerk, proposes a trajectory selection strategy that is based on Pareto solution set by introducing the concept of Fréchet distance, and the strategy enables the end-effector to approach the desired path in Cartesian space. Simulation and experimental results validate the effectiveness and practicability of the proposed methodology on the Sawyer robot manipulator.

1. Introduction

Robot manipulator is widely used in industrial and service fields, and trajectory planning is an important component of robot research. Trajectory optimization techniques can be mainly divided into the following categories: time-optimal, energy-optimal, jerk-optimal, and hybrid-optimal, etc. The basic idea of time optimization is to maximize the speed under motion constraints, which can help improve the work efficiency of robotic arm [1, 2]; energy optimization can help manipulator save energy and improve its economic advantages [3, 4]; jerk optimization can improve the smoothness of the trajectory; moreover, the large jerk has adverse effects on the control performance [5, 6].

The trajectory planning of manipulator can be divided into joint space planning and Cartesian space planning [7, 8], control system acts on the robotic arm joints, it is easier to adjust trajectory according to requirements in joint space. The quintic polynomial curve, seventh-order polynomial curve, Bezier curve, B-spline curve, NURBS curve, etc. are among the commonly used trajectory types for the manipulator [9–11], B-spline curve has benefits of local support and jerk continuity and is highly regarded in the field of trajectory planning [12–15]. In trajectory optimization, the time intervals between adjacent path points can be used as the optimization variables, according to design time intervals for the manipulator to pass through between adjacent path points, traveling time, energy consumption, and jerk optimization can be achieved. In terms of algorithm selection of trajectory optimization, Dorian et al. [16] used the particle swarm optimization, He et al. [17] used genetic algorithm, Jin et al. [18] used chaotic particle swarm algorithm, Gurjeet et al. [19] used hybrid optimization algorithm, Huang et al. [20] adopted NSGA-II, Wang et al. [21] used INSGA-II. Trajectory optimization of the manipulator mainly includes single-objective trajectory optimization and multi-objective trajectory optimization.

In single-objective optimization, the optimization objective is usually a single factor such as time, energy, jerk, etc. In the research on time optimization of manipulator, Elias et al. [22] used genetic algorithm with multiple population to achieve the time-optimal trajectory planning; Yang et al. [23] and Liu et al. [24] adopted particle swarm to achieve time optimal for robot; Zhang et al. [25] proposed an improved sparrow search algorithm to solve the time-optimal trajectory; moreover, there are some related studies [26–28]. In the research on energy consumption optimization of the manipulator, Alexandra et al. [29] proposed a method to minimize energy consumption; Pierre et al. [30] studied the energy-optimal motion planning problem of a planar manipulator with two rotational joints. In the research on jerk optimization of the manipulator, Lu et al. [31] transformed the trajectory planning problem into a constrained nonlinear multivariate optimization problem, achieving minimum jerk trajectory planning. The drawback of employing single-objective optimization is its inability to simultaneously consider multiple aspects of manipulator’s performance.

In multi-objective trajectory optimization of robotic arm, the optimization object is usually a combination function of time-jerk, time-energy, or time-energy-jerk, etc [32, 33]. In this case, there are mainly two optimization strategies. The first strategy is to adopt a weighted method to transform the multi-objective optimization problem into a single-objective problem for solving; however, this strategy requires reasonable setting of weights; otherwise, the desired outcomes may not be achieved, while researchers must rely on their own intuition and proficiency [34, 35]. The second strategy is to utilize multi-objective optimization algorithms, which can yield a set of Pareto optimal solutions encompassing the optimal solutions for each objective function. Researchers can select a solution that favors a certain optimization objective or has a comprehensive advantage in optimization performance from the solution set based on the actual needs. Li et al. [36] proposed an optimal joint trajectory planning model based on time and energy, used quintic B-spline curve and improved NSGA-II for multi-objective optimization; Mathias et al. [37] achieved minimum motion time and energy consumption of industrial machines based on improved S-curve and NSGA-II. Liu et al. [38] proposed a trajectory planning algorithm based on improved multiverse algorithm for time, energy, and jerk optimization; Cao et al. [39] proposed an improved multi-objective particle swarm optimization algorithm, with time, energy consumption, and jerk as optimization objectives. However, they neglected the end-effector path conditions while concentrating on the joint performance of the manipulator. Besides, the dynamic response of the system has a significant influence on the implementation of trajectory planning [40–42].

A novel multi-objective optimization technique is presented to obtain multi-objective optimal trajectory for robot manipulator in this article: the improved non-dominated sorting dung beetle optimizer (INSDBO) algorithm. The dung beetle optimizer(DBO) is a efficient optimization algorithm that was proposed in 2022, which balances global exploration and local exploitation, has characteristics of fast convergence speed and high solution accuracy, received high recognition [43, 44]. In this research, INSDBO and quintic B-spline are utilized to solve the optimal time-energy-jerk trajectory planning problem with complicated constraints. The method proposed in this paper is suitable for kinds of serial manipulators. The main contribution of this paper is presenting a methodology that enables the end-effector of robotic arm to approach the desired path by introducing the concept of Fréchet distance, while focusing on the multi-objective optimization of robot manipulator performance. Meanwhile, this study is the first to apply the DBO algorithm to the robotic manipulator, besides trajectory optimization, the DBO algorithm can be applied to multiple domains for robot manipulators, such as kinematic calibration.

This article is organized as follows. In Section 2, the trajectory optimization problem that needs to be solved in this research is elaborated, quintic B-spline curve is employed to construct trajectory. In Section 3, the improved non-dominated sorting dung beetle optimizer algorithm is proposed and applied to achieve multi-objective trajectory optimization of robot arm. In Section 4, the strategy for selecting a solution in the Pareto solution set that takes into account both the comprehensive performance of trajectory and the end-effector path is presented. In Section 5, the planned trajectory is utilized to conduct simulations and actual experiments, implementing the application on the Sawyer robot manipulator. Finally, this article is summarized in Section 6.

2. Problem statement

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 [45–47].

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. Multi-objective trajectory optimization based on INSDBO algorithm

The DBO algorithm is a novel global optimization algorithm with a strong exploration ability, and it can search for both global and local optima with fast convergence and high accuracy. The DBO algorithm has four categories of dung beetles according to their natural behaviors: rolling, spawning, foraging, and stealing. At each iteration, the dung beetles update their positions based on their own rules and compare their fitness values with all other dung beetles, they keep track of the best dung beetle currently and return the global optimum when they reach the maximum iterations. The four categories of dung beetles are shown in Fig. 2.

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 [48], 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 [49, 50]. Introducing the Lévy flight search mechanism can effectively improve the search diversity and escape from local optima [51, 52]. 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 [53]. 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.

4. An optimal trajectory selection strategy with simultaneous attention to end-effector path limitation

In manipulator trajectory planning, it is common to specify the waypoints that the end-effector of the manipulator needs to follow in Cartesian space and then transform them to the joint space for trajectory planning. Sometimes, it is only required for the end-effector to pass through these waypoints, but other times there is also an expectation for it to move along a certain path. Generally, when performing trajectory planning with stringent constraints on the end-effector’s path, it is carried out in the Cartesian space of the manipulator. However, this approach may lead to certain challenges: inverse solutions may take a lot of computational resources and time, and the manipulator’s trajectory may not be smooth. Therefore, this research proposes a strategy that can take into account the path of the end-effector in Cartesian space while performing trajectory planning in joint space so as to generate an optimal time-energy-jerk trajectory that is close to the desired path.

The robot arm’s trajectory in joint space is represented by a B-spline curve, whose shape depends on the distribution of its nodes. It is worth noting that the nodes are time intervals of the robot arm’s motion. Therefore, the time intervals influence the trajectory in joint space as well as the corresponding trajectory in Cartesian space. Therefore, it can be concluded that the time intervals between adjacent via-points of the manipulator affect not only the performance of the manipulator in terms of time-energy-jerk but also the path of the end-effector in Cartesian space. Different time intervals result in different Cartesian paths. Establish a coordinate system with the base center of the Sawyer manipulator as the origin, and a set of desired via-points in Cartesian space is given in Table V.

Table V. Via-points in Cartesian space.

Figure 8. The position, velocity, acceleration, and jerk of the seven joints in the planned trajectory.

After transforming the given Cartesian space path points to the joint space using the inverse kinematics of the Sawyer robot manipulator, obtaining five path points in joint space. This research conducts two sets of experiments (Exp.1 and Exp.2) with different time intervals (6 s, 4 s, 6 s, 4 s) and (3.5 s, 4.5 s, 6 s, 6 s) between the five path points and performs trajectory planning for the manipulator. Fig. 9 shows the joint angular velocity of the Sawyer robot manipulator in both experiments, and the manipulator has different trajectories in the two experiments due to different time intervals between adjacent via-points. Fig. 10 shows the end-effector paths of the Sawyer robot manipulator in both experiments, obtained by using the forward kinematics of the Sawyer robot manipulator. The Sawyer robot manipulator in Fig. 10 is only used to illustrate the spatial relationship between the end-effector of the manipulator and the planned path and does not represent the actual pose of the manipulator (the same applies to Fig. 13).

Figure 9. The joint angular velocity of two groups of experiments. (a) Exp.1; (b) Exp.2.

Figure 10. The path of the end-effector in the Cartesian space of two groups of experiments. (a) exp.1; (b) exp.2.

This research has verified the following concept through the above example: Different trajectories lead to different Cartesian space paths for the manipulator end-effector, even if the given via points and the total motion time are the same. The time interval between adjacent path points affects the joint space motion trajectory, which in turn affects the Cartesian space path of the manipulator end-effector. Therefore, when expecting the manipulator end-effector to pass through given via points and approach a desired path at the same time, this can be achieved by setting appropriate time intervals when planning the trajectory in joint space.

Based on the above analysis, the following strategy is proposed: Select a solution from the Pareto optimal solution set that has the most similar end-effector path to the desired path and use its corresponding time intervals for trajectory planning in joint space. The solution is chosen from a region where the three objectives of time, energy consumption, and jerk are well balanced. The region excludes 30% of the solutions at both ends of the objective function range and includes 40% of the solutions in the middle. These solutions consider the performance of the manipulator in terms of time-energy-jerk and avoid extreme tradeoffs between different aspects; the performance in all aspects is at a satisfactory level. Fig. 11 describes the candidate solutions (red points set); these solutions are obtained through trajectory optimization using the INSDBO algorithm (via points in Cartesian space are from Table V); each solution corresponds to a set trajectory of the manipulator in joint space as well as a corresponding end-effector path in Cartesian space.

Figure 11. Solutions with relatively superior overall performance.

The Fréchet distance can measure the distance between two curves in space [54, 55]. This research utilizes the Fréchet distance to select a solution from candidate solutions to obtain a trajectory that has the nearest distance from the end-effector path to the desired path. The Fréchet distance is defined as Eq.(21), herein, M and N are two continuous curves in Cartesian space, $\alpha$ and  $\beta$ are two reparametrization functions of the unit interval, and d is a metric function in space. The curve is defined as a continuous mapping $M:[a,b]\rightarrow V$ , where  $a\leq b$ and (V,d) is a metric space, $N : [a',b']\rightarrow V$ has similar properties, and $\alpha (resp.\beta)$ is an arbitrary continuous nondecreasing function from [0,1] onto [a,b] (resp. $[a^{\prime },b^{\prime }]$ ).

(21) \begin{equation} \quad \delta _F(M,N)=\inf _{\alpha,\beta }\max _{t\in [0,1]}d(M(\alpha (t)),N(\beta (t))) \end{equation}

$\delta _F(M,N)$  is employed to measure the distance deviation between the planned path M and the desired path N. When the generated path is closer to the desired path, the $\delta _F$ is smaller; if reversed, the $\delta _F$ is larger. By selecting the solution from the candidate solutions with the smallest $\delta _F$ , a trajectory can be identified whose corresponding end-effector path is closest to the desired path. The shortest  $\delta _F$ value is 0.03, and the time interval between adjacent via-points is (8.8598s, 8.1948s, 7.1527s, 8.3359s); the largest $\delta_F$ value of 0.048 corresponds to the time interval between adjacent via-points of (6.7926s, 8.6074s, 6.3705s, 6.9301s). The trajectories corresponding to these two groups of time intervals between adjacent via-points are defined as Trajectory I and Trajectory II, respectively. The corresponding joint angular velocity of these two groups of trajectories is shown in Fig. 12.

Figure 12. The joint angular velocity of trajectory I and trajectory II. (a) Trajectory I; (b) Trajectory II.

Figure 13. The Cartesian space path of the end-effector. (a) Corresponding Cartesian space path of trajectory I; (b) Corresponding Cartesian space path of trajectory II.

In Fig. 13, the blue dots are the 5 path points selected on the desired path, which are used as the via-points for the manipulator; the green curve is the desired path in Cartesian space; and the red curve is the end-effector’s motion path of the manipulator. It can be seen that in Fig. 13(a), the path of the end of the manipulator corresponding to Trajectory I is closer to the desired path, while in Fig. 13(b), the path of the end of the manipulator corresponding to Trajectory II is farther from the desired path. Trajectory I results in the most similar end-effector path to the desired path; therefore, this research chose trajectory I to conduct an experiment on a Sawyer robot manipulator. The trajectory allows the manipulator to balance traveling time, energy consumption, and jerk performance during the motion and also makes the end-effector follow the desired path closely.

5. Experiments and results

To validate the practicality of the optimized trajectory on the manipulator, this research conducts simulation and experimentation using the 7-DOF Sawyer robot manipulator as experimental platform. Numerical computation is performed on the computer with an Intel (R) Core (TM) i7-4790 CPU @ 3.60 GHz and 16.0 GB RAM. The simulation platform is robot simulator CoppeliaSim (formerly V-REP), while the experimental platform is the actual Sawyer robot (Rethink Robotics Inc.).

To visualize the movement of the Sawyer robot manipulator end-effector, the planned trajectory is tested using CoppeliaSim. The initial configuration is set, the trajectory I that balances traveling time, energy consumption, jerk, and end-effector path simultaneously in Section 4 is used as the desired trajectory, and it is sent to the robot. The motion process of the manipulator in CoppeliaSim and the motion path of the end-effector are shown in Fig. 14.

Figure 14. The motion path of Sawyer robot manipulator end-effector.

The results in Fig. 14 show that the end-effector of the manipulator moves along the path close to the desired path, demonstrating that the optimized trajectory can enable the end-effector to complete the desired motion, and the manipulator motion posture corresponding to the trajectory meets the expectation.

During the experiment, the trajectory that takes into account the Cartesian space path, traveling time, energy consumption, and jerk is generated on the computer. Subsequently, the planned trajectory is sent to the robot controller within the Robot Operating System (ROS) framework. Communication between the computer and the robot controller is established using the standard TCP/IP protocol. The controller employs inverse dynamics feedforward position mode control to track the desired trajectory. The 7-DOF Sawyer robot manipulator is used as the experimental platform.

Based on the desired trajectory, the inverse dynamics feedforward predicts the system’s response and generates corresponding control torque, and the closed-loop control is responsible for the accuracy, stability, and responsiveness of position control. The inverse dynamics feedforward position mode control effectively handles uncertainties and time delays in the system, promptly and precisely tracking the desired trajectory, thus ensuring stable and reliable system responses.

Taking the real application scenario of a manipulator to deliver a tool (electrical tape) as an example, give the initial point and target point of the manipulator end-effector in Cartesian space, where the initial point is in a higher position that the manipulator needs to reach to pick up the tool, and the target point is in a lower position that the manipulator needs to reach to place the tool. The trajectory I is chosen as the motion trajectory for the manipulator, and then, based on the ROS action, the trajectory is sent to the Sawyer robot manipulator for execution, and the change of motion state for the Sawyer robot manipulator is shown in Fig. 15. The measured joint angles of the robot manipulator in the actual motion are shown in Fig. 16.

Figure 15. Robot manipulator tool delivery experiment based on planned trajectory.

Figure 16. Measured joint space trajectory of the Sawyer robot manipulator. (a) Joint1; (b) Joint2; (c) Joint3; (d) Joint4; (e) Joint5; (f) Joint6; (g) Joint7.

As shown in the experiment, the joint trajectory tracking performance is excellent, the experimental profile is basically consistent with the planned profile, and the joint movement is overall stable, with no drastic changes. The experimental trajectory closely aligns with the planned trajectory, exhibiting basically consistent amplitude and trend throughout. The maximum positional errors for the seven joints are around 0.0012 rad, 0.001 rad, 0.0015 rad, 0.0025 rad, 0.0056 rad, 0.0049 rad, and 0.0021 rad, respectively. At the final stage of the trajectory curve, the error of all joints is within 0.0025 rad. The satisfactory optimized trajectory is achieved on the actual robot manipulator; the result validates the applicability of the proposed trajectory optimization methodology on the actual manipulator and simultaneously indicates that the planned trajectory is friendly to the controller.

The experimental results show that the manipulator can perform the task of delivering the tool effectively based on the optimized trajectory by INSDBO, demonstrating that the planned trajectory can be well applied to the actual robot manipulator.

6. Conclusions

This research proposes a novel methodology of trajectory optimization with an improved multi-objective dung beetle optimization algorithm, applies it with quintic B-spline curves to trajectory planning for a robot manipulator, and obtains a trajectory with optimal traveling time, optimal energy consumption, optimal jerk, and a satisfying end-effector path. The main contribution of this methodology is acquiring a superior performance trajectory in joint space while enabling the path of the manipulator end-effector to approach the desired path in Cartesian space. The obtained trajectory is verified through experiments.

1. (1) This study pioneers the application of the DBO algorithm to the trajectory optimization of manipulators, which extends its potential to the domain of robotic systems and offers enlightening implications for future researchers.

2. (2) The performance of the INSDBO algorithm is favorable for traveling time, energy consumption, and jerk in trajectory optimization for manipulators. In terms of the average of optimal solutions in 30 times of experiments, 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.

3. (3) The proposed trajectory selection strategy integrates Fréchet distance and Pareto solution set, which can guarantee manipulator achieves optimal tradeoff among traveling time, energy consumption, and jerk objectives while enabling the end-effector of the manipulator to follow the desired path closely during motion.

4. (4) The actual Sawyer robot manipulator is employed to experimental verification, and the optimized trajectory is applied to tool delivery experiment, showing the applicability of the proposed trajectory planning methodology.

The main limitation of the study is that it focuses primarily on the trajectory planning aspect to achieve the actual path approximating to the desired path, while control accuracy also plays an important role in achieving this objective. In addition, because the methodology plans the trajectory offline and uses it online, the value can be maximized when the planned trajectory needs to be reused. In future work, trajectory tracking with the LQR or SMC controller will be performed to enhance the tracking accuracy. Additionally, the subsequent research focus will delve into the field of online trajectory planning.