2.1. Constraints

2.1.1 Boundary conditions

Without loss of generality, we suppose that the limit configurations ( q ⁱⁿⁱ , q ^fin ) are characterized by null joint velocities and accelerations. Thus, we have

(2a) \begin{equation} \boldsymbol{q}(t = 0) = \boldsymbol{q}^{\boldsymbol{ini}}; \qquad \boldsymbol{q}(t = T) = \boldsymbol{q}^{fin} \end{equation}

(2b) \begin{equation} \dot{\boldsymbol{q}}(t = 0) = 0; \qquad \dot{\boldsymbol{q}}\left(t=T\right)=0\end{equation}

(2c) \begin{equation}\ddot{\boldsymbol{q}}\left(t=0\right)=0;\qquad \ddot{\boldsymbol{q}}\left(t=T\right)=0;\end{equation}

2.1.2 Geometric constraints

First, there are constraints imposed on joint positions:

(3a) \begin{equation}q_{i}^{min}\leq q_{i}(t)\leq q_{i}^{max}\quad i = 1,\ldots,n\quad 0 \leq t \leq T\end{equation}

But, if obstacles are present in the workspace, collisions must be avoided and the following additional constraint will hold during the transfer:

(3b) \begin{equation}B(\boldsymbol{q}(t)) = False \qquad 0 \leq t \leq T\end{equation}

Here, B denotes a Boolean function that indicates whether the robot at configuration q is in collision either with an obstacle or with itself.

2.1.3 Kinematic constraints

During the robot motion and depending on the problem at hand, limitations may be imposed on the following kinematic parameters: for i = 1, …, n

(4a) \begin{equation}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\qquad \textrm{joint velocities:}\;\; | \dot{q}_{i}(t)| \leq \dot{q}_{i}^{max}\end{equation}

(4b) \begin{equation}\!\!\!\!\!\!\!\!\!\!\!\bullet\qquad\textrm{joint accelerations:}\;\; | \ddot{q}_{i}(t)| \leq \ddot{q}_{i}^{max}\end{equation}

(4c) \begin{equation}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\qquad \textrm{joint jerks:}\;\; | \dddot {q}_{i}(t)| \leq \dddot {q}_{i}^{max}\end{equation}

These constraints define bounds ( $\dot{\boldsymbol{q}}^{max}, \ddot{\boldsymbol{q}}^{max}, \dddot {\boldsymbol{q}}^{max}$ ) on the robot kinematics performances due to technological restrictions or to the nature of the assigned task. Of course, nonsymmetrical bounds can be treated also.

2.1.4 Dynamic constraints

In regard to the robot’s actuators’ capacities, limits can be also imposed on the amplitude of the input joint torques, such as:

(5) \begin{equation}| {{\Gamma}} _{i}(t)| \leq {{\Gamma}} _{\mathrm{i}}^{\max } \qquad i = 1, \ldots,n\end{equation}

2.2. Objective functions

The vector F _obj of objective functions regroups the desired m cost functions to be optimized simultaneously during the task achievement. Each component of F _obj represents a significant physical quantity related to the robot’s behavior or the productivity of the robotized task. From the perspective of practical applications of robotic manipulators, the i ^th (i = 1,…,m) component of the vector F _obj might be represented by one of the following expressions:

(6a) \begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\quad \textrm{minimum time}:\quad J_{T}=\int _{0}^{T}dt=T \end{equation}

(6b) \begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\quad\textrm{minimum Jerk}:\quad J_{J}=\int _{0}^{T}\sum _{i=1}^{n}\left(\frac{\dddot {q}_{i}\left(t\right)}{\dddot {q}_{i}^{max}}\right)^{2}dt \end{equation}

(6c) \begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\quad\textrm{minimum efforts}:\quad J_{E}=\int _{0}^{T}\sum _{i=1}^{n}\left(\frac{{{\Gamma}} _{i}\left(t\right)}{{{\Gamma}} _{i}^{max}}\right)^{2}dt \end{equation}

(6d) \begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bullet\quad\textrm{minimum power}:\quad J_{P}=\int _{0}^{T}\sum _{i=1}^{n}\left(\frac{{\dot{\mathrm{q}}_{i}}{{\Gamma}} _{i}\left(t\right)}{\dot{\mathrm{q}}_{i}^{max}{{\Gamma}} _{i}^{max}}\right)^{2}dt \end{equation}

