Multi-Objective Trajectory Planning for Slung-Load Quadrotor System

In this article, multi-objective trajectory planning has been carried out for a quadrotor carrying a slung load. The goal is to obtain non-dominated solutions for path length, mission duration, and dissipated energy cost functions. These costs are optimized by imposing constraints on the slung-load quadrotor system’s endpoints, borders, obstacles, and dynamical equations. The dynamic model of a slung-load quadrotor system is used in the Euler-Lagrange formulation. Although the differential flatness feature is mostly used in this system’s trajectory planning, a fully dynamic model has been used in our study. A new multi-objective Genetic Algorithm has been developed to solve path planning, aiming to optimize trajectory length, mission time, and energy consumed during the mission. The solution process has a three-phase algorithm: Phase-1 is about randomly generating waypoints, Phase-2 is about constructing the initial non-dominated pool, and the final phase, Phase-3, is obtaining the solution. In addition to conventional genetic operators, simple genetic operators are proposed to improve the trajectories locally. Pareto Fronts have been obtained corresponding to exciting scenarios. The method has been tested, and results have been presented at the end. A comparison of the solutions obtained with MOGA operators and MOPSO over hypervolume values is also presented.

The use of quadrotors is increasing day by day. One of the tasks quadrotors are employed recently is load transportation. One can equip a quadrotor with a gripper mechanism or connect the load mass to a quadrotor with some wire (slungload) to transport a load. The gripper mechanism adds to the quadrotor's total mass [2]; that is why the slung-load method is preferred. Besides, the gripper mechanism needs additional controllers and consumes energy while doing its job, which means it is more expensive to carry a load. Alternatively, it is preferred that the quadrotor takes the load with a tied rope or wire. In [1], this type of transportation task was carried out with unmanned helicopters; the helicopter-like systems' mathematical models and controller designs were examined for different wire connections.
It is possible to examine the literature on manipulating slung-load unmanned aerial vehicles, mainly in two categories. Studies in the first category are concerned with designing controllers for tracking desired payload trajectories [5], [7], [8], [9], and [17]. On the other hand, studies in the second category are trajectory planning studies to optimize a particular objective function [3], [6], [15], [18], and [19]. Our study belongs to the second category. For an extensive summary of load transportation studies using quadrotors, a recently published survey [3] can be investigated.
We think that in trajectory planning, usually, there is more than one criterion involved with the problem either directly or indirectly [29], [30]. There are studies to optimize objective functions such as minimum swing angles of load [6], [15], and [16], minimum mission time [24], and minimum energy [21], [22]. In [23], an objective function emphasizing the tradeoff between mission time and dissipated energy were investigated. Interestingly, there is no study aiming to optimize multiple objectives simultaneously to the best of our knowledge. It is one of our research contributions. The trajectory planning for a quadrotor slung-load system is investigated under numerous goals such as mission duration, path length, and dissipated energy.
Although there are many studies [5], [12]- [14] that advocate doing multi-objective path planning, most of them implement scalarization (weighted sum of multiple objectives are considered). For example, [5] aims to minimize the mission completion time, which is the sum of the time to set up a communication and find a target. They do not consider these two times separately; instead, they assign weights to prioritize. An online trajectory planning algorithm is proposed in [16] and [17] to minimize the load swing and target positioning error. Although there are multiple goals, the problem is not approached as multi-objective optimization; instead, the problems are handled separately.
Studies utilize the multi-objective optimization approach to solve trajectory planning problems to optimize path length, path safety, and path smoothness on the grid map without considering the vehicle dynamics constraints [26]. In [3], the firefly algorithm has been proposed to optimize objective functions [26]. Although a different application is discussed in [35] than the one we have covered in this study, a framework is proposed that makes good use of the concept of multiobjective optimization. Our study, conceptually similar to [35], aims to optimize more than one objective simultaneously in the trajectory planning problem for quadrotor slung-load system makes the problem even more difficult.
Quadrotor slung-load system is an under-actuated system with eight degrees of freedom and four control inputs (i.e., it is challenging to obtain a trajectory plan even for only one of the objectives mentioned above). Furthermore, simultaneously optimizing more than one objective function makes the problem even more difficult.
Our aim here is to find out how to create a solution to optimize the length of trajectory, mission time, and energy consumed during the mission simultaneously. Intuitively, it is easier to understand that system will spend more energy on the fastest trajectory. However, it is helpful to make a few explanations regarding the relationship between "mission time and path length" and "path length and energy spent" objectives. The fastest trajectory may not be the shortest one because there are obstacles in the environment. It may be necessary to maneuver too much on the shortest trajectory. Since it may be required to accelerate and slow down during these maneuvers, we cannot say that the shortest trajectory is the least energy expenditure.
Even though these objective functions may not be completely conflicting, they are not entirely supportive of others. For example, as shown in Section 5, Figure 14, minimum length, minimum energy, and minimum time trajectories are different. In addition, a slight change in consumed energy immediately resulted in entirely different paths. Usually, there is no single solution that simultaneously optimizes each objective. Alternatively, it is possible to combine the individual objective functions into a single objective using suitable weights. Here the critical question is how one should select these weights. In this case, we will have a single optimum solution rather than a set of "non-dominated" solutions that can be examined for further tradeoffs (i.e., a set of non-dominated solutions is found corresponding to all possible weight selections). The collection of these nondominated solutions is called the Pareto optimal set or Paretofront [28]. In this set, solutions are non-dominated (i.e., none of the objective functions can be improved in value without degrading some of the other objective values) concerning each other. While moving from one Pareto solution to another, there is always a certain amount of sacrifice in at least one objective(s) to achieve a certain amount of gain in the other(s). Pareto optimal solution sets are often preferred to single solutions because they can be practical when considering reallife problems.
A solution's superiority over other solutions is easily determined in the single-objective optimization problem by comparing their objective function values. In multi-objective optimization problems, the quality of a solution is determined by the importance of the optimization criteria. By investigating the Pareto-front, it is easy to observe how different objectives change as one moves on the Pareto-front. On the other hand, if one optimizes a weighted sum of objectives, only a point on the Pareto-front is found; one does not know what happens to the overall objectives if one of the weights is slightly perturbed. There are many disadvantages to solving multi-objective problems with classical optimization methods [28]. Genetic Algorithms (GA) are frequently preferred in the solution of multi-objective optimization problems [27]. Our study used a GA-based algorithm by customizing the genetic operators for the problem at hand.
As stated in [36], where constrained path planning problem is examined, Evolutionary Algorithm (EA) is of great interest in solving such problems. Moreover, while algorithms such as SQP, IPM, which require numerical gradients, only guarantee convergence to the local optimum, heuristics are generally more likely to find the optimal result [33], [36], [38], [39]. In [33], it is underlined that PSO (Particle Swarm Optimization) requires less number of function evaluations among EAs. In [36], it is stated that direct application of MOPSO [45], [48] to the problem is lacking in finding better values for constraint compliance and objective functions. For this reason, a biased search towards the feasible region is proposed in [36]. Similar to the study we carried out in this article, modifications were made to the classical algorithm. In this context, in our study, not only a biased search application towards a feasible region but also an initial population has been generated in this way. In our study, we used only the main structure of GA as a base since it provides a more suitable structure than other EAs for describing the problem-specific operators. In [29] and [30], it is stated that it is vital to use problem-specific (mutation-like) operators to reach a global optimum when dealing with these challenging problems. We believe that our approach yields many original contributions; 1. Trajectory planning as a multi-objective optimization problem for competing costs (energy, time, and length) was solved using a slung-load quadrotor system for the first time (to the best of our knowledge).

