Multi-objective Trajectory Planning Method based on the Improved Elitist Non-dominated Sorting Genetic Algorithm

Robot manipulators perform a point-point task under kinematic and dynamic constraints. Due to multi-degree-of-freedom coupling characteristics, it is difficult to find a better desired trajectory. In this paper, a multi-objective trajectory planning approach based on an improved elitist non-dominated sorting genetic algorithm (INSGA-II) is proposed. Trajectory function is planned with a new composite polynomial that by combining of quintic polynomials with cubic Bezier curves. Then, an INSGA-II, by introducing three genetic operators: ranking group selection (RGS), direction-based crossover (DBX) and adaptive precision-controllable mutation (APCM), is developed to optimize travelling time and torque fluctuation. Inverted generational distance, hypervolume and optimizer overhead are selected to evaluate the convergence, diversity and computational effort of algorithms. The optimal solution is determined via fuzzy comprehensive evaluation to obtain the optimal trajectory. Taking a serial-parallel hybrid manipulator as instance, the velocity and acceleration profiles obtained using this composite polynomial are compared with those obtained using a quintic B-spline method. The effectiveness and practicability of the proposed method are verified by simulation results. This research proposes a trajectory optimization method which can offer a better solution with efficiency and stability for a point-to-point task of robot manipulators.


Introduction
With the advancement of the times, robotics technology is also developing rapidly, which makes manipulators widely applied in industrial filed. Responding to many practical robotic applications (such as palletizing, labeling, spot welding), the trajectory planning of the manipulator is very importance for accomplishing tasks, which generally involves two key problems, namely trajectory generation and trajectory optimization. The former provides the precondition for the trajectory planning. Meanwhile, the latter is an efficient way to improve the performance of the trajectory and get the most of the manipulator [1].
Trajectory generation is usually to establish a smooth trajectory of a manipulator by means of interpolating between any two given poses. Common interpolation functions include polynomial, spline, Bezier, and NURBS, etc. The joint motion of an industrial robot was divided into accelerated part, constant velocity part and decelerated part, where the accelerated and decelerated trajectories were planned with fourth-order polynomials formed with the property of the root multiplicity [2]. In Ref. [3], the trajectory must pass through a number of given discrete characteristic points. The time-optimal and jerkcontinuous trajectory planning has been implemented under kinematic constraints by combining cubic splines in Cartesian space and septuple B-splines in joint space. Shi et al. [4] adopted a quintic non-uniform rational B-spline (NURBS) to construct a flexible trajectory of a 6-degree-of-freedom (DOF) robot, which can guarantee jerk continuous, and also velocity and acceleration of initial and final point both can be specified. Dinçer and Çevik [5] designed a composite polynomial composed of quadratic Bezier curves and cubic polynomials for the trajectory planning of a 2-DOF parallel mechanism. The composite polynomials provide a smoother transition at the starting and ending points compared to Bezier curves, namely, velocities are zero at the endpoints. Motivated by Ref. [5], a composite polynomial, combined quintic polynomials with cubic Bezier curves, is developed in this article. By the polynomial part, the velocities and accelerations of the actuators at initial and final instance are zero to achieve the stability of start and stop. The better index performance is obtained by adjusting the trajectory constructed by the Bezier curve through optimization algorithm.
In the processing of trajectory optimization, many different technical criteria have been defined to meet the requirements of the task [6][7][8]. For instance, the execution time and jerk are intended for improving the productivity and keeping the trajectories smooth. The energy and torque are aimed at reducing the energy consumption and the load on the actuator of the robot. In addition, with respect to trajectory optimization techniques, evolutionary algorithms, which offer high efficiency, robustness and adaptability, have been widely applied to resolve minimization problems for objective trajectory functions. In Ref. [9], genetic algorithm (GA) was applied to the trajectory planning problem with nonlinear constraints and obstacles to minimize the joint rotation angles of a 2-DOF robot. Lin [10] employed particle swarm optimization (PSO) with K-means clustering to solve the near optimal solution of a minimum-jerk joint trajectory. Only considering a single-objective function may not be suitable for meeting multiple requirements in real-world applications. Currently, in terms of the multi-objective optimization, multi-objective evolutionary algorithms (MOEA) typically utilize non-dominated sorting to provide a number of Pareto solutions for decision-makers rather than converting all objectives into a single-objective function. Thus, this optimization technique has become more popular with researchers. In Ref. [11], the time-jerk trajectory of a robotic manipulator was interpolated in the joint space by means of 5thorder B-splines and then optimized by NSGA-II. In Ref. [12], a multi-objective function, including motion time, dynamic disturbance, and jerk, was addressed by using multi-objective particle swarm optimization (MOPSO) to obtain the high efficiency and safe motion trajectory of a space robot. Marcos et al. [13] combined the closed-loop pseudoinverse method with a multi-objective genetic algorithm (MOGA) to minimize the joint displacement and the positional error of the end-effector. Ramabalan et al. [14] adopted B-spline functions to define the trajectory of a robot manipulator, and the trajectories optimized by NSGA-II and multi-objective differential evolution (MODE) were compared. The results showed that the efficiency obtained by the MODE technique was higher, while the richer the diversity of the Pareto solution was got by NSGA-II.
The main differences among the above trajectory planning methods lie in the processing of the interpolation functions and trajectory optimization models, as well as the selection of interpolation functions and optimization algorithms [15]. Due to the complexity of the trajectory problem of manipulators, it still has improvement space in the accuracy and efficiency of the solution method. Therefore, the proposed composite polynomial is first adopted to construct a point-to-point trajectory in this study. Then, to improve the convergence and diversity of the Pareto optimal front and also the computational efficiency of the traditional NSGA-II, an INSGA-II to obtain the time and torque fluctuation optimal trajectories is proposed.
This article is organized as follows. In Section 2, a composite polynomial curve, by combining cubic Bezier with quintic polynomial, is presented for establishing trajectory optimization model. The three improved genetic operators and INSGA-II are proposed in Section 3. In Section 4, two performance measures are delineated for the Pareto front and the computational efficiency of the algorithms in detail. The numerical simulations are presented with relevant discussion in Section 5. Finally, the main conclusions are outlined in Section 6.

