Multi-objective trajectory planning for industrial robots using a hybrid optimization approach

Abstract

In this paper, a hybrid approach organized in four phases is proposed to solve the multi-objective trajectory planning problem for industrial robots. In the first phase, a transcription of the original problem into a standard multi-objective parametric optimization problem is achieved by adopting an adequate parametrization scheme for the continuous robot configuration variables. Then, in the second phase, a global search is performed using a population-based search metaheuristic in order to build a first approximation of the Pareto front (PF). In the third phase, a local search is applied in the neighborhood of each solution of the PF approximation using a deterministic algorithm in order to generate new solutions. Finally, in the fourth phase, results of the global and local searches are gathered and postprocessed using a multi-objective direct search method to enhance the quality of compromise solutions and to converge toward the true optimal PF. By combining different optimization techniques, we intend not only to improve the overall search mechanism of the optimization strategy but also the resulting hybrid algorithm should keep the robustness of the population-based algorithm while enjoying the theoretical properties of convergence of the deterministic component. Also, the proposed approach is modular and flexible, and it can be implemented in different ways according to the applied techniques in the different phases. In this paper, we illustrate the efficiency of the hybrid framework by considering different techniques available in various numerical optimization libraries which are combined judiciously and tested on various case studies.

1. Introduction

Robotic manipulators play a crucial role in various industrial and service applications, by helping to improve efficiency, productivity, quality control, and innovation. A key step in the exploitation process of these machines is trajectory planning [1]. It consists in determining the time history of various robot kinodynamic parameters for the control unit, ensuring the appropriate execution of the desired task. This process is generally performed by solving an appropriate optimization problem in order to improve specific robot performances.

Early papers dealing with the robot’s trajectory planning problem (TPP) were oriented toward enhancing productivity by minimizing exclusively the execution time [2–4]. They were rooted in the field of optimal control theory, which provided powerful tools to generate optimal trajectories when high-speed motion is desired. Numerous numerical methods have been proposed to solve this class of problems. They can be grouped into two families: direct and indirect methods [5, 6]. In general, indirect methods are numerical solutions of Pontryagin’s maximum principal optimality conditions. However, such techniques suffer from several drawbacks mainly, and they need to solve a difficult nonlinear multi-point shooting problem involving the use of co-state variables that are, in general, quite difficult to estimate [7–10]. Direct methods have been suggested to overcome some of these disadvantages, and they make use of different discretization schemes for state and control variables to transform the original infinite-dimensional TPP to a more tractable finite-dimensional nonlinear parametric optimization problem that can be solved using classical mono-objective optimization techniques [5, 6]. Both nonlinear deterministic programming techniques [9–12] and stochastic optimization techniques [13–17] were successfully applied to compute suboptimal trajectories for the considered robots.

The above-mentioned studies focused mainly on solving TPP for a single objective. However, in real-world applications, multiple conflicting objectives, such as minimizing cycle time, maximizing energy efficiency, reducing vibration, or minimizing joint torques, might be involved in the trajectory planning process [18]. Thus, the TPP should be handled as a multi-objective optimization problem (MOOP) where a set of compromise solutions optimizing simultaneously the considered cost functions are investigated. These solutions constitute in the objective space what we call a Pareto front (PF) [19, 20], from which the decision-maker can select the relative best robot trajectory.

Many techniques and approaches have been proposed over the past four decades to solve MOOPs, and they can be classified according to when the preferences of the decision-maker are inserted in the solution process; thus, we have a priori, interactive, a posteriori method, or no preference methods [19, 21]. Also, they can be categorized, according to the method used for handling the vector of objective functions, into scalarization methods and vectorization methods. Scalarization methods convert the MOOP into a series of parametric single-objective optimization problems which are then solved using standard mono-objective optimization techniques [22]. For example, in refs. [12, 23], authors used the weighted sum method to solve TPP involving different cost functions. Unfortunately, scalarization methods are generally regarded as local optimization approaches and suffer from a loss of information about the relative importance of the different objectives.

Vectorization methods are handling directly the original MOOP and work directly with the PF composed of all nondominated solutions. Pareto-based methods are population-based derivative-free methods, so no-gradient information is needed, and they use in general nature-inspired metaheuristics. They belong to different classes such as swarm intelligence, evolutionary computation, physics- and chemistry-based algorithms, or evolutionary strategy [24, 25]. They are generally simple to implement and regarded as global optimization approaches. This explains the recent interest of researchers in solving the TPP for robotic manipulators using this class of algorithms. For example, Saravanan et al. [26, 27] studied two evolutionary algorithms, multi-objective differential evolution algorithm (MODE) and nondominated sorting genetic algorithm (NSGA-II) to treat the TTP for serial robots. In ref. [28], a constrained multi-objective particle swarm algorithm was applied to solve the TTP in order to minimize the traveling time, the energy, and the distance of the end effector. Also, Shi et al. [29] proposed a method based on NSGA-II to optimize the robot trajectories for three objectives, namely time, energy, and jerk. The performances of MODE and NSGA-II in processing the trajectory optimization problem were also investigated in ref. [30]. In ref. [31], the TPP is solved by using NSGA-II for two objectives, traveling time and mean jerk along the whole trajectory.