The TPP defined by relations (2–6) is a complex multiple objective optimal control problem. It needs to be first transformed into a standard MOOP by adopting an appropriate approximation for the joint trajectory vector q (t). The resulting MOOP can be then treated using a hybrid multi-objective optimization algorithm as explained in the following sections.

3.1. Problem transcription

The first phase of the proposed approach aims to convert the TPP defined by relations (2–6) which is a complex multiple objective optimal control problem, into a standard MOOP. This will be done by adopting an appropriate joint motion approximation. Thus, any trajectory candidate q (t) will be represented using n distinct sets of nodes: one set of N _p nodes for each joint variable (Fig. 2). The first and last nodes are fixed according to boundary conditions (2a), but the free interior nodes can be randomly positioned to produce a trial joint trajectory q _i (t) using any appropriate fitting model that accounts for conditions (2b) and (2c). Indeed, each interior node can be moved horizontally along the timescale and vertically within the admissible range defined by relation (3a).

Figure 2. Approximation of a joint variable q _i (t) using N _p nodes.

The fact that q _i (t) must go exactly through the intermediate nodes means that we have at hand an interpolation problem. One straightforward way to generate smooth curves for a given interpolation point is to use splines. There are two commonly used methods to represent a spline function, the piecewise polynomial form (PP-form) and the basis B-spline form (B-form) [Reference De Boor50, Reference De Boor51]. For example, a k ^th order spline in PP-form, for the strict increasing break sequence t ₀ = 0, t ₁, …, t _{Np + 1} = T, is by definition of any function representing $q_{i}(t)$ , that on each of the intervals [t _i , t _{i + 1}] agrees with some polynomial of degree less than k. The function $q_{i}(t)$ has in consequence continuous derivatives up to order (k-2). The main advantage of using spline functions is that it can provide high-degree accuracy while still being computationally efficient. Also, programming a robot task reduces to finding the best positions of the via nodes of n spline functions in order to build a trajectory candidate q (t), 0 ≤t≤ T, respecting all imposed kinodynamic constraints and minimizing a set of objective functions.

In what follows, we are going to adopt the spline model for representing the robot’s joint trajectories q (t). Thus, the original infinite-dimensional optimal TPP can be cast as a classical MOOP as follows [Reference Marler and Arora19, Reference Deb20]:

(7) \begin{equation}\min _{X} \boldsymbol{F}\left(\boldsymbol{X}\right)\end{equation}

subject to : h ( X ) = 0, g ( X ) ≤ 0, x _min ≤ $\boldsymbol{X}$ ≤ x _max

where

• • $\boldsymbol{X}$ = {x ₁, x ₂, …, x _p } ∈ ℝ ^p is the vector of decision variables inherent to the adopted discretization scheme, and it brings together the coordinates of the spline nodes. Note that the abscissa of the last node is the value of the transfer time and if the nodes are uniformly distributed along the timescale, the size of the problem is reduced considerably. The vectors x _min and x _max define the limits of the search space, and they are defined according to (3a) and a superior born on transfer time T is defined in order to limit the search space.

• • $\boldsymbol{F}(\boldsymbol{X})$ is the vector of m objective functions [f ₁( $\boldsymbol{X}$ ), f ₂( $\boldsymbol{X}$ ), …, f _m ( $\boldsymbol{X}$ )] representing a transcription of the original m cost functions defined by relations (6a-d). The performance vector F ( X ) maps the parameter space of $\boldsymbol{X}$ into the objective functions space. Also, in the previous formulation, pure minimization problems are considered. However, without the loss of generality, an objective which is supposed to be maximized can be multiplied by − 1 and be minimized.

• • The feasible region of any solution candidate X is delimited by a vector of equality constraints h ( X ) and a vector of inequalities constraints g ( X ) corresponding to the transcription of the original problem path constraints defined by relations (3–5).

3.2. Global search phase

When considering multiple conflicting objective functions, the concept of Pareto dominance relationship must be used to look for Pareto-optimal solutions [Reference Logist, Houska, Diehl and Van Impe22]. A Pareto solution is one in which an improvement in one objective requires worsening another objective [Reference Pereira, Oliver, Francisco, Cunha and Gomes21]. This means that a solution candidate X dominates another solution candidate Y for a vector-valued objective function F when f _i ( X ) ≤ f _i ( Y ), for i = 1, …, m and, f _j ( X ) < f _j ( Y ) for some j. The set of nondominated solutions constitutes PF, and any of these solutions is optimal and provides decision-makers with a range of options to choose from based on their preferences and priorities.