Trajectory Optimization Modelling
The trajectory planning is generally carried out in operating space and in joint space. In terms of the trajectory planning in joint space, it can avoid singular configurations for the robotic arm, but its application has been limited due to the nonlinear relationship between operating space and joint space [16,17]. Moreover, the analytical expressions for forward kinematic solutions of most parallel mechanisms are hard to obtain, and only the numerical solutions can be found. Another method is to perform the trajectory planning in operating space. It is intuitive to avoid obstacles and easy to track the end-effector position and posture [18,19], but the problem of kinematic singularity is difficult to address using such a method. Hence, for the serial mechanism, trajectory planning can be carried out in the joint space if there is no need for obstacle avoidance. In light of parallel and hybrid mechanisms, the trajectory planning problem is handled in the operating space to facilitate analyzing the dynamic performances.

Objective Function
In light of the manipulator, it is expected that the joint trajectory is smooth enough to avoid overlarge mechanical vibration, and reduce travelling time as much as possible to improve productivity. In Refs. [11,20], the time integral was set to be a term of the objective function. Li and Wang [21,22] further took the minimum absolute value of torque fluctuation into account. In this study, two objective functions are involved, namely, the travelling time and the torque fluctuation. The objective functions can be mathematically defined as follows. Minimize: where f 1 denotes the total travelling time, f 2 is the variance of the actuator torque, which is to ensure the stability of the manipulator. τ i (t) and τ i (t − 1) denote the torque of the actuator at former and current instance, respectively. n denotes the number of the robotic joint. It is obvious that the two objective functions constrain each other because of the opposite effects. The reduction in travelling time would lead to the larger torque fluctuation and less smooth trajectory, while reducing the torque fluctuation would result into the longer execution time and lower work efficiency. The trajectory planning inevitably encounters a trade-off between these two objective functions. Therefore, by solving the optimization problem with a multi-objective optimization technique, a set of Pareto solutions can be obtained and provided for decision-makers to select. It should be noted that objective functions can be established for different actual needs.