A new and successful Genetic Algorithm was
developed for this problem (details of the algorithm are presented in Section 4). In trajectory planning problems, optimization variables are implicitly dependent, and the degree of dependence increases as one moves along the trajectory. The proposed approach almost completely removes this dependency between optimization variables. 3. Problem-dependent mutation-like operators were developed to avoid obstacles and to improve cost vectors locally.
The rest of the paper is organized as follows; Firstly, the problem definition is given in Section II. Secondly, a dynamic model of a slung-load quadrotor system is described in Section III. Thirdly, the details of the proposed trajectory planner are given in Section IV. Section V contains the resultant trajectories of the slung-load quadrotor system and Pareto Front. Performance comparison with classical GA and PSO algorithms and the analysis of gaps in Pareto Front are given in Section VI. Then, all sections are concluded in Section VII.

II. PROBLEM DEFINITION
Problem definition and associated assumptions are given below. Here we start with the assumptions first:  The load is a point-mass,  The cable is massless, and it is not stretchable,  The cable is attached to the center of mass of the quadrotor,  The cable is under tension throughout the mission.
Next is the optimization problem we are supposed to deal with: where 1 ( ), 2 ( ), 3 ( ) are the objective functions corresponding to the length of the path, total mission time, and energy consumed during the mission.
where is the number of points along a trajectory and can be calculated = 100 * . The time step size of the simulation is 0.01 s.
The duration of the mission time can be determined at the end of the trajectory construction procedure ( ).
The vector u consists of the control inputs (i.e., control forces and moments acting on the quadrotor). The number of time instances these control inputs are supposed to be applied on the quadrotor. The vectors LB and UB are the lower and upper bounds of u. = \ ∈ ℝ 3 is obstacle-free configuration space ( Figure 1). We also impose the constraints on roll ( ) and pitch ( ) angles of the quadrotor. The initial and final points of the trajectory are given as input to the problem. It should be noted the mission duration (i.e., ) is a free variable.

