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 [Reference Craig1]. 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 [Reference Kahn2–Reference Bobrow, Dubowsky and Gibson4]. 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 [Reference Betts5, Reference Bruce6]. 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 [Reference von Stryk and Bulirsch7–Reference Ghasemi, Kashiri and Dardel10]. 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 [Reference Betts5, Reference Bruce6]. Both nonlinear deterministic programming techniques [Reference Fotouhi and Szyszkowski9–Reference Chettibi, Lehtihet, Haddad and Hanchi12] and stochastic optimization techniques [Reference Chettibi and Lehtihet13–Reference Latombe17] 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 [Reference Gasparetto, Boscariol, Lanzutti and Vidoni18]. 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) [Reference Marler and Arora19, Reference Deb20], 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 [Reference Marler and Arora19, Reference Pereira, Oliver, Francisco, Cunha and Gomes21]. 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 [Reference Logist, Houska, Diehl and Van Impe22]. For example, in refs. [Reference Chettibi, Lehtihet, Haddad and Hanchi12, Reference Gasparetto and Zanotto23], 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 [Reference Yang24, Reference Emmerich and Deutz25]. 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. [Reference Saravanan, Ramabalan and Balamurugan26, Reference Saravanan, Ramabalan and Balamurugan27] 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. [Reference Xu, Li, Chen and Hou28], 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. [Reference Shi, Fang and Guo29] 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. [Reference Abu-Dakka, Valero, Suner and Mata30]. In ref. [Reference Huang, Hu, Wu and Zeng31], 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. [Reference Rout, Mahanta, Bbvl and Biswal32] to handle the time–jerk–torque optimal problem for a six-axis robot. In ref. [Reference Lan, Xie, Liu and Cao33], 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. [Reference Wu, Zhao and Zhang34], 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. [Reference Zhang and Shi35], 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. [Reference Cheng, Hao, Wang, Xu and Li36] 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. [Reference Shi, Wang, Ke, Zheng, Zhou, Wang, Fan, Lei and Wu37], 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 [Reference Tian, Cheng, Zhang and Jin38, Reference Blank and Deb39]. 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 [Reference Coello, Abraham, Jain and Goldberg40, Reference Zitzler, Deb and Thiele41]. 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 [Reference Deb20, Reference Sindhya, Miettinen and Deb42–Reference Alotto and Capasso44].
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 [Reference Berger, Bruns, Ehrmann, Haldar, Hafele, Hofmeister, Hübler and Rolfes46], the open-source platform for solving optimization problems PlatEMO [Reference Tian, Cheng, Zhang and Jin38, Reference Tian, Zhu, Zhang and Jin47], or pymoo [Reference Blank and Deb39]. 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 ini to a final location q fin . The equation of motion for the i th joint can be written as follows [Reference Craig1]:
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 ini 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 ini , q fin ) are characterized by null joint velocities and accelerations. Thus, we have
2.1.2 Geometric constraints
First, there are constraints imposed on joint positions:
But, if obstacles are present in the workspace, collisions must be avoided and the following additional constraint will hold during the transfer:
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
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:
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:
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.
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 [Reference Sindhya, Miettinen and Deb42, Reference Talbi43, Reference Jaszkiewicz48, Reference Sindhya, Deb and Miettinen49]. 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 [Reference Tian, Cheng, Zhang and Jin38, Reference Tian, Zhu, Zhang and Jin47]. More precisely, a global–local search strategy based on a high-level relay hybrids model with a sequential execution order will be adopted [Reference Talbi43]. 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).
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 = 0, t 1, …, 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]:
subject to : h ( X ) = 0, g ( X ) ≤ 0, x min ≤ $\boldsymbol{X}$ ≤ x max
where
-
• $\boldsymbol{X}$ = {x 1, x 2, …, 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 1( $\boldsymbol{X}$ ), f 2( $\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]:
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 1( $X$ ), f 2( $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]:
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. 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 [Reference Fotouhi and Szyszkowski9–Reference Massaro, Lovato, Bottin and Rosati11]. According to rel. (1), the inverse dynamic of this robot is as follows:
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.
As an application of the first phase of the proposed approach, the time evolution of the joint position variable (q 1(t), q 2(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:
-
• 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 F1 and F2.
-
• Figure 6 illustrates one selected solution from PF obtained using C3M-PS for which the transfer time is F1 = 1.1092s, and the normalized mean square value of applied torques is about F2 = 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 F1 = 1.002 S (decrease of about 10%), but the corresponding normalized mean square value of applied torques is about F2 = 2.02s (increase of about 90%).
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 [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.
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:
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°/s3, 207°/s3]. 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°/s3, 188.98 °/s3]. 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°/s3 153.41°s3]. 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.
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.
Authors contributions
All the work is done by Taha CHETTIBI.
Financial support
This research received no specific grant from any funding agency.
Competing interests
The author declares no competing interests exist.
Ethical considerations
The submitted work is original and not has been published elsewhere in any form or language.