Constraint Conditions
The kinematic constraints include the limits of angular velocity and acceleration, and the dynamic constraint is mainly the actuator torque. To guarantee the starting and stopping stability of the manipulator, the velocity and acceleration of the actuator are identical zero at the endpoints. The expression in mathematical terms can be written as where sup|·| denotes supremum of parameter. t 0 and t f represent initial and final moment. c d , c v , c a , and c m denote the maximum angle, angular velocity, angular acceleration and torque of each actuator during the entire motion. The equalities describe the initial and final state required for the manipulator, and the inequalities describe the performance of each actuator.
The maximum value of velocity and acceleration in Eq. (2) can be satisfied through determination of the travelling time by the following formula where ψ = dψ/dt , ψ denotes the joint velocity. ψ is the joint acceleration, ψ = dψ/dt.

Composite Curve
In Ref. [5], a new composite polynomial is generated by combining cubic polynomials with Bezier curves based on quadratic Bernstein polynomials. This trajectory planning provides a much lower jerky motion that decreases unwanted vibration. However, the acceleration of the mechanism is not zero at initial and final points by applying the composite polynomial into constructing trajectory, which is unfavorable for the start and stop of the manipulator. Therefore, we develop a new composite polynomial by combining quintic polynomials with Bezier curves based on cubic Bernstein polynomials. The better index performance is obtained by adjusting the trajectory constructed by the Bezier curve through optimization algorithm. By the polynomial part, the velocities and accelerations of the actuators at initial and final movement are zero to improve the stability of start and stop.
A Bezier curve of degree n can be defined in parametric form as where the polynomials B n i (x) are known as Bernstein basis polynomial of order n, n i is a binomial coefficient. P i is the given control point to construct the Bezier curve.
To further improve the start-stop stability of the manipulator, the quintic polynomial can be designed as Hereupon, the composite polynomial can be obtained by substituting the quintic polynomials into Eq. (5) such that x = g (λ). It can be expressed as In Eq. (6), λ denotes the normalized time, for the travelling time T = t f − t s , if we define t = λ·T, the trajectory of the OAJ can be expressed as The trajectory of each actuator can be obtained by applying the inverse kinematics transformation into Eq. (7).
The Bezier curve of Eq. (4) provides a better convergence to the starting and ending points, while the polynomial of Eq. (5) provides a smooth transition in the vicinity of the endpoints. By this method, we tried to exploit the advantages of each polynomial, and the corresponding results would be presented in Section 4.

Proposed Method
The non-dominated sorting genetic algorithm has been established itself as a benchmark algorithm for multiobjective optimization, which was first proposed by Deb et al. [23]. The main contribution is to obtain the Pareto solutions by sorting the dominated relationship among individuals. However, the basic algorithm suffers from a high order of complexity and highly depends on shared parameters. Hereto, in the iterative processing of NSGA-II [24], the shared parameters are replaced with the crowding degree, while the elite strategy is introduced to retain the excellent individuals. It, adopting the fast nondominant sorting method to reduce the computational complexity, has been demonstrated the ability to find a good spread of solutions and converge close to the true Pareto-optimal front. Subsequently, to solve the insufficiency of NSGA-II in dealing with the four or more objectives optimization problems, the reference point method of the NSGA-III [25] was utilized to substitute the crowding degree method in the replace operation, which can perform better in balancing the diversity and convergence of the algorithm. There are only two objective functions in the trajectory optimization problem of this article, so we consider NSGA-II as the benchmark algorithm. However, the selection, crossover and mutation operators of traditional GA are adopted in NSGA-II, which leads to the loss of population diversity and the poor search ability of the algorithm [26]. Moreover, the manipulator is a nonlinear multivariable and strong-coupling system with extremely complex kinematic and dynamic models. To avoid problems such as premature convergence and low convergence speed in the processing of trajectory optimization of the manipulator using the conventional NSGA-II, INSGA-II, integrating three specially designed operators, is proposed to quickly and accurately obtain the optimal trajectory.

