Delivery Drone Routing under Load dependent Flight Speed based on Integer Quadratic Programming Mao Nishira 1, Hiroki Nishikawa 2, Xiangbo Kong 1, and Hiroyuki Tomiyama 1 1 Graduate School of Science and Engineering, Ritsumeikan University, Shiga, Japan 2 Graduate School of Information Science and Technology, Osaka University, Osaka, Japan Abstract Drone delivery has been attracting attention, and many companies are focusing on this as a means to solve logistics issues. In order to deliver packages by drone, it is necessary to optimize the route in advance. So far, drone routing problems that consider changes in flight speed due to load have been proposed. This routing problem is called the FSVRP. An Integer Linear Programming (ILP) approach to the FSVRP was proposed in the past. In this paper, to obtain a more accurate solution than the ILP approach, we propose an Integer Quadratic Programming (IQP) approach. Experiments show that the proposed method is effective compared to state- of-the-art methods. Keywords 1 Delivery drones; Vehicle routing problem; Flight speed-aware vehicle routing problem; Integer quadratic programming; Quadratic Approximation 1. Introduction The number of parcel deliveries in logistics is increasing every year. There are three major issues in logistics. The first is a shortage of drivers and other labor. The second is traffic congestion caused by trucks. The third is the inefficiency of redelivery. As the number of parcels handled increases, so does the number of redeliveries. Recently, unmanned aerial vehicles (UAVs), or drones, have emerged as a promising vehicle for parcel delivery. Compared with delivery by trucks and motorcycles, delivery drones have several advantages in terms of CO2 emission, labor cost and traffic congestion. Over the past few decades, delivery drones have been studied from various aspects, and routing for multi-destination delivery is one of the most important research topics [1,2]. The routing problem for delivery drones is a special case of vehicle routing problems (VRPs) which have been studied for many years. Traditionally, VRP has been studied mainly for delivery trucks [3,4]. However, these studies cannot be directly applied to delivery drones since delivery drones has to take care of additional constraints such as battery capacity, the weight of parcels, weather conditions, and so on. Unlike trucks, most drones are powered by battery, and their flight time (or distance) is severely limited. Moreover, drones are much smaller than trucks and trains in size, which limits the maximum payloads to carry at a time. Recently, a number of researchers have addressed delivery drone routing. For example, Poikonen et al. proposed a drone routing problem which aims to minimize the delivery completion time [5]. The works in [6,7,8] also studied delivery drone routing problems. The common assumption of these works is that the flight speed of drones is a constant independent of the payloads. Funabashi et al. recently proposed a new routing problem for delivery drones, called FSVRP (flight-speed aware vehicle routing problem), that asks the route with the minimum flight time considering payload-dependent flight speed [9]. They also developed a dynamic programming algorithm for FSVRP. The algorithm runs fast, but is not extensible. If new constraints are added, the algorithm may need substantial changes. The authors of [10] proposed an approach to delivery drone routing based on integer linear programming (ILP), where a general-purpose mathematical programming solver IBM CPLEX was used. In fact, the objective The 4th International Symposium on Advanced Technologies and Applications in the Internet of Things (ATAIT 2022), August 24-26, 2022, Ibaraki, Japan © 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) function of the original FSVRP is not linear but very complex. The work in [10] linearizes the objective function by the least-squares approximation method at the cost of degraded accuracy. This paper proposes an integer quadratic programming (IQP) approach to FSVRP. The objective function of this work is more accurate than the one in [10], and thus, this work can find better routes. The rest of this paper is organized as follows. Section 2 describes the FSVRP. Section 3 describes an IQP approach to the FSVRP. Section 4 shows the experimental evaluation. Finally, Section 5 concludes this paper. 2. FSVRP This section describes FSVRP (flight-speed aware vehicle routing problem) presented in [6] and [10]. FSVRP is a variation of routing problems for delivery drones that accounts for flight speed as a function of load. The FSVRP is the problem of finding a route that minimizes the time taken to deliver all parcels to their respective customers and back to the depot. The problem which we are concentrating on requires the shortest route in terms of flight time, taking into account the speed of flight as a function of the load. 2.1. An Example As an example, we depict Figure 1. The node labeled "0" is the depot and the other nodes are the customers. The number in the box indicates the weight of the parcel to be delivered to the customer. The red number on the sides is the distance between the two customers or a depot and a customer. Drones can ignore geographical constraints since the drones can fly, and the drones can visit any customer. Figure 1, if we minimize the total flight distance, this problem is typically solved as a traveling salesman problem. In figure 2, the route shown by the arrow is obtained as the optimal route. The total flight distance of the shortest route in the TSP is 146 (= 43 + 52 + 24 + 27). The total flight time of this TSP optimal route is 102 (= 36 + 36 + 15 + 15) unit of time in Figure 2 (b). This route is the optimal route in terms of the total flight distance. If the flight speed is fixed, the shortest route is optimal in terms of both flight distance and flight time. However, the heavier the load, the slower the drone flies, and the lighter the load, the faster the drone flies. In figure 3, the route shown by the arrows is obtained as the shortest route in terms of flight time, and the shortest route is determined by two factors: the weight of the load and the distance. The total flight distance of the FSVRP-optimal route is 146 (= 27 + 24 + 52 + 43), which is the same as the TSP optimal route. The total flight time of the FSVRP optimal route equals the total flight time of that TSP optimal route, but in general, the total flight time of the FSVRP equals or exceeds the total flight time of the TSP. However, the total flight time of the FSVRP-optimal route is 93 (= 21 + 15 + 33 + 24), which is shorter than the TSP optimal route. In Figure 2, the drone has all the parcels at the depot, and its total weight is 127. It starts delivery with a total weight of 127 and delivers the load to customer 2, whose load is 43. The remaining parcels with a total weight of 84 are then delivered to each customer. In Figure 3, the drone starts the delivery with 127 (units of weight) and delivers the heaviest parcel with a load of 50 to customer 1. It then delivers the remaining parcels with a total weight of 77 to each customer. Then, the drone can fly faster during the rest of the delivery. As shown by this arrow, in the FSVRP with drones, there are cases where it is better to carry the heavy load first, because the drone's flight speed is highly dependent on the load. 2.2. Formulation In this problem, there are N customers to be delivered. To avoid loss of generality, a drone delivers only one parcel per customer. Therefore, the number of parcels is N. The customers are numbered from 1 to N. The parcel to be delivered to customer (1 ) is called parcel . The depot is numbered 0 as shown in Figure 1. In this paper, the drone delivers all parcels to customers in a single delivery. All the parcels are uploaded onto the drone at the depot, which is the delivery base, and the drone begins the delivery. If the Figure 1. An example problem [10] (a) Distance (b) Time Figure 2. Optimal route of TSP [10] (a) Distance (b) Time Figure 3. Optimal route of FSVRP [10] total weight of the parcels to be delivered exceeds the maximum payload of the drone, the parcels must be partitioned into multiple groups and delivered multiple times before delivery, which is not considered in this paper. To begin with 1, 2 represents the distance between points 1 and 2. Also, is the decision variable in the routing problem and indicates the -th customer to visit. The equation (1) represents the drone departs from the depot and begins its delivery. In addition, the drone returns to the depot again when it finishes the delivery. 0 1 0 (1) The drone visits every customer only once. This is defined as follows: 1 1 (2) 1 , , (3) Let is the weight of parcel and is the total payload when the drone leaves -th visited point. In Figure 4, is the weight of the drone itself. The equation (4) represents when the drone begins its delivery from the depot, all parcels are loaded onto the drone. (a) Fight without payload (b) Flight with payload Figure 4. Resultant forces on a drone 0 (4) When the drone arrives at the -th stop (1 ) at point , the drone unloads a parcel of weight for the customer. From the above, the equation (5) represents the total weight loaded on the drone when departing point j. 1 (5) Besides, 1, 2 represents the time it takes for the drone to fly between points 1 and 2. Note that is a function of the load. This equation (6) shows the time between the -th and 1 -th visited points, by dividing the distance by the flight speed. , 1 , 1 / (6) The equation (7) represents objective of the FSVRP. The FSVRP asks the route with the shortest total flight time for a delivery. Min , 1 / (7) The flight speed of a drone is also greatly affected by the weight of the parcel it is carrying, wind, and other factors. When a drone flies, it tilts its fuselage. in the equations (8) to (13) represents lift, and represents the weight of the drone itself. In addition, indicates gravitational acceleration. In order for the drone to fly without falling, (i.e., the vertical component of ) must be equal to the gravity . In equations (8) to (13), is the pitch angle when equals gravity. (i.e., the horizontal component of ) is the force toward the customer making the delivery. It is equal to the air resistance 0 . Here, is a coefficient specific to the drone. cos (8) sin 0 (9) The equations (8) and (9) show the force balance when the drone is not carrying a load, while the equations (10) and (11) show the force balance when the drone is flying with a load. In other words, is bigger than . Hence must be greater than . Then must be less than and must be less than . cos (10) sin (11) From the above, and can be expressed as the equations (12) and (13). arccos / (12) arccos / (13) 3. An IQP Approach to FSVRP Section 3 explains how the problem presented in Section 2 is solved using the IP solver. IP solvers, like CPLEX and Gurobi, utilize a combination of a variety of heuristics on their own to quickly find a solution. Furthermore, many solvers are frequently built to work on multicores. As a result, solvers are able to solve a variety of problems more quickly than an algorithm designed by programmers who are not algorithmic experts. As the objective of the problem in Section 2 is not linearized and included trigonometric functions, it is not possible to solve the problem with ordinary IP solvers. This section provides an approximated objective by quadratic approximation to solve the problem using an IP solver. Based on Section 2, we can obtain equation (14), which is a sort of approximation, but it is still too complex to be solved as an integer programming. sin arc cos / (14) The equation (15), which is the reciprocal of the flight speed, is used for the objective of the FSVRP. 1/ 1/ sin arc cos / (15) The vertical axis in figure 5 is the reciprocal of the flight speed, and the horizontal axis is the weight of loads delivered by the drone. Currently, drones are being developed to carry cargo, and the weights of the cargo that drones can carry vary. In this figure, the maximum payload a drone can carry is based on SkyLift [11]. SkyLift is the name of the cargo drone. The literature in [11] states that the maximum payload of this cargo drone is 30 kg, but we experimented with a maximum payload of 27 kg to give ourselves a margin of 10%. The total weight of the parcels are determined in the range of 0 27000. For the black curve, a quadratic approximation using the least-squares method gets a red curve in the figure. The R-squared value, which is the error from equation (16) due to the approximation, is 0.9647. 1/ 3 10 3 10 0.1126 (16) Substituting equation (16) into the expression for velocity in equation (7) in section 2, the problem can be solved with the IP solver. Min , 1 3 10 3 10 0.1126 (17) Figure 5. Quadratic regression based on least squares method 4. Experiments The effectiveness of our proposal is evaluated through experiments. Six routing algorithms shown below are compared in terms of total flight time and the runtime of the algorithms. DP-FSVRP: Solve the FSVRP by dynamic programming. It calculates the route that minimizes the total flight time. DP-TSP: Solve the TSP by dynamic programming. It calculates the route that minimizes the total flight distance. CPLEX-ILP-single: Solve the linearized FSVRP with CPLEX on a single thread. CPLEX-ILP-multi: Solve the linearized FSVRP with CPLEX on multiple threads. CPLEX-IQP-single: Solve the FSVRP based on an IQP approach with CPLEX on a single thread. This is our proposal. CPLEX-IQP-multi: Solve the FSVRP based on an IQP approach with CPLEX on multiple threads. This is also our proposal. In our experiments, we have prepared five problem instances for each method with 5 to 20 customers, which are randomly generated. We set the runtime limited up to 3600 seconds in wall-clock time. The experiments are conducted on AMD Ryzen 7 PRO 4750G (8 cores, 16 threads) and 64GB of main memory. In addition, to compute the solution, we used Python 3.8.5 for the existing method and IBM ILOG CPLEX Optimization Studio 20.1 for the proposed method. 4.1. Runtime Evaluation In Figure 6, it shows comparison in terms of runtime of each method using a logarithmic scale. The vertical axis is runtime, and the horizontal axis is the number of customers. The optimal solutions of CPLEX-IQP-single and CPLEX-IQP-multi can be found because it takes so much runtime to compute a quadratic function. From 5 to 16 customers, DP-FSVRP is the fastest in finding the optimal route compared to other algorithms since CPLEX incurs the overheads by generating the model and parallelizing it into multiple threads. We think that the reason why there are variations of the solution time of CPLEX is because the number of instances is as small as 5. More than the number of customers 17, there are cases in which CPLEX-ILP can require to solve faster than DP-FSVRP. Due to the multi- threaded execution, CPLEX-ILP-multi and CPLEX-IQP-multi can achieve shorter runtime than CPLEX-single. These results show that linear approximations can be solved faster than quadratic approximations. Since the runtime of CPLEX-IQP is about 1 minute for up to the number of customers 9, we do not consider this to be a particular problem in practical use. 4.2. Routing Results Figure 7 shows the total flight time of each method. The vertical axis is the total flight time normalized to DP-FSVRP, and the horizontal axis is the number of customers. The experimental results show that for many problem instances, CPLEX gets the same optimal routes as DP-FSVRP. Figure 7 shows DP-TSP can hardly find good routes in terms of delivery time since DP-TSP solves the problem that minimize total flight distance. Actually, CPLEX-ILP was worse than DP-FSVRP by an average of 0.73%. CPLEX-IQP was worse than the DP-FSVRP by an average of 0.83%. As mentioned above, the runtime is limited to 3600 seconds, so more than the number of customers 11, CPLEX-IQP cannot obtain optimal solutions. One possible reason why CPLEX-IQP showed worse results than CPLEX- ILP is that the drone tries to move as long as possible between customers at the time when the drone is flying fastest. Since CPLEX-IQP selects a route different from the optimum route obtained by DP- FSVRP and CPLEX-ILP, it is considered that the optimum solution cannot be obtained. Figure 6. Logarithmic scale of runtime DP-FSVRP DP-TSP CPLEX-ILP CPLEX-IQP 1.15 1.1 1.05 1 0.95 0.9 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Number of customers Figure 7. Routing results of total flight time 5. Conclusion This paper presented an IQP approach to delivery drone routing under load-dependent flight speed. Up to the number of customers 9, the runtime of CPLEX-IQP is at most 1 minute, and since CPEX- IQP provides a better solution than CPLEX-ILP, we think that our proposed method is useful in practical use. We also think that CPLEX-IQP is superior to dynamic programming for FSVRP based on the fact that IP solvers can easily extend the problem. In the future, we will take into account multiple deliveries to routing problems of a drone considering load-dependent flight speed. 6. Acknowledgements This work is partly supported by JSPS KAKENHI 19H04081, 20H00590 and 21K19776. 7. References [1] Zheng Wang and Jiuh-Biing Sheu, “Vehicle routing problem with drones,” Transportation Research part B: methodological, vol. 122, pp. 350-364, 2019. [2] Chun Cheng, Yossiri Adulyasak, and Louis-Martin Rousseau, “Drone routing with energy function: Formulation and exact algorithm,” Transportation Research Part B: Methodological, vol. 139, pp. 364-387, 2020. [3] Chase C. Murray and Amanda G. Chu, “The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery,” Journal of Transportation Research, vol. 54, pp. 86-109, 2015. [4] Kevin Dorling, Jordan Heinrichs, Geoffrey G. Messier, and Sebasian Magierwoski, “Vehicle routing problems for drone delivery,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 47, no. 1, pp. 70-85, 2016. [5] Stefan Poikonen, Xingyin Wang, and Bruce Golden, “The vehicle routing problem with drones: Extended models and connections,” Networks, vol. 70, no. 1, pp. 34-43, 2017. [6] Yusuke Funabashi, Ittetsu Taniguchi, and Hiroyuki Tomiyama, “Work-in-progress: Routing of delivery drones with load-dependent flight speed,” IEEE Real-Time Systems Symposium (RTSS), pp. 520-523, 2019. [7] Patchara Kitjacharoenchai, Mario Ventresca, Mohammad, Mohammad Moshref-Javadib, Seokcheon Leea, Jose M. A. Tanchocoa and Patrick A. Brunese, “Multiple traveling salesman problem with drones: Mathematical model and heuristic approach,” Computers & Industrial Engineering, vol. 129, pp. 14-30, 2019. [8] Patchara Kitjacharoenchai and Seokcheon Lee, “Vehicle routing problem with drones for last mile delivery,” Procedia Manufacturing, vol. 39, pp. 314-324, 2019. [9] Ho Young Jeong and Seokcheon Lee, “Optimization of vehicle-carrir routing: Mathematical model and comparison with realted routing models,” Procedia Manufucturing, vol. 39, pp. 307-313, 2019. [10] Mao Nishira, Satoshi Ito, Hiroki Nishikawa, Xiangbo Kong and Hiroyuki Tomiyama, “An ILP- based Approach to Delivery Drone Routing under Load-dependent Flight Speed,” International Conference on Electronics, Information, and Communication (ICEIC), pp. 1-4, 2022. [11] SkyDrive Inc., SkyLift transports your payload to wherever you need it most., 2020. [Online]. Available: https://skydrive2020.com/eng/cargo-drone/spec.html#spec-04 [Accessed On Jun., 2022].