An Integer Programming Based Approach to Delivery Drone Routing under Load-Dependent Flight Speed

: Delivery drones have been attracting attention as a means of solving recent logistics issues, and many companies are focusing on their practical applications. Many research studies on delivery drones have been active for several decades. Among them, extended routing problems for drones have been proposed based on the Traveling Salesman Problem (TSP), which is used, for example, in truck vehicle routing problems. In parcel delivery by drones, additional constraints such as battery capacity, payload


Introduction
The number of parcel deliveries in logistics has been increasing every year.As a study of logistics, Lu et al. proposed a new optimization algorithm based on the combination of an ant colony system and gray wolf optimization to solve the fourth-party logistics routing problem [1].The coronavirus pandemic has further increased the number of parcel deliveries via e-commerce, especially small parcel deliveries.As a result, there are three major challenges in logistics: First, there is a labor shortage of drivers and other workers [2].This problem is called the 2024 problem [3].The second is traffic congestion caused by trucks.As the number of parcel delivery services increases, the number of trucks will inevitably increase.The third is the inefficiency of redelivery.Many people are unable to pick up their packages when they are delivered, so the delivery company must visit the customer again.As the number of parcels increases, the number of redelivery attempts also increases.
Drones 2023, 7, 320 2 of 14 Delivery drones are attracting attention as a promising solution to these problems.Compared to deliveries by truck, drones deliver packages while flying, thereby alleviating traffic congestion.In addition, drones emit less CO 2 since they are battery-powered.Furthermore, the unmanned flight of drones can reduce the labor costs required for delivery.Thus, delivery drones have many advantages.
However, because drones are battery-powered, they cannot fly for long periods of time.Therefore, to efficiently deliver packages by drone, the route must be determined in advance.There is a routing problem called TSP.TSP is the problem of finding the route with the shortest total travel distance among the routes on which a drone starts from the depot, visits every customer exactly once, and returns to the depot.Hoffman et al. defined TSP, formulated it, and presented a comprehensive overview of TSP, including the various mathematical models and algorithms developed to solve it [4].
In drone delivery, the drone's flight speed is greatly affected by the load; therefore, a delivery plan that takes changes in flight speed into account is necessary.Funabashi et al. proposed drone delivery planning (FSTSP), the authors have defined their problem as FSVRP, but it is actually FSTSP since they assumed a single drone routing).that takes into account changes in flight speed due to load and finds a route that minimizes the total flight time required for delivery [5].Their proposed FSTSP is not linearized and cannot be solved by a mathematical optimization solver because it contains trigonometric functions.In this paper, based on the fact that the problem can be easily extended, we propose integer quadratic programming and integer cubic programming for FSTSP to make it solvable by a general-purpose mathematical programming solver.In the experiments, the proposed method is compared with the conventional method in terms of solving time and total flight time.
The structure of this paper is as follows: Section 2 describes the literature review.Section 3 describes the Traveling Salesman Problem (TSP).Section 4 describes FSTSP.Section 5 describes an IQP approach and an ICP approach to the FSTSP.Section 6 describes the experimental method and the results of the evaluation experiments.Section 7 describes the overall summary and future work.

Literature Reviews
Research on package delivery has been active for several decades.There is research that focuses on the problem of truck-drone hybrid delivery routing, taking into consideration the dependency between the payload and the energy of the drone as well as the existence of No-Fly zones [6].This paper proposes a hybrid metaheuristic algorithm that aims to find the most efficient delivery route while considering these constraints.Other studies on algorithms include: Wen et al. proposed a new optimization algorithm based on colony search and global optimization to solve complex optimization problems [7].Lu et al. presented a bilevel optimization algorithm based on the whale optimization algorithm to schedule information technology projects considering outsourcing and risk management [8].Yan et al. developed a hybrid metaheuristic algorithm to solve the multi-objective locationrouting problem in the early post-disaster stage [9].
Research on drones for delivery has also been active.Drone delivery planning is one of the most important problems in realizing package delivery by drones [10,11].Masone et al. described a method for routing drones in a manner that involves multiple visits to certain locations while considering the launch and landing locations of the drone [12].The authors presented an iterative approach to solving this problem that combines both discrete and continuous improvements.We will describe drone delivery planning.
The Traveling Salesman Problem (TSP) has been proposed as a route planning problem.However, it is difficult to apply TSP to the drone route planning problem because additional constraints such as battery capacity, package weight, and weather conditions need to be considered for package delivery by drones.Drones in general are battery-powered, making it difficult to fly for long periods of time.Therefore, it is necessary to obtain efficient paths from aspects such as distance, energy, and time during routing.
Drones 2023, 7, 320 3 of 14 The drone needs to minimize its flight distance.In addition, the maximum payload capacity of drones is relatively small compared to that of trucks, trains, and other vehicles.Based on these considerations, Poikonen et al. proposed a path-planning problem to minimize the time required for drone delivery [13].In this problem, the flight speed is assumed to be constant in advance, and it is assumed that the flight speed does not change during delivery.
The drone needs to minimize energy consumption.To develop a more accurate energy consumption model, D'Andrea developed an energy consumption model for stable flight based on the lift-drag ratio [14].This model is characterized by the presence of two types: one with headwinds and the other without consideration of headwinds.Stolaroff et al. developed an energy consumption model that takes into account the forces to which the drone is subjected due to parasitic and induced drag [15], and Kirchstein developed an energy consumption model for takeoff, level flight, and delivery with landing [16].A feature of the model is the inclusion of downwash coefficients.Tseng developed the model based on horizontal and vertical speeds and accelerations, payload, mass, and wind speed data obtained from empirical experiments [17].This model is a function of package weight and airspeed.Carlos et al. used deep learning with Long Short-Term Memory (LSTM) to improve accuracy over previous methods [18].LSTM is a type of Recurrent Neural Network (RNN) that can learn long-term dependencies.
The drone needs to minimize flight time.In drone delivery, the drone's flight speed is greatly affected by the load; therefore, a delivery plan that takes changes in flight speed into account is necessary.There are also some studies that assume that the flight speed varies with the load.Fontaine presented a solution to a common logistics problem: optimizing the routes for delivery vehicles, specifically cargo bicycles [19].The author takes into account the fact that the travel time for a cargo bicycle may depend on the weight of the load it is carrying.Qian et al. focused on optimizing the routing of vehicles in order to minimize the fuel emissions produced by the vehicles [20].The problem being addressed takes into account not only the standard parameters in vehicle routing problems (such as the distance between stops, the capacity of the vehicles, etc.) but also the effect of varying speeds on the emissions produced by the vehicles.Rosati et al. focused on optimizing the routing of unmanned aerial vehicles (UAVs) in an ad-hoc network by considering the speed of the UAVs [21].Funabashi et al. proposed a drone delivery planning (FSTSP) system that takes into account changes in flight speed due to load and finds a route that minimizes the total flight time required for delivery [5].Their proposed FSTSP cannot be solved by a mathematical optimization solver because it contains trigonometric functions.In this paper, based on the fact that the problem can be easily extended, we propose integer quadratic programming and integer cubic programming for FSTSP in order to make it solvable by a general-purpose mathematical programming solver.In the experiments, the proposed method is compared with the conventional method in terms of solving time and total flight time.