Ranking Group Selection
The roulette wheel selection and tournament selection are generally used as GA selection operators. Although their operation mechanism is simple, the process is complicated and requires repeated comparison of the fitness [24,26]. Motivated by Ref. [26], a ranking group selection is used to replace the conventional selection.
The procedure of the RGS is shown in Figure 1. First, the parent population P 0 of size N is randomly initialized based on the constraints of the designed variables, where N is set to a multiple of four. Then the initialized population is sorted into several ranks based on the non-domination sorting. The solutions are assigned fitness equal to the corresponding non-domination levels. Individuals in the first front are assigned fitness value of 1, and individuals in second are given a fitness value 2 and so on. Afterward, the sorted population is uniformly divided into 4 elements in sequence, namely X 1 , X 2 , X 3 , and X 4 . Using the basic concept of combinatorics, there are six cases in which two elements are selected from four elements for pairing, including (X 1 , X 2 ), (X 1 , X 3 ), (X 1 , X 4 ), (X 2 , X 3 ), (X 2 , X 4 ), and (X 3 , X 4 ). The paired groups of individuals formed by RGS are I A = (X 1 , X 1 , X 1 , X 2 , X 2 , X 3 ) and I B = (X 2 , X 3 , X 4 , X 3 , X 4 , X 4 ). In the iterative process, I A is responsible for guiding the population towards the optimal region while I B is responsible for increasing the population diversity.
In Figure 1, P t and Q t represent parent population and offspring population, and P t+1 represents the parent population in next generation. F 1 , F 2 and F 3 denote the different ranks of the population. The individuals in I A and I B that are paired in turn to participate in the crossover can improve the gene diversity of the population and avoid inbreeding, which can promote the generation of high-quality individuals in the procession of gene recombination. Additionally, the RGS is a way to directly calculate the values of the objective functions instead of contrasting looping manner, so the time complexity of the method is small and such the method is easy to implement.

Direction-based Crossover
The simulated binary crossover is commonly adopted in GA, which uses a random way to carry out the gene exchange between individuals. Although the operation mechanism is simple, the method causes degree of blindness. Based on the principle that the better the objective function is, the closer the individual is to the optimal region, and a DBX operator is designed.
Taking two objective functions and two-dimensional variables as an example is shown in Figure 2. The distribution of the Pareto solution set is obtained according to Section 3.1, and the corresponding individuals are assumed to be X 1 and X 2 . The DBX takes X 1 as the center and uses the direction vector d 11 or d 12 as the crossover direction to generate new individuals along random steps. The DBX operator can be mathematically expressed as where i denotes the ith individual, j denotes the variable dimension, the parameter r ij is a unifomly distributed random number in the interval [− 1, 1]. Different from the traditional fixed step crossover, the step size in DBX is randomly generated by the parameter r ij , which expands the search range of the algorithm.
Meanwhile, different rectangular areas are produced by different paired individuals, which indicates that the difference between paired individuals can increase the population diversity to improve the search ability of the algorithm and increase the generation probability of high-quality individuals. It is noted that if the population generated by Eq. (8) crosses the boundary, it will be limited to the boundary to ensure the rationality of the population genes.
where p it denotes the value of the ith individual in tth iteration. p imin and p imax represent the minimum and maximum value of designed variables.

Adaptive Precision-Controllable Mutation
The purpose of introducing mutation in GA is twofold: One is to make genetic algorithm have local random search ability. The second is to maintain population diversity of the algorithm to avoid immature convergence. In Ref. [27], a simple and efficient precision-controllable (8)  mutation (PCM) operator is proposed for exploration and exploitation. On the basis of the Ref. [27], a selfadaptive mechanism is incorporated into the PCM to improve the convergence speed of the algorithm in this article.
where �α = 1 10 Random (p) +1 × (Random(9) + 1) The variable p is the parameter to control the precision in decision space. Function Random (p) can generate a pseudorandom number in the range of 0 to p − 1. If the required search precision is 0.001, the parameter p can be set to 3. The value of random number Random (3) should be in the set of {0, 1, 2}, then the corresponding value ranges of �α from 0.001 to 0.9.
Eqs. (10) and (11) are intended for exploitation, while Eqs. (12)-(15) are designed for exploration. The operator can effectively explore and exploit the decision space, and its computation process is simple, and precision is controllable. However, the mutation operator in Ref. [27] has not sufficiently utilized the potential information of the contemporary population, which can be used for the adaptive selection of exploitation or exploration.
A common adaptive adjustment method is to use the information of the objective function value to adjust the mutation strategy. For the trajectory planning problem of the manipulator, the real Pareto front cannot be obtained in reality. Therefore, the distribution value of population instead of objective function can be used as a parameter to select exploitation and exploration. In this study, the ratio of contemporary population space to decision space is used to determine the mutation strategy of the individual, which can promote rapid and stable evolution of the population.
In the early iterations, the differences between the individuals are larger, so the exploration is selected to ensure the diversity of the population and avoid the algorithm falling into local optimality. In the later iterations, the population gradually tends toward the region of the optimum, and the differences between the individuals are smaller, hence the exploration is selected to keep the excellent individuals to improve the search effective. The APCM operator can be expressed as where x C i max and x C i min are maximum and minimum value of the individuals in contemporary population, x i max − x i min denotes the magnitude of decision space.
Compared with the mutation operator in Ref. [27], adding the adaptive adjustment of the mutation operator can promote the balance of local search and global exploration capabilities, thereby making the Pareto boundary distribution better.