More recently, an improved multi-objective ant-lion algorithm is presented in ref. [32] to handle the time–jerk–torque optimal problem for a six-axis robot. In ref. [33], the trajectory competitive multi-objective particle swarm optimization algorithm is used to search for the optimal PF for a collaborative robot. Three performance criteria were considered: total motion time, average acceleration, and average jerk. In ref. [34], authors treated also the TPP for serial manipulators passing through waypoints by minimizing execution time, energy consumption, and joint jerks using NSGA-II. In ref. [35], authors presented a method to solve the minimum time–jerk TPP for serial manipulators in the presence of obstacles using a competitive mechanism-based multi-objective particle swarm optimizer and the Z-type fuzzy membership function is proposed to select the best solutions. Cheng et al. [36] utilized NSGA-II to solve the problems of inefficient transcranial magnetic stimulation manipulator execution and uncontrollable head safety collision, by optimizing the arm trajectories, which effectively improved the operation efficiency and protected the human head from injury. In ref. [37], authors proposed a many-objective trajectory optimization method based on response surface methodology and NSGA-III to optimize the motion of a masonry robot taking into account numerous objectives and constraints.

The analysis of the previous papers indicates that most multi-objective robot TPPs were solved using population-based nature-inspired metaheuristics and proposed approximations of the PF regrouping compromise solutions. However, these algorithms do not guarantee convergence to the true PF because of the inherent stochastic- and heuristic-based search procedures. In fact, these algorithms, in general, mimick the natural process of evolution and operate by using stochastic operators such as crossover, breeding, mutation, and selection to improve the nondominated solutions [38, 39]. Thus, most of these algorithms require high computation time to converge and have a slow convergence speed to the true PF. Moreover, there is a lack of theoretical convergence proof to the true optimal PF, and it is hard to determine the appropriate stopping criterion [40, 41]. Researchers often define the iterations as the termination condition. In consequence, in order to get a satisfactory result, the number of iterations is always defined as large, which means that computation time and efficiency are unacceptable. In order to overcome these drawbacks and make the overall search procedure faster in a more guaranteed manner, it is strongly recommended to use hybrid algorithms combining population-based metaheuristics with mathematical optimization techniques having better convergence properties [20, 42–44].

In this paper, we propose to use a hybrid approach to solve the multi-objective optimization TPP. It combines population-based nature-inspired algorithms with deterministic algorithms to find a better approximation of the true optimal PF. It is organized in four successive phases. The first phase is a transcription step where the original problem is transformed into a standard MOOP by adopting an adequate parametrization scheme for the continuous robot configuration variables. In the second phase, at least one population-based metaheuristic is used to perform a global search in the solution space and reach the region near the optimal PF in a relatively small number of generations. In the third phase, a deterministic algorithm is used to enhance the quality of found solutions by performing a local search starting from each solution of the approximated PF according to an adequate scalarization method. Finally, in the fourth phase, results of the global and local searches are gathered and postprocessed using a multi-objective direct search method to enhance the quality of compromise solutions and to converge toward the true optimal PF. This judicious combination of vectorization and scalarization methods balances their exploration and exploitation capabilities to converge toward a diverse and high-quality PF. Indeed, by using various optimization techniques having different mechanisms to explore and exploit the search space, the resulting hybrid algorithm ensures not only improvement of the overall search mechanism but also the resulting high-level relay hybrid algorithm should keep the robustness of the nature-inspired algorithm while enjoying the theoretical properties of convergence of the deterministic component. It is worth noting that the proposed approach can be implemented in different ways according to the techniques being used at the successive phases. In order to illustrate the flexibility of the proposed approach, we propose to use, in the first phase, spline functions of different orders to interpolate a set of control nodes in order to generate trajectory candidates. The coordinates of these nodes become the optimization variables of the transformed problem. In the following phases, different optimization algorithms can be used, and they are available in various numerical optimization libraries such as the Matlab Global Optimization Toolbox [45], the object-oriented framework for engineering optimization EngiO [46], the open-source platform for solving optimization problems PlatEMO [38, 47], or pymoo [39]. The efficiency of the proposed hybrid approach is illustrated by solving a set of test cases. The remainder of this paper is organized as follows. In Section 2, the problem formulation of the original TPP is detailed. In Section 3, the proposed approach is developed with its four phases. Section 4 regroups different numerical examples treating the cases of two-degrees-of-freedom (dof) and six-dof robots. Finally, Section 5 concludes this work.