In the past 20 years, a large number of metaheuristic programming methods have been proposed to solve problems with multiple conflicting objectives and to construct efficient good approximations of PF. Among these methods, two main classes seem to be more wildly used, namely evolutionary and swarm-based techniques [Reference Sharma and Kumar52]. They are simple and robust population-based search metaheuristics, attempting to find multiple Pareto-optimal solutions in a single simulation by emphasizing multiple nondominated and isolated solutions belonging to the PF. A large number of these algorithms are regrouped in various numerical optimization libraries. These algorithms have been successfully applied to a wide range of real-world problems in various fields, such as engineering, finance, and scheduling, where multiple objectives need to be considered simultaneously. Consequently, the second phase can be performed using a large diversity of metaheuristics. Hereafter, we are going to test the following techniques from PlatEMO [Reference Tian, Cheng, Zhang and Jin38, Reference Tian, Zhu, Zhang and Jin47]:

• • NSGAII: A fast and elitist multi-objective genetic algorithm [Reference Deb, Pratap, Agarwal and Meyarivan53].

• • ANSGAIII: An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach [Reference Jain and Deb54].

• • CCMO: A coevolutionary framework for constrained MOOPs [Reference Tian, Zhang, Xiao, Zhang and Jin55].

• • C3M: A multistage algorithm for solving MOOPs with multi-constraints [Reference Sun, Zou, Liu, Yang and Zheng56].

Unfortunately, the search operators used in these metaheuristic techniques are generic and do not guarantee that Pareto optimal solutions will be found in a finite number of solution evaluations for an arbitrary problem. Thus, it is suggested to postprocessing the result of this phase by using another deterministic search strategy having better convergence properties.

3.3. Local search phase

In this third phase, a deterministic algorithm is used to search for new nondominated solutions by examining the neighborhood of solutions generated in the second phase. However, a scalarization method is necessary for converting the MOOP into a single-objective optimization problem and subsequently solved it using an adequate local deterministic search algorithm. In this paper, we are going to use two methods: the goal attainment method (GAM) [Reference Rao57, Reference Messac58] and the distance to a reference objective method (DROM) [Reference Collette and Siarry59, Reference Eichfelder60].

GAM is a variant of goal programming method and solves the MOOP defined in rel. (7) by converting it into the following nonlinear programming problem [Reference Rao57]:

(8) \begin{equation}\min _{\lambda,\boldsymbol{X}} \lambda\end{equation}

subject to:

$f_{i}(\boldsymbol{X})-w_{i}\lambda \leq f_{i}^{*}$ , i = 1,…, m

h ( X ) = 0, g ( X ) ≤ 0, x _min ≤ $\boldsymbol{X}$ ≤ x _max

here, $f_{i}^{*}$ , i = 1,…, m, represent a set of designated targets which are associated with the m objective functions [f ₁( $X$ ), f ₂( $X$ ), …, f _m ( $X$ )] corresponding to the original cost functions (6a-d). The coefficients $w_{i}\geq 0$ , i = 1,…, m, are weights specified by the user. The minimization of the attainment factor λ leads to finding a nondominated solution which under- or over-attains the specified goals to a degree represented by the quantities $w_{i}\lambda$ . The MatLab function “fgoalattain” implements the GAM method.

DROM allows also to transform a MOOP into a mono-objective one. It is based on the minimization of a sum quantity corresponding to a distance function $L_{r}$ which is defined as follows [Reference Collette and Siarry59]:

(9) \begin{equation}L_{r}=\left(\sum _{i=1}^{m}\left| f_{i}\left(\boldsymbol{X}\right)-f_{i}^{*}\right| ^{r}\right)^{\frac{1}{r}}\end{equation}

where $1\leq r\leq \infty$ and $f_{i}^{*}$ , i = 1,…, m, represent a set of specified objective targets. In case $r=2,$ we deal with the Euclidean distance, whereas, $r=\infty$ , the method is called the min-max method or the Tchebychev method [Reference Collette and Siarry59]. Note that normalization and weight coefficients can be introduced in the distance formula (9) in order to better orient the search process. The MatLab function “fminimax” implements the DROM method.