Overall Algorithm
The flowchart used the INSGA-II for the trajectory planning of the manipulator is shown in Figure 3. Initializing randomly the parameters within the threshold value gains the initial trajectory curves (Eq. (7)), and the storage of these trajectory is performed, then calculating the objective functions for each chromosome. The first generation population performs non-dominated sorting to find a set of Pareto front (PF), and the population is sorted by the crowding distance. Afterwards, a new parent population is generated by RGS, DBX, and APCM operator, and the parents and offspring are combined to form a population of N individuals according to the elite strategy. It is continuously judged whether it reaches the number of iterations, and the objective function of each trajectory is compared. Finally, the Pareto solution of the objective functions is obtained after the iteration and the corresponding designed parameters are output.
Compared with conventional NSGA-II, the combination principle to construct a selection operator is used in the proposed INSGA-II, which can avoid repeated comparison of the fitness between individuals to improve the convergence speed of the algorithm. (16)  DBX can expand the search space and increase the generation probability of high-quality individuals, thereby improving the search ability and convergence speed of the algorithm. The local random search ability of the APCM can accelerate the convergence to the optimal solution, and the exploration strategy of the operator can expand the search space to keep the population diversity.

Performance Measures for INSGA-II
Given a set of solutions by Section 3, but in some cases, the weight cannot be determined by a decision maker due to insufficient information related to the different criteria. In that situation, we offer a strategy. First, the performances of the INSGA-II are evaluated according to the convergence, diversity and computational efficiency, and then the fuzzy comprehensive evaluation of the solution set is adopted to determine the optimum solution for decision markers.

Performance Evaluation Index
As for the performances of the multi-objective optimization algorithms, inverted generational distance (IGD) and hypervolume (HV) are very popular for comprehensively measuring the convergence and diversity of algorithms [28]. Meanwhile, the proportion relation between total number of evaluations and total CPU time is used to test the algorithm efficiency [29]. The three metrics can be mathematically expressed as follows where P* indicates a set of points uniformly sampled over the true PF, and S is the set of solutions obtained by an MOEA. dist(x, S) denotes the Euclidean distance between the closest individual from x to S. |P*| is the cardinality of set P*. The smaller IGD value indicates that the set S is closer to the entire PF, and thereby the convergence and diversity are better.
where r * = r * 1 , r * 2 , · · · , r * m is a reference point in the objective space that is dominated by all solutions in a PF approximation S. VOL(•) is the Lebesgue measure. HV metric measures the size of the objective space dominated by the solutions in S and bounded by r * . The larger the HV value, the closer S is to the entire PF.
where OO stands for the optimizer overhead. T Total denotes the total CPU time taken, and T PFP denotes the time taken for pure function evaluations. The lower OO metric corresponds to the higher efficiency of the algorithm.
It should be noted that in IGD, an average minimum distance is calculated from each point in the true PF to those obtained by an MOEA. In the processing of calculating IGD, since without any priori PF shape knowledge, all the non-dominated solutions are used as the reference points [30].