Description of TSP
In recent years, delivery drones have been expected in many aspects, such as labor costs, traffic congestion, and environmental conservation, and many companies are aiming for practical use [22][23][24][25].To deliver packages to multiple customers within a set time using drones, it is necessary to determine delivery routes in advance.
There is a routing problem called TSP.TSP is the problem of finding the route with the smallest total travel distance among the routes in which a drone starts from the depot, visits every customer exactly once, and returns to the depot.Jünger et al. described the various mathematical models and algorithms developed to solve TSP, including exact algorithms, approximation algorithms, and heuristics [26].TSP is known as one of the NP-hard problems, and many heuristics, meta-heuristics, and exact algorithms have been developed over the decades.The simplest yet most efficient heuristic algorithm is the Nearest Neighborhood (NN) search algorithm.NN is an approach that selects the closest customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.
Drones 2023, 7, 0 4 of 13 customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.
Drones 2023, 7, x FOR PEER REVIEW 4 of 14 hard problems, and many heuristics, meta-heuristics, and exact algorithms have been developed over the decades.The simplest yet most efficient heuristic algorithm is the Nearest Neighborhood (NN) search algorithm.NN is an approach that selects the closest customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.Next, we describe the formulation of TSP.First, the  parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also .Packages are numbered from 1 to , and the customer to whom package  (1 ≤  ≤ ) is to be delivered is called customer .In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let (1, 2) be the distance between customers 1 and 2.Furthermore, () is the decision variable of the routing problem, which denotes the visiting order at the j-th The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.
Drones 2023, 7, x FOR PEER REVIEW 4 of 14 hard problems, and many heuristics, meta-heuristics, and exact algorithms have been developed over the decades.The simplest yet most efficient heuristic algorithm is the Nearest Neighborhood (NN) search algorithm.NN is an approach that selects the closest customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.Next, we describe the formulation of TSP.First, the  parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also .Packages are numbered from 1 to , and the customer to whom package  (1 ≤  ≤ ) is to be delivered is called customer .In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let (1, 2) be the distance between customers 1 and 2.Furthermore, () is the decision variable of the routing problem, which denotes the visiting order at the j-th Next, we describe the formulation of TSP.First, the N parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also N. Packages are numbered from 1 to N, and the customer to whom package i (1 ≤ i ≤ N) is to be delivered is called customer i.In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let d(i1, i2) be the distance between customers i1 and i2.Furthermore, x(j) is the decision variable of the routing problem, which denotes the visiting order at the j-th customer.For example, x(j) = 3 indicates that the drone visits the j-th customer at the third point.Since the start and end points of the route are the depot, this is defined as follows: The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.
Drones 2023, 7, 0 4 of 13 customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.
Drones 2023, 7, x FOR PEER REVIEW 4 of 14 hard problems, and many heuristics, meta-heuristics, and exact algorithms have been developed over the decades.The simplest yet most efficient heuristic algorithm is the Nearest Neighborhood (NN) search algorithm.NN is an approach that selects the closest customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.Next, we describe the formulation of TSP.First, the  parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also .Packages are numbered from 1 to , and the customer to whom package  (1 ≤  ≤ ) is to be delivered is called customer .In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let (1, 2) be the distance between customers 1 and 2.Furthermore, () is the decision variable of the routing problem, which denotes the visiting order at the j-th The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.
Drones 2023, 7, x FOR PEER REVIEW 4 of 14 hard problems, and many heuristics, meta-heuristics, and exact algorithms have been developed over the decades.The simplest yet most efficient heuristic algorithm is the Nearest Neighborhood (NN) search algorithm.NN is an approach that selects the closest customers one by one until all customers have been visited.Algorithms for solving TSP include dynamic programming developed by Bellman [27] and Held et al. [28].
Here, Figure 1 shows an example of the problem.A node written as "0" indicates the depot, and the other three nodes indicate customers.The number inside the box indicates the weight of the package requested by the customer, and the number on the edge connecting each node indicates the distance between two points.Since the drone moves by flying, it is assumed that there are no geographical restrictions and that it can always move between arbitrary customers.The optimal route of TSP in this problem is shown in Figure 2a,b.The route indicated by the arrow is the shortest route.In this route, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 10 at customer 1, and travels to customer 2 with a load of 80 (= 60 + 20).After that, the drone unloads a load of 60 at customer 2 and goes to customer 3 with a load of 20.The drone unloads a load of 20 at customer 3 and returns to the depot with no load.The total distance with TSP is 164 (= 28 + 42 + 54 + 40), and the total flight time is 36 (= 7 + 10 + 11 + 8).If the flight speed is constant, it can be said that the shortest route is optimal in terms of both flight distance and flight time.Next, we describe the formulation of TSP.First, the  parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also .Packages are numbered from 1 to , and the customer to whom package  (1 ≤  ≤ ) is to be delivered is called customer .In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let (1, 2) be the distance between customers 1 and 2.Furthermore, () is the decision variable of the routing problem, which denotes the visiting order at the j-th Next, we describe the formulation of TSP.First, the N parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also N. Packages are numbered from 1 to N, and the customer to whom package i (1 ≤ i ≤ N) is to be delivered is called customer i.In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let d(i1, i2) be the distance between customers i1 and i2.Furthermore, x(j) is the decision variable of the routing problem, which denotes the visiting order at the j-th customer.For example, x(j) = 3 indicates that the drone visits the j-th customer at the third point.Since the start and end points of the route are the depot, this is defined as follows: Next, we describe the formulation of TSP.First, the N parcels to be delivered are given.Without loss of generality, no more than one package is delivered to the same customer.Therefore, the number of customers is also N. Packages are numbered from 1 to N, and the customer to whom package i (1 ≤ i ≤ N) is to be delivered is called customer i.In this paper, we assume that the drone can deliver all packages in one flight.When the drone leaves the depot, it loads all the packages and departs.If the total weight of the package exceeds the maximum payload of the drone, it is necessary to divide the package into groups before delivery planning, but the method of dividing is not covered in this paper.
Let d(i1, i2) be the distance between customers i1 and i2.Furthermore, x(j) is the decision variable of the routing problem, which denotes the visiting order at the j-th customer.For example, x(j) = 3 indicates that the drone visits the j-th customer at the third point.Since the start and end points of the route are the depot, this is defined as follows: Drones 2023, 7, 320 5 of 14 Since every customer must be visited only once, it is defined as follows: Since general TSP seeks the route with the shortest total flight distance, the objective function is defined as follows:

Description of FSTSP
As a derivative of TSP described in Section 2, there is a problem called the Flight Speed-Aware Traveling Salesman Problem (FSTSP) [5].In this problem, given a set of packages to be delivered, the optimal flight path is found to depart from a delivery location, deliver all packages to customers, and return to the delivery location.Here, we focus on the problem of finding the flight path with the shortest flight time, taking into account changes in flight speed due to load.

An Example
Figure 3a,b show the optimal paths of the FSTSP for the problem in Figure 1.The path indicated by the arrows is the shortest path in time.In this path, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 60 at customer 2, and travels to customer 3 with a load of 30 (= 10 + 20).It then unloads a load of 20 at customer 3, travels to customer 1 with a load of 10, unloads a load of 10 at customer 1, and returns to the depot with a load of 0. At this time, the total distance flown by FSTSP is 168 (= 22 + 54 + 64 + 28).However, the total flight time for FSTSP is 35 (= 5 + 11 + 13 + 6), which is shorter than the 36 (= 7 + 10 + 11 + 8) total flight time for the TSP optimal path.In general, the heavier the load carried, the slower the drone flies, and the lighter the load, the faster the drone flies.In Figure 2a, the drone first visits customer 1 and unloads a load of 10.The drone then continues flying with a load of 80 (= 60 + 20).On the other hand, in Figure 3a, the drone first visits customer 2 and unloads the heaviest load, 60.The drone then continues flying with a load of 30 (= 10 + 20).This allows the drone to fly faster than the TSP during subsequent deliveries because of the lighter weight of the onboard luggage.This example shows that the optimal path for TSP is not always the same as the optimal path for FSTSP.
Drones 2023, 7, 0 5 of 13 Since every customer must be visited only once, it is defined as follows: Since general TSP seeks the route with the shortest total flight distance, the objective function is defined as follows:

Description of FSTSP
As a derivative of TSP described in Section 2, there is a problem called the Flight Speed-Aware Traveling Salesman Problem (FSTSP) [5].In this problem, given a set of packages to be delivered, the optimal flight path is found to depart from a delivery location, deliver all packages to customers, and return to the delivery location.Here, we focus on the problem of finding the flight path with the shortest flight time, taking into account changes in flight speed due to load.

An Example
Figure 3a,b show the optimal paths of the FSTSP for the problem in Figure 1.The path indicated by the arrows is the shortest path in time.In this path, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 60 at customer 2, and travels to customer 3 with a load of 30 (= 10 + 20).It then unloads a load of 20 at customer 3, travels to customer 1 with a load of 10, unloads a load of 10 at customer 1, and returns to the depot with a load of 0. At this time, the total distance flown by FSTSP is 168 (= 22 + 54 + 64 + 28).However, the total flight time for FSTSP is 35 (= 5 + 11 + 13 + 6), which is shorter than the 36 (= 7 + 10 + 11 + 8) total flight time for the TSP optimal path.In general, the heavier the load carried, the slower the drone flies, and the lighter the load, the faster the drone flies.In Figure 2a, the drone first visits customer 1 and unloads a load of 10.The drone then continues flying with a load of 80 (= 60 + 20).On the other hand, in Figure 3a, the drone first visits customer 2 and unloads the heaviest load, 60.The drone then continues flying with a load of 30 (= 10 + 20).This allows the drone to fly faster than the TSP during subsequent deliveries because of the lighter weight of the onboard luggage.This example shows that the optimal path for TSP is not always the same as the optimal path for FSTSP.
customer.For example, x(j) = 3 indicates that the drone visits the j-th customer at the third point.Since the start and end points of the route are the depot, this is defined as follows: Since every customer must be visited only once, it is defined as follows: Since general TSP seeks the route with the shortest total flight distance, the objective function is defined as follows: ∑ ((), ( + 1))  =0 (4)

Description of FSTSP
As a derivative of TSP described in Section 2, there is a problem called the Flight Speed-Aware Traveling Salesman Problem (FSTSP) [5].In this problem, given a set of packages to be delivered, the optimal flight path is found to depart from a delivery location, deliver all packages to customers, and return to the delivery location.Here, we focus on the problem of finding the flight path with the shortest flight time, taking into account changes in flight speed due to load.