The application of GAM or DROM leads in general to a nonlinear optimization problem that can be solved using various nonlinear mono-objective optimization techniques; hereafter, we are particularly concerned by deterministic ones in order to acquire the inherent local convergence features of these techniques. In fact, the choice for an appropriate solution algorithm is based on the inherent characteristics of the problem at hand (i.e., linear, nonlinear, differentiable, nondifferentiable, and so on). In general, the deterministic algorithms can guarantee finding at least local optimum [Reference Collette and Siarry59, Reference Martins and Ning61]. Given the same initial conditions and inputs, a deterministic algorithm, such as the sequential quadratic programming (SQP), Nelder–Mead, pattern search (PS), or DIRECT algorithms, will always give the same solution. These deterministic algorithms are well documented and are also available in various numerical optimization toolbox with multiple languages.

3.4. Postprocessing phase

The optimization process in the third phase aims to generate new solutions in the neighborhood of each solution of the PF approximation obtained during the global search phase. However, this local search phase might generate non-Pareto solutions and not well-distributed solutions on the entire PF. Indeed, we may lose the solutions diversity of the initial PF approximation because the local search does not try to preserve it. As a consequence, the new-generated nondominated neighbors should be compared with existing solutions, and only nondominated solutions are saved.

An easy way to do that, without losing the local convergence characteristics of the third phase, is to postprocess points obtained by the local search (phase 3) and those constituting the approximation of PF from the global search (phase 2) by using a multi-objective direct search method such as the multi-objective pattern search algorithm (MOPSA) [Reference Custòdio, Madeira, Vaz and Vicente62].

MOPSA is a deterministic multi-objective technique based on the PS algorithm [Reference Hooke and Jeeves63], which is a well-established local derivative-free optimization algorithm. MOPSA starts the search from a set of candidate solutions delivered by both the local search process and the global search process. A set of patterns is defined based on which the solutions will be updated. These patterns determine the search directions in the objective space. The algorithm performs the space exploration according to specific search and poll steps and updates solution archives until a convergence criterion is met. Theoretically, the algorithm converges to points near the true PF [Reference Custòdio, Madeira, Vaz and Vicente62]. The MatLab function "paretosearch" proposes an interesting implementation of MOPSA.

4.1. Multi-objective trajectory planning for a planar two-dof robot

We consider in this section the case of a two-link planar manipulator executing a point-to-point task. This robot served as a benchmark for many references dealing with the minimum time TPP [Reference Fotouhi and Szyszkowski9–Reference Massaro, Lovato, Bottin and Rosati11]. According to rel. (1), the inverse dynamic of this robot is as follows:

(10) \begin{equation}M_{11}\left(\boldsymbol{q}\right)\ddot{q}_{1}+M_{12}\left(\boldsymbol{q}\right)\ddot{q}_{2}+C_{1}\left(\boldsymbol{q},\dot{\boldsymbol{q}}\right)=\mathbf{\Gamma }_{\mathbf{1}}\left(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}}\right)\end{equation}

\begin{equation*} M_{12}\left(\boldsymbol{q}\right)\ddot{q}_{2}+M_{22}\left(\boldsymbol{q}\right)\ddot{q}_{2}+C_{2}\left(\boldsymbol{q},\dot{\boldsymbol{q}}\right)=\mathbf{\Gamma }_{\mathbf{2}}\left(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}}\right) \end{equation*}

such as

• • $M_{11}(\boldsymbol{q})=I_{1}+I_{2}+m_{1}b_{1}^{2}+m_{2}(l_{1}^{2}+b_{2}^{2}+2l_{1}b_{2}\cos q_{2})+M(l_{1}^{2}+l_{2}^{2}+2l_{1}l_{2}\cos q_{2})$

• • $M_{12}(\boldsymbol{q})=I_{2}+m_{2}(b_{2}^{2}+l_{1}b_{2}\cos q_{2})+M(l_{2}^{2}+l_{1}l_{2}\cos q_{2})$

• • $M_{22}(\boldsymbol{q})=I_{2}+m_{2}b_{2}^{2}+Ml_{2}^{2}$