III. DYNAMIC MODEL of a SLUNG-LOAD QUADROTOR SYSTEM
In this section, the 3D dynamic model for the slung-load quadrotor system will be described. We start with the coordinate system definitions in the mathematical model, as shown in Figure 2. The mathematical model used in this study is borrowed from [9], and the detailed derivation of this model can be found there. This model is differentially flat, and the flat outputs are the position of the load and the yaw angle of the quadrotor. Unlike [9], the full 3D dynamic model is preferred over the flat model since we imposed constraints on the quadrotor's location and attitude. Additionally, different from previous studies [5], [6], and [18], we do trajectory planning for both the load and the quadrotor to avoid obstacles. The state vector is; 18 where is the position vector of load and is the unit vector from the quadrotor to load in the inertial frame. The input vector is = [ ] where is the lifting force acting on the system and , , are roll, pitch, and yaw moments, respectively (see Figure 2).
where , ( = 1, 2, 3, 4) is the angular velocity of rotor j, aerodynamic force constant, aerodynamic moment constant. The mathematical model equations are given below: Where R is the rotation matrix from Body Frame to Inertial Frame, is the angular velocity of the quadrotor. I is the inertia matrix, and and are the mass of the quadrotor and the load, respectively. ̂ is the skew-symmetric matrix described in [41]. The energy equation, a function of rotor speeds and accelerations given in [20], has been used to calculate the quadrotor's energy consumed during a mission.

IV. TRAJECTORY PLANNER
As stated earlier, the trajectory planning problem discussed in this study is expressed as a multi-objective optimization problem. A set of "non-dominated" [28] solution vectors are to be obtained. It is advisable to represent the Pareto front by a high number of non-dominated vectors. On the other hand, we expect them to be evenly distributed throughout the front. The nature of Genetic Algorithms (GA) is very suitable for solving multi-objective optimization problems [2], which is why a modified GA is developed and used in this study. As emphasized in [32], it is very effective to handle the problem using several stages in constraint trajectory planning problems.
The structure of a chromosome is shown in Figure 3. It consists of slots, each of which defines a maneuver. Each slot VOLUME XX, 2017 9 contains control inputs (u = [v NoS]) and a variable (NoSi i = 1, 2, …, N) that keeps information about how many time steps these control inputs are to be applied. The maximum number of maneuvers is constant. Thus, while the chromosome length remains constant, the total mission duration is not persistent. We present some of the parameters we use in the scenarios below;  The size of the Population is 150.  The length of each Chromosome is at most 500 (i.e., 100 slots x 5 variables).  The step size in our simulations is 0.01 seconds.  The maximum number of generations is 1000. Reasonable initial estimates can produce better solutions with faster convergence if we have prior knowledge of the problem or can be produced at a low computational cost [46], [47]. The method of using good initial estimates is called seeding. In [46], five different algorithms were tested on 48 optimization problems using different seeding methods, and their performance was examined. It has been shown that the seeding process based on prior knowledge significantly reduces the computational cost. We did not randomly generate the initial population. Instead, it is obtained by using waypoints obtained by removing some randomly generated points out of the obstacle zone that does not violate the obstacles. In Phase 2, a structure was created between these waypoints so that the chromosome length gets longer and longer to reach the endpoint goal. In this way, when we come to Phase 3, not all individuals in the initial population have the same fitness value. In cases where the initial population is randomly generated, as noted in [28], nearly all MOEAs evaluate all solutions of the first non-dominant front equally, assigning a similar or similar fitness to all individuals in the initial population. Neither of these solutions has a choice advantage. Naturally, it will take longer for the EAs to advance the search towards the Pareto optimal region. This search is difficult to achieve in the problem we consider, especially in complex scenarios, due to compelling dynamic constraints. Fortunately, the initial population generation technique in the proposed method greatly reduces the probability of not finding a solution in cases where the solution is very difficult to find.