An Example
Figure 3a,b show the optimal paths of the FSTSP for the problem in Figure 1.The path indicated by the arrows is the shortest path in time.In this path, the drone departs from the depot with a load of 90 (= 10 + 60 + 20), unloads a load of 60 at customer 2, and travels to customer 3 with a load of 30 (= 10 + 20).It then unloads a load of 20 at customer 3, travels to customer 1 with a load of 10, unloads a load of 10 at customer 1, and returns to the depot with a load of 0. At this time, the total distance flown by FSTSP is 168 (= 22 + 54 + 64 + 28).However, the total flight time for FSTSP is 35 (= 5 + 11 + 13 + 6), which is shorter than the 36 (= 7 + 10 + 11 + 8) total flight time for the TSP optimal path.In general, the heavier the load carried, the slower the drone flies, and the lighter the load, the faster the drone flies.In Figure 2a, the drone first visits customer 1 and unloads a load of 10.The drone then continues flying with a load of 80 (= 60 + 20).On the other hand, in Figure 3a, the drone first visits customer 2 and unloads the heaviest load, 60.The drone then continues flying with a load of 30 (= 10 + 20).This allows the drone to fly faster than the TSP during subsequent deliveries because of the lighter weight of the onboard luggage.This example shows that the optimal path for TSP is not always the same as the optimal path for FSTSP.

FSTSP Formulation
Let w(i) be the weight of package i and W(j) be the total weight when the drone leaves the j-th customer.When the drone begins its delivery, all packages are loaded, therefore the following Equation holds.

FSTSP Formulation
Let w(i) be the weight of package i and W(j) be the total weight when the drone leaves the j-th customer.When the drone begins its delivery, all packages are loaded, therefore the following Equation holds.
Drones 2023, 7, 320 6 of 14 When the drone visits the j-th (1 ≤ j ≤ N) point x(j), it unloads a load w(x(j)).Thus, the total weight when it finishes visiting the point x(j) is defined as Let t(i1, i2) be the flight time between points i1 and i2.The time taken by the drone to fly from i1 to i2 is the distance between i1 and i2 divided by the drone's flight speed v.The drone's flight speed is a function of the load variable, and the flight time between the j-th and (j + 1)-th drone is defined as FSTSP is the problem of finding the path with the shortest total flight time.Therefore, the objective function is defined as Here is a brief approximation of the effect of the weight of the parcel carried by the drone on the flight speed v.The function of flight speed, v, within the problem definition depends on the drone.v is assumed to be given, and the exact v is beyond the scope of this paper.
Figure 4a shows a drone flying horizontally without a parcel, where P represents the lift of the drone and W d represents the weight of the drone itself.The horizontal component P x of P is the force toward the destination and is equal to the aerodynamic drag k•v(0).Therefore, the following Equation holds, where k is a coefficient specific to the drone.
To prevent the drone from falling, the vertical component of P, P y , must be equal to the gravity, W d •g.When θ is the pitch angle and P y is equal to the gravity, the following Equation holds.
Figure 4b shows the drone as it carries the load w.To avoid crashing, the pitch angle θ must be less than θ.Therefore, the following Equations hold for the lift forces P x and P y of the drone.
Therefore, the drone's flight speed v(w), pitch angle θ, and θ can be defined by the following Equations, respectively.

Approximation Approaches by Integer Programming
This section describes how to solve the FSTSP described in Section 3 using a generalpurpose mathematical optimization solver.In this paper, we ascribe FSTSP to integer quadratic programming and integer cubic programming problems.The significance of using a mathematical optimization solver to solve the problem is that it facilitates the extension of the problem.When designing algorithms such as dynamic programming, it is necessary to redesign the algorithm from scratch when expanding the problem.However, by using mathematical optimization solvers, even if the problem needs to be expanded, it can be easily extended by adding only one or two lines of code.Nuzhet et al. proposed a mixed-integer linear programming approach to address the problem of task allocation for multiple heterogeneous robots in an unknown environment, subject to predefined constraints [29].The authors introduced several modifications to both the objective function and constraints to illustrate the flexibility of the mixed-integer linear programming formulation.

Approximation Approach
From Equation ( 8), if 1/ can be expressed by a quadratic or cubic expression, FSTSP becomes an integer quadratic programming problem or an integer cubic programming problem.In other words, we want to approximate Equation ( 17) with a linear or quadratic function.