• • $C_{1}(\boldsymbol{q},\dot{\boldsymbol{q}})=-l_{1}\dot{q}_{2}(2\dot{q}_{1}+\dot{q}_{2})(m_{2}b_{2}\sin q_{1}+Ml_{2}\sin q_{2})$

• • $C_{2}(\boldsymbol{q},\dot{\boldsymbol{q}})=l_{1}\dot{q}_{1}^{2}(m_{2}b_{2}\sin q_{1}+Ml_{2}\sin q_{2})$

The associated data are reported in Table 1, and two scenarios are discussed for this robot. For both cases, the robot is required to move from an initial configuration $(q_{1i},q_{2i})=(0,0)$ to the final configuration $(q_{1f}, q_{2f})=(1,-0.5)rad$ with null limit velocities. Also, we are interested in solving the MOOP for a TPP involving minimizing simultaneously $F_{1}$ , the transfer time (rel. (6a)), and $F_{2}$ , the mean value of applied joint torques (rel. (6c)). The motion is to plan under constraints on joint speed and torques amplitudes. Two scenarios are treated: the first one involves transfer in free space and the second one corresponds to a transfer in the presence of one obstacle.

Table I. Data of a two-dof planar robot.

As an application of the first phase of the proposed approach, the time evolution of the joint position variable (q ₁(t), q ₂(t)) is approximated using clamped cubic spline functions interpolating a set of two uniformly distributed control points (N _p = 2). A superior-born T _max of five seconds is fixed for the transfer time.

Optimization techniques are exploited from PlatEMO and the Matlab optimization toolbox. In particular, four different population-based methods (NSGAII, ANSGAIII, CCMO, and C3M) have been selected from PlatEMO [Reference Tian, Cheng, Zhang and Jin38, Reference Tian, Zhu, Zhang and Jin47] to be used in the second phase. For these techniques, the optimization process has been performed using the following parameters: population size = 100, number of generations = 50, and the maximum number of available function evaluations which is thus fixed to 5000. For the achievement of the local search of phase 3, both “fgoalattain” and “fminimax” MatLab functions are used to perform, respectively, GAM- and DROM-based optimizations. Additionally, the final phase is executed using the “paretosearch” MatLab function. All these MatLab functions are executed using their default parameters. It is important to recall that the various optimization techniques applied in phases 2 and 3 are presented in order to illustrate the flexibility of the proposed framework and not for comparison between themFootnote 1.

Mainly two metrics are used to evaluate the quality of obtained PF at the end of phases 2 and 4, namely hypervolume indicator “HI” and average distance “AD.” HI measures the size of the dominated space corresponding to the total volume of polytope defined by Pareto points in the objective function space. The closer points move to the PF, and the more they distribute along PF, the more space gets dominated; thus, HI increases [Reference Emmerich and Deutz25]. AD is a measure of the closeness of an individual of the PF to its nearest neighbors, and it measures distance among individuals of the same rank in objective function space, so low values of AD are better.

4.1.1 Scenario 1: Free point-to-point trajectory planning task

As a first scenario, the robot is moving in a space free of obstacles. Figures (3), (4a-d), and (5) illustrate the solutions in the objective function space obtained at the end of phases 2, 3, and 4. Table II regroups the main performances recorded at the end of phases 2 and 4. The analyze of these results shows that:

Table II. Main performances recorded at the end of phases 2 and 4.

Figure 3. PF obtained at the end of phase 2 using four different population-based metaheuristics.

• • The four techniques applied during the second phase delivered different approximations of the PF with different metrics as indicated in Table II and as it is noticed in Fig. 3.

• • Figure 4 illustrate the fact that the local search of phase (3) based on GAM or DROM was able to enhance the majority of solutions obtained in phase (2). However, sometimes the optimization process delivers dominated solutions that should be filtered later. Also, we noticed that each point of the approximated PF of the second phase involved about two seconds to be treated using GAM or DROM. The duration of phase 2 depends on the size of PF obtained at the end of phase 1.

• • Additionally, as illustrated in Fig. 5, the fourth phase was able to eliminate dominated solutions and enhance the diversity of the final constructed PF. Also, although the applied optimization techniques during the second phase delivered different approximations of the PF, these differences disappeared in the final PF approximations obtained at the end of the fourth phase; they all converge to the same area in the objective function space. The difference resides in the covered range of the PF as indicated in Table II by the min/max values of the objective functions F₁ and F₂.