Fitness Evaluation
Fuzzy comprehensive evaluation is one of the effective decision-making methods for objectives affected by various factors, which adopts a fuzzy membership function to describe the fitness factor of an objective function [11,29]. For the objective function minimization problem, the fuzzy membership function can be expressed as where f i (j) represents the objective function, i denotes the number of the objective function, and j denotes the jth solution at the PF. f max and f min are the the maximum and minimum value of the objective function. According to Ref. [11], the synthetical membership value can be written as The larger synthetical membership value indicates the better fitness of the Pareto solution. The highest synthetical membership value is 1, which can be considered as the most satisfactory solution for the decision-maker.

Numerical Example
The purpose of conducting simulation is to verify the search capability and convergence speed of the proposed INSGA-II as well as validity and competency of the composite polynomial approach for creating trajectory. In this section, taking a serial-parallel hybrid manipulator as instance [31], as shown in Figure 4, the symbols in the figure are given a detailed introduction in the Ref. [32]. To facilitate analyzing the dynamic performances of a hybrid manipulator, the trajectory planning problem is handled in the space of the output angle of joint moving platform (OAJ) [32]. It is noted that the mapping of the trajectory between the joint space and the OAJ space can be obtained by applying the inverse kinematics transformation of each joint of the hybrid manipulator, and the OAJ space can be transferred to the operating space through the forward kinematics analysis.
The maximum value of each OAJ velocity and acceleration can be obtained based on the search method [33] by combining the actuator velocity and acceleration boundaries in Eq. (2) with the workspace of the manipulator [21].
where ψ,ψ,ψ ∈ R n , ψ denotes the n-vector of OAJ, n is the number of DOF of OAJ, j denotes each OAJ. J o max represents Jacobian matrix of the manipulator corresponding to the maximum value of each OAJ velocity.
Then, we start by analyzing the performance of the several algorithms for the trajectory optimization of a point-point motion, including the proposed INSGA-II, MO-INSGA-II [34], success history-based adaptive multi-objective differential evolution with whale optimization (SHAMODE-WO) [35], IMOPSO [36], many-objective evolutionary algorithm based on decomposition with random and adaptive weights (MOEA/ D-URAW) [37] and IMODE [38]. In a second phase, the composite polynomial with the quintic B-splines approach are compared to evaluate its effectiveness.

Comparison with MO-NSGA-II, SHAMODE-WO, IMOPSO, MOEA/D-URAW and IMODE
Given position and posture of the hand at initial and final instance, including one starting point and eight ending points, the OAJ trajectories are parameterized here by composite polynomial functions with four nodes uniformly distributed along time scale ( Figure 5). There are sixteen unknown parameters, where two intermediate adjustable nodes contain fourteen, ending point only includes one by applying inverse kinematics transformation, and the travelling time is one of them. These To validate the superiority of the proposed INSGA-II, its performances with some state-of-the-art representatives are compared from different categories of multiobjective algorithms. All non-dominated solutions of the trajectory optimizations, from the starting point to the eight ending points, offered by INSGA-II over 100 runs, are compared to that of MO-NSGA-II, SHAM-ODE-WO, IMOPSO, MOEA/D-URAW and IMODE in terms of IGD, HV and OO, and then these experiment results are gathered for statistical analysis. It is noted that all objective functions are normalized by adopting the min-max standardization method to have a same range, which can avoid the function with largest range would dominate selection.
The initialization parameters for NSGA-II are as follows: the population size is 100, the generation number is 80. Mutation probability is 1/16, which is selected as 1/n (where n represents the number of variables) proposed by Deb [25]. For constraint-optimization problems, the distribution indexes for real coded crossover and mutation operators are 20 and 100, respectively. The values of the parameter that have been used in SHAMODE-WO technique are as follow: the population size is 100, the generation number is 80, the historical memory of scaling factor is 0.5, the historical memory of crossover ratio is 0.5, the memory index is 1 and the memory size is 5. In overall tested experiments, IMOPSO was run using the parameters as follows: the population size is 100, the generation number is 80, the jump improved operation mechanism number is 100, the disturbance rate range is [0.1, 0.3]. For MOEA/D-URAW, the population size is 100, the generation number is 80, the historical memory of scaling factor is 0.5, the historical memory of crossover ratio is 0.5, the memory index is 1 and the memory size is 5. For IMODE, the population size is 100, the maximum number of iterations is 80, the crossover probability is 0.1, the scaling factor is 0.5, the size of initial Pareto front approximation is 100, the number of points desired by the decision maker is 100 and the selection parameter is 0.1.
The mean and standard deviation values (SD) of all the instances are shown in Table 1 (The best results for each index are marked in bold). Demonstrated in mean and SD of the IGD evaluation results, INSGA-II finds better solutions, which has superior values in all test problems when compared to related works while its IGD remains approximate to zero. In most of the test problems, the INSGA-II performs better than other related methods in the HV evaluation results. The two evaluation results indicate that the convergence and diversity of the nondominated solutions obtained by INSGA-II performs better than that of other related methods. However, the experiment results of the INSGA-II are not satisfactory in the OO evaluation. As for the OO evaluation results, the calculation efficiency of all the test instances addressed by IMODE is the best.
In order to intuitively reflect the performance of each algorithm, the corresponding boxplots ( Figure 6) is drawn by synthesizing the results in Table 1. The IGD and HV evaluation results clearly show that the INSGA-II can perform exceptionally in solving the problem of manipulator's trajectory planning, which demonstrates its convergence and diversity are better than other algorithms. MO-NSGA-II takes second place in the convergence and diversity performances, but the computational efficiency of the INSGA-II and MO-NSGA-II are undesirable. Additionally, although the convergence and diversity performance of IMODE technique is not good, it gives minimum OO thereby it is the better one for a multicriterion to obtain a best optimal solution trade-off very quickly.