A. Phase-1 Randomly Generated Waypoints
To find non-dominated trajectories more quickly, we found waypoints between the starting and the final points without obstacle violation (Figure 4). These waypoints have been found by generating random points in polar coordinates. In particular, first, a radial line making an arbitrary angle to the radial line connecting the starting point [0 0 3] T and final point [40 40 10] T is generated. Next, points are created randomly along this line. If it exists, obstacles can be handled by moving these randomly generated points into the feasible region to have a pre-defined margin with the obstacles (Figure 4).
Although about 10 -15 solutions are usually sufficient in real applications, the population size must be chosen much larger due to the nature of EA algorithms. Population size should be determined according to the complexity of the problem and the number of decision variables [44]. Population sizing is an iterative process; however, the graph based on the randomly generated population, given in [28 pp.403], can be used as a good starting point to determine the population size. The initial population in our study has been formed based on the data obtained at the former stages, not completely random. Therefore, the population size has been chosen as 150, aiming to find a nominal number of trajectories between 10 and 15.

B. Phase-2 Initial non-dominated pool construction
One of the main difficulties in trajectory planning is that any change in a specific point significantly impacts the trajectory's remaining parts. In other words, optimization variables are implicitly dependent, and the degree of dependence increases as one moves along the trajectory [29]. The following idea can reduce this dependency; the trajectory is constructed stepwise. At each step, new waypoints are considered ( Figure 5 and Figure 6). Equivalently, new slots are added to each Chromosome's end at each step, depending on the number of waypoints considered. Even though the solution obtained in previous steps (except the first step) is allowed to change within the algorithm's structure, it is observed that these changes are minimal. The operations outlined above can be repeated for all trajectory candidates, for all the waypoints on each trajectory are repeated, and finally, one obtains the initial non-dominated chromosome pool.

C. Phase 3 Multi-objective optimization
There are several multi-objective GA optimization codes in the literature. The most well-known of these are NSGA and NSGA-II [2]. NSGA-II is better than NSGA in terms of computational cost and elitism. NSGA-II uses a fast nondominated sorting technique, an elitist-keeping technique, and a new niching genetic operator without parameters [2]. In this study, we have implemented a multi-objective optimization algorithm. This algorithm uses only the sorting structure of the NSGA-II to find a non-dominated solution. An important contribution of this study is that the chromosome length is variable, and the introduction of Phase-1 and Phase-2 will result in a very fast algorithm.

D. Problem specific operators
One should be careful in applying the crossover operation; the crossover points are selected between the slots. Regularization Operation (RegOp) -Since the ranges of force and moment inputs are different (between 20 N and 35 N for force, ±0.5N. m for moments), it is not meant to change the force and moment creating values in the same range to mutate these inputs. For this reason, we perform the mutation operations for force and moments separately. The force and VOLUME XX, 2017 9 moment values have been changed per the maneuver to be made.

The number of
Step Mutation (NoSMut) -As the NoS values are integers, a different mutation operator is also defined. By mutating these values, one can change the total mission duration. Again, this operator affects the dissipated energy by the quadrotor.
Combined Mutation (ComMut) -In this mutation operation, NoS is increased/decreased, and the objective functions are minimized by increasing or decreasing the force and moment values accordingly. Figure 7 presents this operator's effect; the slung-load quadrotor positions are plotted during that particular time interval. The operator decreases the NoS and increases the applied force simultaneously. While the time spent in this part of the trajectory decreases, consumed energy increases, but the path length remains constant for this part.