• • Figure 6 illustrates one selected solution from PF obtained using C3M-PS for which the transfer time is F₁ = 1.1092s, and the normalized mean square value of applied torques is about F₂ = 1.1371s. As an indication, the minimum time solution under dynamic constraint obtained in ref. [Reference Massaro, Lovato, Bottin and Rosati11] using a nonlinear optimal control approach is bang-bang in terms of actuator torques and the transfer is achieved in a duration F₁ = 1.002 S (decrease of about 10%), but the corresponding normalized mean square value of applied torques is about F₂ = 2.02s (increase of about 90%).

Figure 4. Results of local search versus approximations of PF obtained by the global search.

Figure 5. Superposition of Pareto fronts obtained using various combinations of optimization techniques at the end of phase 4.

Figure 6. Joint’s trajectory corresponding to the minimum time solution selected from the PF obtained using the combination C3M-DROM-MOPSA.

Table III. Main performances recorded at the end of phase 4 in the presence of one obstacle.

Figure 7. Results of phases 2, 3, and 4 using various number of free control points: (a) 2 points, (b) 3 points, (c) 4 points, and (d) 5 points.

4.2. Multi-objective trajectory planning for a six-dof robot

In this case, the goal is to plan optimal time–jerk trajectories for the PUMA 560 robot [Reference Corke64] and to compare the results with those given in Ref. [Reference Huang, Hu, Wu and Zeng31, Reference Zhang and Shi35]. The manipulator is required to move through a set of via points which are reported in Table IV while respecting limited kinematic performances reported in Table V.

Figure 8. Time evolution of joint trajectories corresponding to minimum time solution obtained using four nodes.

Table IV. Knot angles of the via points defining the transfer task.

Table V. Kinematic constraints.

Figure 9. Stroboscopic robot motion corresponding to the minimum time solution: (a) scenario 1 and (b) scenario 2.

Huang et al. [Reference Huang, Hu, Wu and Zeng31] adopted the fifth-order B-spline interpolation method to construct the trajectory in the joint space with null limit conditions on velocity, acceleration, and jerk. The vector of optimization variables includes only time intervals ( $h_{i}=t_{i}-t_{i-1}$ ) defining the horizontal positions of the via points along the timescale. Two objective functions are considered: the transfer time (rel. (6a)) and the global mean jerk which is a modified expression of rel. (6b). Thus, the vector $\boldsymbol{F}_{\boldsymbol{obj}}$ to be minimized is as follows:

(11) \begin{equation}\boldsymbol{F}_{\boldsymbol{obj}}=\left[\begin{array}{c} T=\sum _{i=1}^{m-1}h_{i}\\[5pt] JM=\sum _{j=1}^{n}\sqrt{\frac{1}{T}\int _{0}^{T}\dddot {q}_{j}^{2}\left(t\right)dt} \end{array}\right],(n=6, m=4)\end{equation}

This problem is treated according to the proposed approach as follows. In the first phase, a fifth-order spline function interpolation technique is used to construct joint trajectories taking into account limit conditions and via points angles of Table IV. The transfer duration T is limited to 25 s. In the second phase, the C3M technique is applied using the same parameters of the previous example. In the third phase, the local search is applied according to the DROM approach using the “fminimax” MatLab function. The final phase is executed using the “paretosearch” MatLab function. The obtained PF is illustrated in Fig. 10. The range of the final PF is [8.496s, 25s] × [8.168°/s³, 207°/s³]. In ref. [Reference Huang, Hu, Wu and Zeng31], the same problem is solved by using NSGA-II, and the obtained PF range is [9.058s, 13.96s]×[55.55°/s³, 188.98 °/s³]. Also in ref. [Reference Zhang and Shi35], this example has been treated using the constrained multi-objective particle swarm algorithm, and the range of obtained PF is [8.722s,15.072s] × [29.26°/s³ 153.41°s³]. Of course, the shortest execution time corresponds to the largest jerk index and vice versa. The quality of our final PF is better, and the minimum time solution found using our approach is the best one, and it is illustrated in Fig. 11. We can see that generated trajectories are smooth in terms of position, velocity, acceleration, and jerk for all joints. Also, boundary and intermediate conditions are respected.