2. Problem formulation

Let us consider a serial manipulator with n dof required to move from an initial location defined by q ⁱⁿⁱ to a final location q ^fin . The equation of motion for the i ^th joint can be written as follows [1]:

(1) \begin{equation}\sum _{j=1}^{n}M_{ij}\left(\boldsymbol{q}\left(t\right)\right) \ddot{q}_{j}\left(t\right)+C_{i}\left(\boldsymbol{q}\left(t\right),\dot{\boldsymbol{q}}\left(t\right)\right)+G_{i}\left(\boldsymbol{q}\left(t\right)\right)=\mathbf{\Gamma }_{\boldsymbol{i}}\left(t\right)\end{equation}

where $\dot{\boldsymbol{q}}$ and $\ddot{\boldsymbol{q}}$ are, respectively, vectors of joint velocity and acceleration. M ( q ) is the inertia matrix and C ( $\boldsymbol{q}$ , $\dot{\boldsymbol{q}}$ ) is the vector of centrifugal and Coriolis forces. $G(\boldsymbol{q})$ is the vector of potential forces and $\mathbf{\Gamma }(t)$ is the vector of input actuators efforts.

The multi-objective TPP for a robotic manipulator aims to program the robot to move from an initial configuration q ⁱⁿⁱ to a final configuration q ^fin , with optimum performances, while a set of kinematic and dynamic constraints are respected. Thus, the problem consists of determining the transfer time T, the joint trajectory vector q (t) and its time derivatives as well as the corresponding actuator efforts $\mathbf{\Gamma }(t)$ such that all constraints are respected and a vector of cost functions F _obj is optimized.

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. Proposed approach

The proposed approach includes four phases (Fig. 1). The first phase is a transcription step where the original TPP is transformed into a standard MOOP by adopting an adequate parametrization scheme for the robot configuration variables. Thus, instead of looking for a set of continuous functions representing the robot state or its input controls, the optimization process searches for a set of discrete parameters. After that and based on the transformed TPP formulation, the optimization process is performed in the three following phases that constitute together a hybrid multi-objective algorithm. A hybrid multi-objective algorithm aims to capitalize on the strengths of different algorithms, leading to an improved convergence, a better exploration of the solution space, and a more effective exploitation of the inherent problem characteristics. Thus, in the second phase of the proposed approach, a population-based search technique is used to globally explore the solution space and reach the region near the targeted optimal PF in a relatively small number of generations. Then, in the third phase, a deterministic optimization technique is used to search iteratively for new nondominated points by examining the neighborhood of solutions obtained in the second phase. Finally, solutions obtained in the second and third phases are gathered in a unique population and postprocessed in the fourth phase in order to remove possible non-Pareto optimal points and to enhance the quality of the final PF.

Figure 1. A four-phase multi-objective trajectory planning approach.

In fact, different schemes of hybridization, according to low-level schemes or high-level schemes, have been proposed in the literature and showed to be more effective than using pure metaheuristics [42, 43, 48, 49]. In a low-level scheme, a given function of a metaheuristic is replaced by another algorithm. In a high-level scheme, there is no composition of the applied metaheuristics, retaining their own identities and working as a relay, where one metaheuristic takes as its inputs the output of the precedent metaheuristic, working in series and using cooperative optimization models.

In this paper, we are going to use different optimization techniques which are available in the Matlab Global Optimization Toolbox [45] and PlatEMO [38, 47]. More precisely, a global–local search strategy based on a high-level relay hybrids model with a sequential execution order will be adopted [43]. We propose to use first, at least, one population-based algorithm for obtaining a diversified approximation of PF. Then, this approximation is improved by applying a local deterministic search on each solution of the PF approximation. Finally, non-Pareto optimal points are eliminated and the quality of the final PF is enhanced by applying a direct search multi-objective algorithm. In the proposed hybrid optimization framework, the population-based algorithm plays the role of a global optimizer by searching the entire search space to find the most promising regions with a population of individuals, while a local search module locally improves the individuals of a population. The final algorithm acts like a Pareto filter.

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) [50, 51]. 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 [19, 20]:

(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 [22]. A Pareto solution is one in which an improvement in one objective requires worsening another objective [21]. 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 [52]. 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 [38, 47]:

• • NSGAII: A fast and elitist multi-objective genetic algorithm [53].

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

• • CCMO: A coevolutionary framework for constrained MOOPs [55].

• • C3M: A multistage algorithm for solving MOOPs with multi-constraints [56].

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) [57, 58] and the distance to a reference objective method (DROM) [59, 60].

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 [57]:

(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 [59]:

(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 [59]. 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 [59, 61]. 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) [62].

MOPSA is a deterministic multi-objective technique based on the PS algorithm [63], 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 [62]. The MatLab function "paretosearch" proposes an interesting implementation of MOPSA.

4. Numerical simulations

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 [9–11]. 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 [38, 47] 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 them1.

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 [25]. 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. [11] 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.1.2 Scenario 2: point-to-point trajectory planning in the presence of one obstacle

For this scenario, the robot is required to link the same extreme configurations of the first scenario, but this time, in the presence of one circular obstacle of radius 0.1m located at the point of coordinates (x = 0.45m, y = 0.25m).

This scenario is treated exclusively using the combination of optimization techniques C3M-DROM-MOPSA but using different number of free control nodes. The optimizations are conducted using the same parameters as in the previous scenario. Table III regroups some characteristic values of the obtained PFs, and Fig. 7 illustrates the results of different phases in the objectives space. We can see that the number of free control points influences the quality of the constructed PF and the computing efforts. Low AD values are obtained for low number (2 and 3) of free nodes, whereas high HI values are obtained using relatively high number (4 and 5) of free nodes. This is the consequence of modifying the size of the search space: increasing the number of free nodes makes the joint trajectory model more flexible and enables us to represent a large spectrum of solutions but, in the same time, augments the search space and makes its exploration more difficult and time-consuming.

Globally, versus the previous scenario, the motion of the robot becomes slower and involves more efforts. This is due to the presence of the obstacle. Indeed, joining the same extreme configurations involves maneuvers in order to avoid collisions. As an illustration, we have presented in Fig. 8 the fastest transfer recorded in the PF obtained using four nodes. We can see that all imposed constraints are respected and the motion is executed in a smooth way. Figure 9b illustrates the corresponding robot motion in the presence of the obstacle, and it can be compared with the previous scenario 1 (Fig. 9a).

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 [64] and to compare the results with those given in Ref. [31, 35]. 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. [31] 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. [31], 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. [35], 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.

5. Conclusion

Multi-objective trajectory planning is a powerful tool that can be used to improve the performance of industrial robots. By optimizing multiple objectives simultaneously, it can help to increase productivity, reduce energy consumption, and improve the quality. In fact, there is no single best solution to a MOOP, as the different objectives may conflict with each other. Instead, the goal is to find a set of compromise solutions constituting the Pareto-optimal front.

The recent development of powerful multi-objective optimization techniques and their availability to the research community in various digital libraries encouraged us to propose in this paper a hybrid optimization framework to generate optimal reference trajectories for robotic manipulators. These trajectories are compromise solutions of the intricate original multi-objective optimal control problem. This framework includes four phases. The first one is a transcription phase where an adequate parametrization scheme is adopted for the robot configuration variables. In the second phase, a population-based search technique is used to globally explore the solution space and reach the region near the targeted optimal PF in a relatively small number of generations. Then, in the third phase, a deterministic optimization technique is used to search iteratively for new nondominated points by examining the neighborhood of solutions obtained in the second phase. Finally, solutions obtained in the second and the third phases are gathered in a unique population and postprocessed in the fourth phase in order to remove possible non-Pareto optimal points and to enhance the quality of the final PF.

Figure 10. Pareto front for time–jerk optimal trajectory planning problem.

The proposed framework is flexible and has been implemented according to a high-level relay hybrid model using different optimization techniques. It has a robust search mechanism thanks to the population-based component, and it is having the theoretical properties of convergence of the deterministic component. It has been successfully applied to solve problems involving two- and six-dof robots under various kinodynamic constraints and optimizing various cost functions. These results have demonstrated the efficiency of our proposal and opened the door for further applications. Indeed, based on this work, other studies can be envisaged in the future. For example, the influence of the discretization scheme parameters on the final PF can be investigated in deep. Also, instead of using only one metaheuristic, the global search phase can be also performed using a combination of various population-based metaheuristics, and this idea can also be tested. As well, it is interesting to study the possibility of introducing clustering techniques in order to reduce the number of solutions proposed to the decision-maker for selecting the best compromise solution. Finally, extending the use of the proposed approach for other classes of robots (e.g., mobile robots and aerial robots) taking into account their particularities is an interesting perspective.

Figure 11. Joint’s trajectory for the PUMA 560 robot corresponding to the minimum time solution selected from the PF obtained using the combination C3M-DROM-MOPSA: (a) position (b) velocity, (c) acceleration, and (d) Jerk.