34:
Call CheckForConstraintViolation #Discard paths if there is obstacle or border violation.

35:
Call CheckAngleConstraints(v) #check for the angle of the p vector and penalize if the angle values exceed the defined limit.

50:
Call CheckForConstraintViolation #Discard paths if there is obstacle or border violation.

51:
Call CheckAngleConstraints(v) #check for the angle of the p vector and penalize if the angle values exceed the defined limit.
The reason for changing the moments in these intervals is the dynamic constraints of the vehicle. The leaning angles of the vehicle must be within certain limits. Therefore, even if we give higher torque values, the vehicle will not fulfill these commands. We could have given it in lower steps, but it takes longer to find the result this time. For Force, choosing the steps too small increases the time required to achieve results. We could have given larger force values, but in this case, when there is a velocity component in the horizontal, it increases the oscillation of the load by performing sudden accelerations, thus resulting in a long time to get the result. It is not very important to choose 4% or 8% of the proportional values of the operators. Thanks to the robustness that the operators we define provide to GA, the effect of these choices on performance is almost equal. Increasing the ratio of these operators too much would be against the GA structure.

V. RESULTS
The simulations of case studies investigated in our approach were made using MATLAB/Simulink. The state equations (differential equations) and constraints have been taken care of implicitly by the Simulink solver. The simulations run on the personal computer, with processor configuration: 3 GHz, Memory (RAM): 32 GB. The overall algorithm takes 3612.145 s for the 1 st Scenario and 3686.104 s for the 2 nd Scenario. But the dynamic simulation of the quadrotor slungload system's model takes about 75% of the total time. The method described in the previous sections has been tested and associated with Pareto Fronts for two exciting scenarios. The trajectories corresponding to some interesting points on Pareto Fronts have been given. The trajectories shown in these figures are essentially just one element of a trajectory set. Trajectories, which are in the same set, are similar in general shape. Since there are minor differences between trajectories, we prefer to present only one of them; giving the others would not make much sense.
The region for the first Scenario is shown in Figure 8. The borders of this region are defined as; [0 0 0] ≤ [ ] ≤ [50 50 32] . The region boundary bounds the obstacles' height; the vehicle cannot fly over the obstacles. We try to construct a relatively complicated region of interest to get more interesting Pareto Fronts at the end.
The quadrotor's starting point is [0 0 3] T, and the rope's length is 1 meter. The final point of the quadrotor is [40 40 10] T in meters. In Figures 8 and 9, obstacles are shown in blue, borders of the region are shown in yellow, the quadrotor is in black, the load is in blue, and the rope is in red.
In Figure 9, one can observe the enlarged view of the trajectory of the slung-load quadrotor system. One can observe in Figure 10 that the Pareto Front is a disjoint set of points in ℝ 3 consisting of two space curves. Please note that the objective function values in these figures are normalized to visualize the Pareto Front better.
In Figure 11, the projection of the Pareto Front is taken onto the Time vs. Energy plane. When the Path Length is removed from the multi-objective cost, the resultant Pareto Front in ℝ 2 becomes as shown in Figure 12.
In Figures 10 to 15, the Pareto Front in ℝ 3 associated with the multi-objective optimization problem in Scenario-1 and its projections on different coordinate planes in ℝ 2 can be observed. This Pareto Front has been obtained according to path length, total mission time, and dissipated energy minimizations. We have performed a similar analysis when the Pareto Front's projection is taken onto the Energy vs. Path Length plane ( Figure 13). If the scenario duration is removed from the multi-objective cost, the resultant Pareto Front in ℝ 2 becomes as shown in Figure 14. It is interesting to observe the jump on the projected Pareto Front in this figure. The actual trajectories corresponding to the points (Pareto Point-1, Pareto Point-2 in Figure 14) at the beginning and the end of this jump are shown in Figures 15 and 16, respectively. It is important to note that the solution curves in Region-1 in Figure 14 are similar to those in Figure 15. Similarly, the solution curves in Region-2 in Figure 14 are similar to the curve in Figure 16.  The 3D plot of the Pareto Front points is shown in Figures  10, 11, and 13. Note that they are not Pareto Front surfaces (contrary to their representations in those figures). These figures have been obtained using the MATLAB function "scatteredinterpolant".