1/𝑣(𝑤) = 1/ (𝑘 • 𝑃 • sin (arccos ( (𝑊
Different drones have different values of , , and   . is a constant.  is obvious. and  need to be determined.We will explain how to find them. The range of w also varies from drone to drone.Within the range of w for each individual drone, we can find an approximation of Equation (17).For a linear approximation, define Equation ( 18) and use the least-squares method to obtain the values of  and .
For quadratic approximation,

Approximation Approaches by Integer Programming
This section describes how to solve the FSTSP described in Section 3 using a generalpurpose mathematical optimization solver.In this paper, we ascribe FSTSP to integer quadratic programming and integer cubic programming problems.The significance of using a mathematical optimization solver to solve the problem is that it facilitates the extension of the problem.When designing algorithms such as dynamic programming, it is necessary to redesign the algorithm from scratch when expanding the problem.However, by using mathematical optimization solvers, even if the problem needs to be expanded, it can be easily extended by adding only one or two lines of code.Nuzhet et al. proposed a mixedinteger linear programming approach to address the problem of task allocation for multiple heterogeneous robots in an unknown environment, subject to predefined constraints [29].The authors introduced several modifications to both the objective function and constraints to illustrate the flexibility of the mixed-integer linear programming formulation.

Approximation Approach
From Equation ( 8), if 1/v can be expressed by a quadratic or cubic expression, FSTSP becomes an integer quadratic programming problem or an integer cubic programming problem.In other words, we want to approximate Equation ( 17) with a linear or quadratic function.
Different drones have different values of k, P, and W d .g is a constant.W d is obvious.k and P need to be determined.We will explain how to find them.
The range of w also varies from drone to drone.Within the range of w for each individual drone, we can find an approximation of Equation (17).For a linear approximation, define Equation ( 18) and use the least-squares method to obtain the values of a and b.For quadratic approximation, define Equation (19) and use the least-squares method to obtain the values of a, b and c.

Case of AR Drone 2.0
This section describes the flight speed approximation for AR Drone 2.0 [30].First, how to obtain k and P will be explained.P can be obtained from the force balance when the maximum load capacity is loaded.

Approximation Approaches by Integer Programming
This section describes how to solve the FSTSP described in Section 3 using a generalpurpose mathematical optimization solver.In this paper, we ascribe FSTSP to integer quadratic programming and integer cubic programming problems.The significance of using a mathematical optimization solver to solve the problem is that it facilitates the extension of the problem.When designing algorithms such as dynamic programming, it is necessary to redesign the algorithm from scratch when expanding the problem.However, by using mathematical optimization solvers, even if the problem needs to be expanded, it can be easily extended by adding only one or two lines of code.Nuzhet et al. proposed a mixedinteger linear programming approach to address the problem of task allocation for multiple heterogeneous robots in an unknown environment, subject to predefined constraints [29].The authors introduced several modifications to both the objective function and constraints to illustrate the flexibility of the mixed-integer linear programming formulation.

Approximation Approach
From Equation (8), if 1/v can be expressed by a quadratic or cubic expression, FSTSP becomes an integer quadratic programming problem or an integer cubic programming problem.In other words, we want to approximate Equation ( 17) with a linear or quadratic function.
Different drones have different values of k, P, and W d .g is a constant.W d is obvious.k and P need to be determined.We will explain how to find them.
The range of w also varies from drone to drone.Within the range of w for each individual drone, we can find an approximation of Equation (17).For a linear approximation, define Equation ( 18) and use the least-squares method to obtain the values of a and b.For quadratic approximation, define Equation ( 19) and use the least-squares method to obtain the values of a, b and c.

Case of AR Drone 2.0
This section describes the flight speed approximation for AR Drone 2.0 [30].First, how to obtain k and P will be explained.P can be obtained from the force balance when the maximum load capacity is loaded.
Drones 2023, 7, 320 8 of 14 Here, W d is 490 g and w max is 250 g.The acceleration of gravity g is assumed to be 9.8 m/s 2 .Substituting these two values into the Equation ( 20), P is 7252 N. Furthermore, k is the value obtained by dividing the velocity when the load is 0 by P sin(θ).
v(0) is 5 m/s.Furthermore, θ was derived using Equation (14).Substituting the values of v(0), P, and θ into Equation ( 21), k is 0.00092.Regarding the weight w of the luggage, AR Drone 2.0 is within the range of 0 w 200.
The black curve in Figure 5 represents the reciprocal of Equation ( 16).The vertical axis of this figure represents the flight speed, and the horizontal axis represents the weight of the payload.The green line in the figure is linearly approximated using the method of least squares.The red curve in the figure is quadratically approximated using the method of least squares.The R-squared value indicating the error from Equation ( 22) by approximation was 0.9019, and the R-squared value indicating the error from Equation ( 23) was 0.9895.
Finally, we redefine the objective function of Equation ( 8) in Section 2 by transforming it into a quadratic or cubic approximation and solving it with a mathematical optimization solver.
Drones 2023, 7, 0 8 of 13 Here, W d is 490 g and w max is 250 g.The acceleration of gravity g is assumed to be 9.8 m/s 2 .Substituting these two values into the Equation ( 20), P is 7252 N. Furthermore, k is the value obtained by dividing the velocity when the load is 0 by P sin(θ).
v(0) is 5 m/s.Furthermore, θ was derived using Equation (14).Substituting the values of v(0), P, and θ into Equation ( 21), k is 0.00092.Regarding the weight w of the luggage, AR Drone 2.0 is within the range of 0 w 200.
The black curve in Figure 5 represents the reciprocal of Equation ( 16).The vertical axis of this figure represents the flight speed, and the horizontal axis represents the weight of the payload.The green line in the figure is linearly approximated using the method of least squares.The red curve in the figure is quadratically approximated using the method of least squares.The R-squared value indicating the error from Equation ( 22) by approximation was 0.9019, and the R-squared value indicating the error from Equation ( 23) was 0.9895.
Finally, we redefine the objective function of Equation ( 8) in Section 2 by transforming it into a quadratic or cubic approximation and solving it with a mathematical optimization solver.

Case of SkyLift
This section describes the flight speed approximation for SkyLift [31].First, how to obtain  and  will be explained. can be obtained from the force balance when the maximum load capacity is loaded.Here,   is 55,000 g and   is 30,000 g.Substituting these two values into Equation ( 20),  is 833,000 N. Furthermore,  is the value obtained by dividing the velocity when the load is 0 by  sin().(0) is 10. was derived using Equation (14).Substituting the values of (0), P, and  into the Equation ( 21

Load (g)
The black curve is the reciprocal of Equation ( 16) The red curve The green line

Case of SkyLift
This section describes the flight speed approximation for SkyLift [31].First, how to obtain k and P will be explained.P can be obtained from the force balance when the maximum load capacity is loaded.Here, W d is 55,000 g and w max is 30,000 g.Substituting these two values into Equation ( 20), P is 833,000 N. Furthermore, k is the value obtained

Case SkyLift
This section describes the flight speed approximation for SkyLift [31].First, how to obtain k and P will be explained.P can be obtained from the force balance when the maximum load capacity is loaded.Here, W d is 55,000 g and w max is 30,000 g.Substituting Drones 2023, 7, 320 9 of 14 these two values into Equation ( 20), P is 833,000 N. Furthermore, k is the value obtained by dividing the velocity when the load is 0 by P sin(θ).v(0) is 10.θ was derived using Equation (14).Substituting the values of v(0), P, and θ into the Equation ( 21), k is 0.000015.Regarding the weight w of the luggage, SkyLift is within the range of 0 w 27, 000.
The black curve in Figure 6 represents the reciprocal of Equation ( 16).The vertical axis of this figure represents the flight speed, and the horizontal axis represents the weight of the payload.The green line in the figure is linearly approximated using the method of least squares.The red curve in the figure is quadratically approximated using the method of least squares.The R-squared value indicating the error from Equation ( 26) by approximation was 0.8214, and the R-squared value indicating the error from Equation ( 27) was 0.9647.
As described in Section 5.2, we redefine the objective function of Equation ( 8) by transforming it into a quadratic or cubic approximation and solving it with a mathematical optimization solver.
Using these Equations ( 1)-( 6) and Equation (28) or Equation (29), it is possible to solve the quadratic or cubic problem FSTSP with a general-purpose mathematical optimization solver.
Drones 2023, 7, 0 9 of 13 by dividing the velocity when the load is 0 by P sin(θ).v(0) is 10.θ was derived using Equation (14).Substituting the values of v(0), P, and θ into the Equation ( 21), k is 0.000015.Regarding the weight w of the luggage, SkyLift is within the range of 0 w 27, 000.The black curve in Figure 6 represents the reciprocal of Equation ( 16).The vertical axis of this figure represents the flight speed, and the horizontal axis represents the weight of the payload.The green line in the figure is linearly approximated using the method of least squares.The red curve in the figure is quadratically approximated using the method of least squares.The R-squared value indicating the error from Equation ( 26) by approximation was 0.8214, and the R-squared value indicating the error from Equation ( 27) was 0.9647.
As described in Section 5.2, we redefine the objective function of Equation ( 8) by transforming it into a quadratic or cubic approximation and solving it with a mathematical optimization solver.
Using these Equations ( 1)-( 6) and Equation (28) or Equation (29), it is possible to solve the quadratic or cubic problem FSTSP with a general-purpose mathematical optimization solver.

Evaluation
In this section, we demonstrate the effectiveness of the proposed method through experiments.In this paper, the proposed method and the algorithms to be compared are implemented by Python and CPLEX, and the solution time and total flight time of each algorithm are compared.The experiments are conducted on AMD Ryzen 7 PRO 4750G (8 cores, 16 threads) with 64 GB of main memory.For the solution calculation, the existing method uses Python 3.8.5, and the proposed method uses IBM ILOG CPLEX Optimization Studio 20.1.

Experimental Setups
This section describes the six methods used in the experiments.Each method is de- The black curve is the reciprocal of Equation ( 16) The red curve The green line

Evaluation
In this section, we demonstrate the effectiveness of the proposed method through experiments.In this paper, the proposed method and the algorithms to be compared are implemented by Python and CPLEX, and the solution time and total flight time of each algorithm are compared.The experiments are conducted on AMD Ryzen 7 PRO 4750G (8 cores, 16 threads) with 64 GB of main memory.For the solution calculation, the existing method uses Python 3.8.5, and the proposed method uses IBM ILOG CPLEX Optimization Studio 20.1.

6.
In this section, we demonstrate the effectiveness of the proposed method through experiments.In this paper, the proposed method and the algorithms to be compared are implemented by Python and CPLEX, and the solution time and total flight time of each algorithm are compared.The experiments are conducted on AMD Ryzen 7 PRO 4750G (8 cores, 16 threads) with 64 GB of main memory.For the solution calculation, the existing method uses Python 3.8.5, and the proposed method uses IBM ILOG CPLEX Optimization Studio 20.1.

Experimental Setups
This section describes the six methods used in the experiments.Each method is described below.

•
DP-TSP: This method solves TSP using dynamic programming.

•
CPLEX-IQP-single: This method solves FSTSP based on an IQP approach by quadratic approximation of the objective with CPLEX on a single thread.

•
CPLEX-IQP-multi: This method solves FSTSP based on an IQP approach by quadratic approximation of the objective with CPLEX on multiple threads.

•
CPLEX-ICP-single: This method solves FSTSP based on an ICP approach by cubic approximation of the objective with CPLEX on a single thread.

•
CPLEX-ICP-multi: This method solves FSTSP based on an ICP approach by cubic approximation of the objective with CPLEX on multiple threads.
In the experiments, we assumed two types of drones (e.g., AR Drone 2.0 and Skylift), which are very different in size and maximum payload.Since two drones, AR Drone 2.0 and SkyLift, were used in this experiment, the name of the drone is given at the beginning of the name of each proposed method, for example, ARDrone-CPLEX-IQP-single.
The runtime is limited to 3600 s in wall-clock time.In addition, benchmarks with five to twenty customers were prepared.A total of 320 problems were prepared, 20 for each number of customers, and the average of these was used as the experimental result.The coordinates of each destination and the weight of the package requested by each destination were set randomly.The total weight of the package was randomly set within the range not exceeding the maximum payload capacity for each drone.For example, if 20 packages are to be delivered, the weight of each package is randomly assigned so that the total weight of the 20 packages does not exceed the maximum load capacity.

Experimental Results
Figure 7 shows a comparison of the runtimes of the six algorithms of AR Drone 2.0 using a logarithmic function.The vertical axis is runtime, and the horizontal axis is the number of customers.In this experiment, DP-TSP is not used for comparison because we focus on the solution time of the method that solves FSTSP.From the experimental results with five to twelve customers, DP-FSTSP can find solutions faster than other methods.This is because CPLEX incurs overheads such as problem generation and parallelization across multiple threads.In addition, the proposed methods, CPLEX-IQP-single and CPLEX-IQPmulti, have shorter solution times than DP-FSTSP in cases where the number of delivery customers is 14 or more.Since CPLEX-IQP-multi is executed with multiple threads, the solution is shorter than CPLEX-IQP-single, which is executed with a single thread when the number of customers is 14 or more.
Figure 8 shows the total flight time results normalized by DP-FSTSP when using AR Drone 2.0.The vertical axis is the result of total flight time normalized by DP-FSTSP, and the horizontal axis is the number of customers.In the proposed method, whether the problem is solved using a single thread or multiple threads does not affect the total flight time results.ARDrone-CPLEX-IQP-single and ARDrone-CPLEX-IQP-multi are integrated into ARDrone-CPLEX-IQP, and ARDrone-CPLEX-ICP-single and ARDrone-CPLEX-ICP-multi are integrated into ARDrone-CPLEX-ICP.The experimental results show that DP-TSP has a longer total flight time than other methods.This indicates that DP-TSP deals with the problem of optimizing the flight distance, so it is difficult to obtain the optimal solution in terms of flight time.result shows that the minimization of the flight distance is not the minimization of the flight time.In addition, the solutions of the four proposed algorithms, ARDrone-CPLEX-IQP was only 0.28% worse on average than the solutions of DP-FSTSP.ARDrone-CPLEX-ICP was only 1.26% worse on average than the solutions of DP-FSTSP.
From this result, it can be said that the degradation of the solution is hardly seen even if the problem is approximated.
focus on the solution time of the method that solves FSTSP.From the experime with five to twelve customers, DP-FSTSP can find solutions faster than other me is because CPLEX incurs overheads such as problem generation and paralleliza multiple threads.In addition, the proposed methods, CPLEX-IQP-single and C multi, have shorter solution times than DP-FSTSP in cases where the number customers is 14 or more.Since CPLEX-IQP-multi is executed with multiple th solution is shorter than CPLEX-IQP-single, which is executed with a single th the number of customers is 14 or more.
with five to twelve customers, DP-FSTSP can find solutions faster than other methods.This is because CPLEX incurs overheads such as problem generation and parallelization across multiple threads.In addition, the proposed methods, CPLEX-IQP-single and CPLEX-IQP-multi, have shorter solution times than DP-FSTSP in cases where the number of delivery customers is 14 or more.Since CPLEX-IQP-multi is executed with multiple threads, the solution is shorter than CPLEX-IQP-single, which is executed with a single thread when the number of customers is 14 or more.Figure 8 shows the total flight time results normalized by DP-FSTSP when using AR Drone 2.0.The vertical axis is the result of total flight time normalized by DP-FSTSP, and the horizontal axis is the number of customers.In the proposed method, whether the problem is solved using a single thread or multiple threads does not affect the total flight time results.ARDrone-CPLEX-IQP-single and ARDrone-CPLEX-IQP-multi are integrated into ARDrone-CPLEX-IQP, and ARDrone-CPLEX-ICP-single and ARDrone-CPLEX-ICPmulti are integrated into ARDrone-CPLEX-ICP.The experimental results show that DP-TSP has a longer total flight time than other methods.This indicates that DP-TSP deals with the problem of optimizing the flight distance, so it is difficult to obtain the optimal solution in terms of flight time.This result shows that the minimization of the flight distance is not the minimization of the flight time.In addition, the solutions of the four proposed algorithms, ARDrone-CPLEX-IQP was only 0.28% worse on average than the solutions of DP-FSTSP.ARDrone-CPLEX-ICP was only 1.26% worse on average than the solutions of DP-FSTSP.From this result, it can be said that the degradation of the solution is hardly seen even if the problem is approximated.Drones 2023, 7, 0 Figure 8 shows the total flight time results normalized by DP-FSTSP when u Drone 2.0.The vertical axis is the result of total flight time normalized by DP-FSTSP horizontal axis is the number of customers.In the proposed method, whether the is solved using a single thread or multiple threads does not affect the total flig results.ARDrone-CPLEX-IQP-single and ARDrone-CPLEX-IQP-multi are integra ARDrone-CPLEX-IQP, and ARDrone-CPLEX-ICP-single and ARDrone-CPLEX-IC are integrated into ARDrone-CPLEX-ICP.The experimental results show that DPa longer total flight time than other methods.This indicates that DP-TSP deals problem of optimizing the flight distance, so it is difficult to obtain the optimal so terms of flight time.This result shows that the minimization of the flight distance i minimization of the flight time.In addition, the solutions of the four proposed alg ARDrone-CPLEX-IQP was only 0.28% worse on average than the solutions of D ARDrone-CPLEX-ICP was only 1.26% worse on average than the solutions of DP From this result, it can be said that the degradation of the solution is hardly seen the problem is approximated.Figure 9 shows a comparison of the runtimes of the six algorithms of SkyLift using a logarithmic function.The vertical axis is runtime, and the horizontal axis is the number of customers.In the case of the number of customer(s) 5, the method using a single thread can find the solution faster than the method using multiple threads.The reason is that CPLEX incurs overhead due to parallelization across multiple threads.Due to the multithreaded execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a shorter runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These results show that quadratic approximations can be solved faster than cubic approximations.Since the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of customer(s) 9, we do not consider this to be a particular problem in practical use.  Figure 9 shows a comparison of the runtimes of the six algorithms of SkyLift logarithmic function.The vertical axis is runtime, and the horizontal axis is the nu customers.In the case of the number of customer(s) 5, the method using a single th find the solution faster than the method using multiple threads.The reason is tha incurs overhead due to parallelization across multiple threads.Due to the multi-t execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These show that quadratic approximations can be solved faster than cubic approximation the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of custom we do not consider this to be a particular problem in practical use. Figure 9 shows a comparison of the runtimes of the six algorithms of SkyLift using a logarithmic function.The vertical axis is runtime, and the horizontal axis is the number of customers.In the case of the number of customer(s) 5, the method using a single thread can find the solution faster than the method using multiple threads.The reason is that CPLEX incurs overhead due to parallelization across multiple threads.Due to the multithreaded execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a shorter runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These results show that quadratic approximations can be solved faster than cubic approximations.Since the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of customer(s) 9, we do not consider this to be a particular problem in practical use.  Figure 9 shows a comparison of the runtimes of the six algorithms of SkyLift using a logarithmic function.The vertical axis is runtime, and the horizontal axis is the number of customers.In the case of the number of customer(s) 5, the method using a single thread can find the solution faster than the method using multiple threads.The reason is that CPLEX incurs overhead due to parallelization across multiple threads.Due to the multi-threaded execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a shorter runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These results show that quadratic approximations can be solved faster than cubic approximations.Since the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of customer(s) 9, we do not consider this to be a particular problem in practical use.
Figure 10 shows the total flight time results normalized by DP-FSTSP when using SkyLift.The vertical axis is the result of total flight time normalized by DP-FSTSP, and the horizontal axis is the number of customers.In the proposed method, whether the problem is solved using a single thread or multiple threads does not affect the results.SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-IQP-multi are integrated into SkyLift-CPLEX-IQP, and SkyLift-CPLEX-ICP-single and SkyLift-CPLEX-ICP-multi are integrated into SkyLift-CPLEX-ICP.SkyLift-CPLEX-IQP was only 0.53% worse on average than the solutions by DP-FSTSP.SkyLift-CPLEX-ICP was only 2.54% worse on average than the solutions by DP-FSTSP.
customers.In the case of the number of customer(s) 5, the method using a single thread can find the solution faster than the method using multiple threads.The reason is that CPLEX incurs overhead due to parallelization across multiple threads.Due to the multithreaded execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a shorter runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These results show that quadratic approximations can be solved faster than cubic approximations.Since the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of customer(s) 9, we do not consider this to be a particular problem in practical use. Figure 10 shows the total flight time results normalized by DP-FSTSP when using SkyLift.The vertical axis is the result of total flight time normalized by DP-FSTSP, and the horizontal axis is the number of customers.In the proposed method, whether the problem is solved using a single thread or multiple threads does not affect the results.SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-IQP-multi are integrated into SkyLift-CPLEX-IQP, and SkyLift-CPLEX-ICP-single and SkyLift-CPLEX-ICP-multi are integrated into SkyLift-CPLEX-ICP.SkyLift-CPLEX-IQP was only 0.53% worse on average than the solutions by customers.In the case of the number of customer(s) 5, the method using a single th find the solution faster than the method using multiple threads.The reason is tha incurs overhead due to parallelization across multiple threads.Due to the multi-t execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.Thes show that quadratic approximations can be solved faster than cubic approximatio the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of custo we do not consider this to be a particular problem in practical use.
customers.In the case of the number of customer(s) 5, the method using a single thread can find the solution faster than the method using multiple threads.The reason is that CPLEX incurs overhead due to parallelization across multiple threads.Due to the multithreaded execution, SkyLift-CPLEX-IQP-multi and SkyLift-CPLEX-ICP-multi can achieve a shorter runtime than SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-ICP-single.These results show that quadratic approximations can be solved faster than cubic approximations.Since the runtime of SkyLift-CPLEX-ICP is about 1 min for up to the number of customer(s) 9, we do not consider this to be a particular problem in practical use. Figure 10 shows the total flight time results normalized by DP-FSTSP when using SkyLift.The vertical axis is the result of total flight time normalized by DP-FSTSP, and the horizontal axis is the number of customers.In the proposed method, whether the problem is solved using a single thread or multiple threads does not affect the results.SkyLift-CPLEX-IQP-single and SkyLift-CPLEX-IQP-multi are integrated into SkyLift-CPLEX-IQP, and SkyLift-CPLEX-ICP-multi are integrated into SkyLift-CPLEX-ICP.SkyLift-CPLEX-IQP was only 0.53% worse on average than the solutions by  Drones 2023, 7, 0 12 Figure 10 shows the total flight time results normalized by DP-FSTSP when u SkyLift.The vertical axis is the result of total flight time normalized by DP-FSTSP, and horizontal axis is the number of customers.In the proposed method, whether the prob is solved using a single thread or multiple threads does not affect the results.Sky CPLEX-IQP-single and SkyLift-CPLEX-IQP-multi are integrated into SkyLift-CPLEX and SkyLift-CPLEX-ICP-single and SkyLift-CPLEX-ICP-multi are integrated into Sky CPLEX-ICP.SkyLift-CPLEX-IQP was only 0.53% worse on average than the solut by DP-FSTSP.SkyLift-CPLEX-ICP was only 2.54% worse on average than the solut by DP-FSTSP.

Conclusions
In this work, we have presented IQP-based and ICP-based approaches for delivering drone routing under load-dependent flight speeds.Experimental results show that the IQP approach is superior to the ICP approach in terms of the runtime and total flight time as the number of customers increases.However, if the number of customers is small, the ICP approach can obtain a better solution within one minute than the IQP approach.Therefore, we believe that the ICP approach is also valid in practical use.
In future work, we will consider multiple deliveries to routing problems of a drone considering load-dependent flight speeds.
In the experiments in this paper, we have limited the number of customers to a maximum of 20.In the future, when visiting more than 20 customers in a single trip is possible, we plan to develop more efficient approaches to delivering drone routing.

Conclusions
In this work, we have presented IQP-based and ICP-based approaches for delive drone routing under load-dependent flight speeds.Experimental results show that the approach is superior to the ICP approach in terms of the runtime and total flight tim the number of customers increases.However, if the number of customers is small, the approach can obtain a better solution within one minute than the IQP approach.There we believe that the ICP approach is also valid in practical use.
In future work, we will consider multiple deliveries to routing problems of a d considering load-dependent flight speeds.
In the experiments in this paper, we have limited the number of customers maximum of 20.In the future, when visiting more than 20 customers in a single tr possible, we plan to develop more efficient approaches to delivering drone routing.

Conclusions
In this work, we have presented IQP-based and ICP-based approaches for delivering drone routing under load-dependent flight speeds.Experimental results show that the IQP approach is superior to the ICP approach in terms of the runtime and total flight time as the number of customers increases.However, if the number of customers is small, the ICP approach can obtain a better solution within one minute than the IQP approach.Therefore, we believe that the ICP approach is also valid in practical use.
In future work, we will consider multiple deliveries to routing problems of a drone considering load-dependent flight speeds.
In the experiments in this paper, we have limited the number of customers to a maximum of 20.In the future, when visiting more than 20 customers in a single trip is possible, we plan to develop more efficient approaches to delivering drone routing.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 2 .
Figure 2. Optimal routes of TSP.(a) Total flight distance and (b) total flight time.

Figure 3 .
Figure 3. Optimal routes of FSTSP.(a) Total flight distance and (b) total flight time.

Figure 4 .
Figure 4. Forces exerted on a drone.(a) Flight without payload and (b) flight with payload .

Figure 4 .
Figure 4. Forces exerted on a drone.(a) Flight without payload and (b) flight with payload w.

Figure 4 .
Figure 4. Forces exerted on a drone.(a) Flight without payload and (b) flight with payload w.
),  is 0.000015.Regarding the weight w of the luggage, SkyLift is within the range of 0