Comparison of Trajectory Planning Methods
Taking the point-to-point mission with FPH1 as an example, to gain the better designed parameters of trajectory optimization, the non-dominated solutions obtained by the four algorithms are all taken as candidates for decision-makers. As shown in Figure  In order to be able to compare the results yielded by the two trajectory planning methods, the travelling time is consistent with the optimized results. In addition, the velocity and acceleration of the starting and ending points are set as zero. The profiles of velocity and acceleration for each actuator created using quintic B-splines are shown in Figures 8 and 9.
Likewise, the trajectory of each OAJ created by using the composite polynomials can be obtained by Table 1 Mean and standard deviation values of GD, SSM and OO between INSGA-II, MO-NSGA-II, SHAMODE-WO, IMOPSO, MOEA/D-URAW and IMODE over 100 Table 2, and it can be noticed that the results yielded by the approach described in this article are better than those provided by the approach [21] with respect to the maximum values of velocity and acceleration. It is well known that the lower maximum velocity provides an advantage because lower velocity extends the life of the actuator. Meanwhile, lower acceleration profiles decreases the noise in the mechanism and increases the mechanical life by reducing wear.
To further demonstrate the superiorities of the proposed method, the torque fluctuations are calculated by the quintic B-spline and composite polynomial approach. The reduction of the torque fluctuation is 33.47%, as shown in Figure 12. Consequently, the manipulator can work with higher stability via the proposed method.

Conclusions
A new methodology for optimal trajectory planning has been described in this article. The composite polynomials are adopted to construct the trajectory of each OAJ and the trajectory is optimized with INSGA-II technique. The objective functions take into account both the travelling time and the torque fluctuation along the whole trajectory.
(1) A new composite polynomial is created by combining quintic polynomials with Bezier curves based on cubic Bernstein polynomials. By the Bezier curve part, the convergence to the starting and end- The simulation results demonstrated that well-converged and well-diversified non-dominated solutions can be obtained by INSGA-II, but the efficiency is lower than that of IMODE. (4) Using the synthetical fuzzy membership function to obtain a trade-off for decision-users, the trajectory of the OAJ constructed by composite polynomials compared in the velocity and acceleration with quintic B-splines. The former velocity and acceleration are lower, which increases the mechanical Figure 11 The acceleration profile of each joint using composite polynomials  life by reducing wear. Moreover, the reduction of the torque fluctuation is 33.47%, thereby ensuring higher motion stability of the manipulator.
Future work will aim to reduce the time complexity of the INSGA-II to improve the calculation efficiency so that the optimization method can be used in the realtime trajectory planning for the manipulator.