FIGURE 17. Resulting trajectory corresponding to the non-dominated solution in terms of path-length (3D view).
We do not present the trajectory corresponding to Pareto Point-3 explicitly because it is very similar to the one corresponding to Pareto Point-2. The trajectory corresponding to Pareto Point-4 is illustrated in Figure 18 (its 3D version is shown in Figure 16).
In Figures 20 and 22, one can observe the global optimum solution corresponding to energy and time, respectively. The 3D representations of the paths in Figures 20 and 22 are given in Figures 19 and 21. Scenario-2 aims to construct obstacles so that jumps in the resultant Pareto Front and the corresponding differences in the solution trajectories will be more distinctive. For example, in the Pareto Front for the 2 nd Scenario (Figure 23), the jump points are indicated as has been done earlier. In  Figures 25, 26, and 27, one can observe the global-optimum trajectories corresponding energy (Pareto Point-1), pathlength (Pareto Point-2), and scenario duration (Pareto Point-3), respectively. Additionally, the 3D representation of the solution corresponding to the minimum length is shown in Figure 24.  Since we normalized the objective function values between 1 and 2, it does not give the maximum and minimum values.
Still, the mean and standard deviation values are given in Table 2 for the second scenario.   It can be seen that the quadrotor slows down near the obstacles (for example, see Figures 26 and 27). The effect of this slowdown on the time objective function is obvious. However, due to the possible accelerations and deceleration before and after this slowdown, the effect of the energy objective function may be different. We can compare this to the fuel consumption in traffic; while the short-term red light will slightly prolong our total driving time, the main effect will be on fuel consumption due to stop-go or slow-speed movements. The positions of the obstacles with respect to each other are also factors that directly affect the path-length objective function.  There is no study in the literature regarding the region of interest complexity. We have proposed a method to measure this complexity in our previous work [29], but this problem is still open to solving. From this point of view, we can say that obstacles are the most important element in the trajectory planning problem. While determining the goal function, either a knowledge about the complexity of the region in it should be produced or, as we suggest in our study, VOLUME XX, 2017 9 different goal functions should be looked at together. This study aimed to examine how different objective functions change when different regions are in question.

FIGURE 28 Resulting trajectories obtained using MOGA (yellow), MOPSO (blue), and the proposed method (pink).
We think it is necessary to emphasize Scenario-1 and Scenario-2 covered in this study are pretty challenging. Therefore, finding a Pareto Front with sufficient points and good spread with classical GA operators (i.e., MOGA) or MOPSO. This situation is investigated in a less challenging scenario (Scenario-3) using classical MOGA and MOPSO. We also have tried to solve multi-objective problems given in Scenario-1 and Scenario-2 by increasing the maximum number of generations to ten times that of our algorithm. They could only find trajectories given in Figure 28; they could not find the paths advancing through the obstacles as given in Figure 29. It cannot find the paths advancing through obstacles as given in Figure 26 and Figure 27. We tried to exemplify this situation through a less challenging scenario (Scenario-3) using MOGA and MOPSO. In Figure 28, one can see from this figure; the slung-load quadrotor system reaches the final point by passing as close as possible to the obstacles thanks to the problem-specific operators. So, this solution dominates the solutions obtained using MOGA and MOPSO. In MOGA and MOPSO, it can find the trajectories that do not require too many maneuvers. Still, these two classical methods cannot solve difficult trajectories, such as the route given in Figure 29.
We optimized the controller design for the Quadrotor slung-load system in the previous study [49]. Within the scope of this study, it is checked that the results are insensitive to the quadrotor controller's parameters as long as they are within 20% of the nominal parameters.
Underlining the situation that the MOGA and MOPSO cannot find some trajectories that our proposed method can find, we compare the computation time required for an iteration of our method with the computation time required for an iteration of other algorithms over the 3 rd scenario. Since our method performs extra calculations such as combined mutation, one iteration takes 1.6418 s, excluding model runs. For MOGA, this time is 1.4321 s; for MOPSO, it is 1.2123 seconds. Although the single iteration times are shorter, the number of iterations required to find similar trajectories is higher in MOGA and MOPSO. Hence, the total times are higher for these methods. Essentially, for this comparison to be more meaningful, all algorithms should have been able to find similar route families.

VI. PERFORMANCE EVALUATIONS
In single-objective optimization, it is reasonably easy to measure the quality of a particular solution; if the problem under consideration is a minimization problem, the smaller the objective function value, the better. However, evaluating the quality of a Pareto set of solutions also requires examining several criteria, and 57 performance indicators are examined [42]. In our study, hypervolume [36], [43] is preferred to compare the proposed method with other methods. The hypervolume metric is sensitive to the scales of the objective functions; it is suggested to evaluate the [28] metric using normalized objective function values. Actually, in the Pareto Fronts, objective functions are evaluated as normalized values. Hypervolumes have been calculated using (Eq. 12) with normalized values [42].
where is the Pareto front approximation, is the reference point dominated by all points in , and is the Lebesgue measure. The gaps in the Pareto Fronts need special attention, as stated in [37]. We need to determine whether these gaps are due to the nature of the problem or an error in the algorithm. As mentioned earlier, there are two goals in the Multiobjective optimization problem: finding a well-converged and well-distributed set of Pareto-optimal solutions and choosing a single preferred solution from this set. Therefore, it must be determined that there is no solution (i.e., trajectory) associated with these gaps. A method for this issue is proposed in [37]. In the first stage of this method, it is investigated whether there is a gap or not. An Automated Gap finding procedure described in [37] has been used to find gap points ( j=1, 2, …,M). At the second stage, using the following two equations, we can find whether these gaps are formed due to the nature of the problem [37]. The method given in [40] has been used to find R near-PO point (i.e., ) The metric is the average normalized distance of gappoints from their closest points on Pareto Optimal front to the points on the Pareto Optimal. If this metric value is greater than one (i.e., > 1 ), the obtained gap-point can be confirmed to exist.
The metric is the average distance points to the gap-point compared to their distance from points on Pareto Optimal front. If > 1 there exists a natural gap in the Pareto-set, and it means it is impossible to find any solution in the gap region. If < 1 the observed gap is not natural; it is due to the problems in the algorithm. The metric is calculated as 1.151 and 1.782 for the first and the second scenarios, respectively.

VII. CONCLUSIONS
In this study, multi-objective trajectory planning has been implemented for a slung-load quadrotor system. The path length, scenario duration, and the dissipated energy are taken as objectives (costs). In the literature, multi-objective optimization problems are frequently solved by transforming them into a single-objective optimization problem for a weighted sum of costs (scalarization). However, determining proper weights for scalarization is an arduous task. Also, the scalarization approach prevents the investigation of the problem strategically from the perspective of objectives. In our practice, we can easily choose the proper weights by observing, especially the points corresponding to jumps on the Pareto Front (whenever hops exist).
We think that Pareto Front must be obtained before converting from a multi-objective optimization problem to a single objective problem in the type we have discussed here, where more than one objective is essential, and the importance of these objectives may vary according to the mission to be performed. As we emphasized in the results section, there may be points where we can't solve the problem by giving different weights to the single objective function. Therefore, before looking for a solution to the problem, we need to see a solution, if any.
Using the Pareto front, one may determine a point with short mission duration, sufficiently small energy expenditure, and not too long path length. Then, we can find the associated trajectory together with planned input. Finally, the physical vehicle is run with these inputs (Actually, assuming the vehicle has an autopilot, the obtained trajectory will be used as a guidance command for the autopilot).
One could add safety to these objective functions (can be considered the load's swing angle). However, we added this as a constraint in our problem. In the future, trajectory planning can be performed by evaluating the security objective over the scenarios that can pass over the obstacles or as a constraint in the vertical direction. Besides, we plan to extend the method and the solution algorithm to cases where more than one, two, or three quadrotors carry the load. Our algorithm will work very well as long as the obstacle locations are known. However, we plan to extend the present algorithm to cases with uncertainties in obstacle locations.