Drones Routing with Stochastic Demand

: Motivated by the increasing number of drones used for package delivery, we ﬁrst study the problem of Multiple drOne collaborative Routing dEsign (MORE) in this article. That is, given a ﬁxed number of drones and customers, determining the delivery trip for drones under capacity constraint with stochastic demand for customers such that the overall expected traveling cost is minimized. To address the MORE problem, we ﬁrst prove that MORE falls into the realm of the classical vehicle routing problem with stochastic demand and then propose an effective algorithm for MORE. Next, we have a scheme of resplitting customers into different individual delivery trips while the stochastic demands are determined. Moreover, we consider a variety of MORE, MORE-TW, and design an effective algorithm to address it. We conduct simulation experiments for MORE to verify our theoretical ﬁndings. The results show that our algorithm outperforms other comparison algorithms by at least 79.60%.


Introduction
Boosted by the rapid development of technology in communication and artificial intelligence, the market for drones is growing rapidly. The global commercial drone market size was valued at 13.44 billion dollars in 2020; it is expected to expand at a compound annual growth of 57.5% from 2021 to 2028 [1]. Due to the benefit of drones, such as small size, lower construction cost, and lower fuel consumption, drones have been utilized in many areas, such as aerial photography, search and rescue operations, agriculture, shipping and delivery, engineering, 3D mapping, safety surveillance, wireless internet access, research and nature science, and so on [2][3][4][5][6][7][8][9][10][11][12][13][14].
Recently, drones have been widely used in parcel delivery areas [15]. Large companies such as Amazon, Walmart, UPS, and Google have invested in drone delivery projects. The US Federal Aviation Administration (FAA) has approved Amazon's drone delivery program. The drone delivery of Amazon can connect rural populations and deliver packages in 30 min or less via Amazon Prime Air [16]. UPS was granted the FAA approval to operate the UPS Flight Forward drone airline [17]. Zipline delivers blood and then expands to PPE when COVID-19 hits [18]. Google's drones have finished burritos and dog food delivery in Australia, baked goods in Finland, and recently groceries in the United States. According to the new market research report, the drone package delivery market grew from USD 528 million in 2020 to USD 39,013 million [19]. With the growth of the market, an increasing number of drones are used for package delivery.
Drone routing is a basic problem in package delivery, and drone routing has the following advantages [20]: • Drone routing for package delivery can save labor as no drivers are needed • Drone routing for package delivery can provide faster and on-time delivery service for customers • Drone routing for package delivery is not limited by the complex road network Motivated by the above advantages of drone routing for package delivery, we embark in this paper on a systematical formulation and analysis of the Multiple drOne collaborative Routing dEsign (MORE) to increase delivery efficiency as shown in Figure 1. Specifically, we consider a generic drone delivery scenario with a set of customers distributed on a 2D Euclidean plane having different demands every day, termed stochastic demands. Our focus on stochastic demands for customers is motivated by the practical parcel delivery scenarios including but not limited to: • Drone contactless delivery for nucleic acid test samples, plasma, medicines, and other medical supplies; • Fresh food and daily groceries delivery for customers every day; • Takeaway delivery.
In all these applications, customer demands vary from day to day and are volatile. It is also impossible to forecast customer demands before a few days. Hence, we seek to determine the adjusted delivery routes for most m drones to meet the actual demands of customers with the goal of minimizing the total expected delivery cost subject to the constraint that the total demands of served customers in one delivery trip of a drone should be no more than its capacity.
To the best of our knowledge, the above multiple drone collaborative delivery route design with stochastic demand problems has never been studied in the literature despite its practical importance. More specifically, there is literature [20][21][22][23][24][25][26][27] that considers multiple drone delivery problems, but none of them consider the stochastic issues for customer demands. There have emerged some works studying truck and drone collaborative delivery problems [28][29][30][31][32][33][34]. However, most of them consider the collaborative scheduling between drones and trucks, which is fundamentally different from ours and cannot be applied to address MORE. In addition, some other work considers the problem of vehicle routing with stochastic demand [35][36][37][38], which is shown to be mostly related to ours. Nevertheless, none of them can be directly used to address MORE. The root reason is that we consider the property of flexibility and capacity constraint for each drone.
To fill the gap in drone delivery routing design with the stochastic demands of the state-of-the-art, we make the following contributions in our paper:

•
We first give the definition of MORE-the collaborative delivery route design with stochastic demands for customers.

•
We develop a comprehensive algorithm to solve MORE and give the performance analysis for the proposed algorithm.
Our MORE problem has three main technical challenges. The first challenge is that MORE is NP-hard. MORE is NP-hard because MORE can be polynomially reduced to the realm of the vehicle routing problem with stochastic demand, which has proved to be NP-hard. The second challenge is to design an effective algorithm to find delivery trips for drones with stochastic and capacity constraints. The third challenge is that it is difficult to determine the number of individual delivery trips for a drone with limited capacity due to the uncertain property of demands.
To address the first challenge, we reduce our MORE to the realm of the classical vehicle routing problem with stochastic demand. To address the second challenge, we design a new algorithm for MORE. Most importantly, we prove the performance of our proposed algorithm for MORE. To address the third challenge, we propose an algorithm scheme to reduce the coupling of capacity constraint and stochastic demand.
To verify the performance of our proposed algorithm for MORE, we conducted extensive simulations. The results show that our algorithm outperforms other comparison algorithms by at least 79.60%.

Related Work
This review focuses on problems involving the coordinated use of trucks and drones for making deliveries, the use of only drones for parcel delivery, and the vehicle routing problem with stochastic demand [29].
Truck-drone delivery problems. The truck and drone collaborative delivery problem is closely related to ours [39][40][41][42][43][44][45][46][47][48][49]. However, it is mainly mixed integer linear programming as it considers the collaborative relationship between drones and trucks, which is fundamentally different from our problem. Chase C. Murray et al. [28] consider the scenario in which an unmanned aerial vehicle works in collaboration with a traditional delivery truck to distribute parcels. They present mixed integer linear programming formulations for the proposed flying sidekick traveling salesman problem (FSTSP) and parallel drone scheduling TSP (PDSTSP), along with two simple yet effective heuristic solution approaches. Later, Chase C. Murray et al. [29] extend to the situation in which a delivery truck and a heterogeneous fleet of drones coordinate to deliver small parcels. Chase C. Murray et al. proposed a heuristic solution approach for mFSTSP that consists of solving a sequence of three subproblems. The authors in [30] address the tradeoffs between speed and range in a variant of a delivery truck and the drone's coordinated delivery problem where the speeds of the drone are treated as decision variables. In [31][32][33][34], the authors consider using multiple trucks and multiple drones to deliver parcels. The authors in [50] model two different types of drone tasks (drop and pickup) scheduling problems as an unrelated parallel machine scheduling with sequence-dependent setup, precedence-relationship, and reentrant. In [34], the authors propose the Drone Assignment and Scheduling Problem (DASP) to look for an optimal assignment and schedule of drones with the goal that the makespan is minimized. In [49], the authors consider the uncertain set of pickup requests with deadlines to maximum working hours using trucks and drones. They propose a heuristic approach to solve their problems. However, our study considers the uncertain demands of customers to minimize the delivery distances using only drones. We propose a greedy algorithm to solve our problem. In [51], the authors consider the impact of the wind to plan minimum-energy trajectories for the drone to provide service for the customers from starting and returning to the truck. In [52], the authors introduce the docking hub to improve the service coverage of load-dependent drones. They design a branch-and-price-and-cut algorithm to solve the mixed-integer problem.
Drone-only delivery problems. There are many works that consider multiple drone delivery problems, but none of them take into account the stochastic issue regarding customer demands. Kevin et al. [21] propose two multi-trip VRPs for drone delivery that address the issues of the effect of battery and payload weight on energy consumption.
They design a simulated annealing heuristic to solve the problem. However, they do not consider the stochastic demands. The work [22] is to create a decision-making tool for the design of a drone fleet in the case of forecast deliveries over a time horizon under operational constraints. The difference between our work with the work of [22] is that our work considers the demands of customers to be stochastic, while the work of [22] regards the demands of the customers as constant. The objective of [22] consists of minimizing the total delivery distance and the total batteries, while our work focuses on minimizing the total delivery cost. Figliozzi et al. [23] present a novel analysis of lifecycle UAV and ground commercial vehicles CO 2 emissions. The difference between our study with study [23] is that the [23] focuses on the analysis of lifecycle UAV and ground commercial vehicles in terms of CO 2 emissions, while our study makes more attention to solving the combinational problem.
Study [24] shows that the long traveling distances of drones per package greatly increase the life-cycle impacts. They focus on minimizing extra warehousing and limiting the size of drones, while our objective is to plan delivery routes for drones to reduce the cost. Moreover, there is also literature [20,25] considering drone energy consumption. They focus on designing exact algorithms, while we can not propose an exact algorithm due to the uncertainty of the demands of customers. In [53], authors want to find the distribution point for drones in the scenario that drones must follow up the open space above the road while our scenario is that drones can freely fly on straight lines between customers. In [54], authors jointly optimize the route and sensing task to minimize the energy consumption of the drone and maximize sensing reward. The authors propose the Bellman-Ford algorithm to solve their problem. In [55], the authors design a centralized framework to deal with enormous collisions and can obtain a collision-free path for fast delivery. However, our study deals with flying collisions by means of splitting customers into different groups; only one drone can serve the customers in the same group. In [56], the authors investigate the attention-based pointer network method to obtain multitrips for a drone for parcel delivery.
Vehicle routing problem with stochastic demand. The related problem of vehicle routing with stochastic demand has been well studied [35][36][37][38][57][58][59]. The proposed problem in this paper distinguishes itself from this one by taking into account flexible drones as well as capacity constraints for each drone. There is literature [35,36] studying this problem with the assumption that the stochastic customer demand can only take values 0 or 1. Lei et al. [37] consider the case that routes can be paired and the customer demand can be split at the shared location. The authors in [38] combine ideas from vehicle routing and manufacturing process flexibility to design an overlapped routing strategy before vehicles are dispatched.

Assumptions
We make some assumptions to help us to build a reasonable formulation about drone routing with stochastic demand.

•
All drones start flying from the depot, deliver several parcels for customers, and need to return to the depot to load new parcels with recharging. • The depot has enough batteries for drones to provide energy. • All drones are homogeneous; that is, all drones have the same limited capacity, flight speed, and so on. • The sum of demands of served customers by a drone in one trip can not exceed the capacity of the drone. • The stochastic demand for each customer is independently and identically distributed. • The flight speed of all drones is constant, and we do not consider the weather factor affecting flight speed. • The time to replace batteries, load parcels at the depot, and unload parcels for customers is ignored.
• The takeoff and landing processes of all drones for different customers are the same.

Drone Routing Cost Model
Suppose there are n customers denoted as C = {c 1 , · · · , c n } in a 2D area with fixed known positions. We use a directed graph G = (C 0 , E) to represent the drone routing cost where C 0 = {c 0 } ∪ C with c 0 representing the depot of all drones, and We use the travel distance from customer c i to customer c j as a routing cost u(c i , c j ) associated with the edge (c i , c j ) ∈ E. We consider the drone flying from customer c i to customer c j as shown in Figure 2 where h i is the end point of the take-off process and h j is the start point of the landing process. The routing cost u(c i , c j ) is the sum of distances of the take-off process, flying directly from h i to h j , and the landing process. Figure 2. Drone routing cost model.
Let D i ∼ D(µ, σ), (i = 1, · · · , i = n) be the random daily demand of customer c i . We use µ and σ to denote the mean and standard deviation of the daily demand for each customer. We assume all customer daily demands are revealed before the drones are sent out on each day.
There are m drones to be sent out to be responsible for meeting the daily demands of all customers. We define a delivery trip of the drone i denoted as T =< c 0 , c i , c i+1 , · · · , c j , c 0 > where a drone departs from the depot c 0 , serves a sequence of customers from customer c i to customer c j , and then returns the depot c 0 . The delivery trip cost for the drone i can be represented as follows.
where u(c 0 , c i ) and u(c j , c 0 ) are the traveling distance between the depot c 0 and the served customer, u(c m , c m+1 ) denotes the delivery distance of the neighboring customer served by the drone i. In our model, we do not restrict the number of the delivery trip that the drone i can fly. We define Γ i as the set of all the delivery trips for the drone i. We assume each drone has limited capacity Q where it can only serve a part of customers with a total demand of Q in one delivery trip before returning to the depot c 0 for refilling the capacity. Due to the limited capacity for one delivery trip, the customer demand can be met by multiple drones. Table 1 lists the main notations and symbols used in this paper. Table 1. Main notations and symbols used in the paper.

Problem Formulation
We consider the problem of determining the delivery trip for at most m drones towards minimizing the total expected delivery trip cost of all of the drones under the constraint that the total demand of all the served customers in one delivery trip should be no more than the limited capacity Q of the drone. We formulate the Multiple drOne collaborative Routing dEsign problem for meeting the daily stochastic demands of all customers under capacity constraint (MORE) as follows.
whereD j is a kind of realization of the random daily demand for customer c j . Our optimization goal is to determine the delivery trip set Γ i for each drone while meeting the daily random demand of all customers to minimize the total expected delivery trip cost of all of drones. The MORE problem (P1) is NP-hard.
Proof. The MORE problem can be polynomially reduced to the realm of the classical vehicle routing problem with stochastic demand, which has proved to be NP-hard.
Therefore, our MORE problem is NP-hard.

Algorithm Overview
In this subsection, we describe our algorithm via a toy example. Our algorithm as shown in Algorithm 1 is mainly divided into four steps. Figure 3 shows the distribution of customer positions and drone positions; all of the drones are located in depot station c 0 , waiting to deliver commodities to customers. The demands of all customers are not revealed for drones. First, we approximately solve the traveling salesman problem (TSP) with the goal of minimizing the total traveling cost. We relabel the customer according to the order in the TSP tour by picking an arbitrary customer and orientation in the TSP as the starting customer and orientation. For simplicity, we pick the starting customer and its orientation in the TSP tour as our first customer and orientation in the relabeled. We obtain the relabeled customer set {c 1 , · · · , c 12 }. Second, we divide the customers into three groups where the number of groups is determined by the number of drones. We calculate the number of customers in each group N = 4. Then, the three group sets of customers are Φ 1 = {c 1 , · · · , c 4 }, Φ 2 = {c 5 , · · · , c 8 }, and Φ 3 = {c 9 , · · · , c 12 }, respectively. We assign the three groups customers to the drones as their primary serving customer set as shown in Figure 4. Third, we set the overlapped number of the serving customer for the drone as K = 1. We can obtain the extended serving customer set for each drone as shown in Figure 5. Fourth, we use Algorithm 2 to find the delivery trip for each drone to meet the demand of all customers as shown in Figure 6.

Approximation Algorithm for MORE
As MORE is NP-hard and difficult to be directly tackled using the precise algorithm, in this section we first design an approximation algorithm to solve MORE. Then, we theoretically analyze the effectiveness of our proposed algorithm for MORE.

Solution for MORE
Our approximation algorithm described in Algorithm 1 for MORE consists of four steps. First, we find the shortest path that visits all of the customers and returns to the depot by solving the classical Travelling Salesman Problem (TSP). we relable the customers according to the order and the orientation in the TSP tour. We call this step relabelling customers. Second, we divide the relabelled customers into disjoint groups. Each group is served by one drone. We call this step primary customers assignment. Third, we assign the overlapped served customers for each drone, which we call extended customers assignment. Finally, we design a greedy algorithm to determine the delivery trips for each drone, which we call finding delivery trips. Next, we give a description in detail of these four steps one by one in the following four subsections.

Relabelling Customers
In this step, we only consider that a drone starts from the depot c 0 , passes through all the customers in customer set C once, and finally returns to the depot c 0 with the goal of minimizing the total traveling distance, which can be transformed into the TSP problems. Therefore, we formulate the TSP problem as follows.
Algorithm 1 General Algorithm for MORE.
Input: Customer C, drone capacity Q, observed daily demandD 1 , · · · ,D n . Ouput: A sequence of serving customers for individual delivery trips. Solve TSP problem with the instance C 0 as input to obtain new customer sequence; Divide new customer sequence into m disjoint groups to obtain primary customer set; Assign extended customers to obtain extended customer set; Call Algorithm 2 with the instance ( Φ 1 · · · Φ m , Ψ 1 · · · Ψ m , Q,D 1 · · ·D n ) as input to obtain the delivery trip; Output the individual delivery trip.
We use the cutting-plane method, iteratively solving linear programming relaxations of the problem P2. As shown in Figure 7, we can obtain the TSP tour by solving problem P2. We pick an arbitrary customer and its orientation in the TSP tour as our starting customer and orientation. Then, we relabel all customers according to their order in the TSP tour from the starting customer. For simplicity, we pick the first customer and its orientation in the TSP tour as our starting customer and orientation. Then, we relabel all customers and obtain the new number for each customer. Note that we use the new number of all customers in the next section.

Primary Customers Assignment
At this step, we divide n customers into m disjoint groups. Each group is served by one drone. That is, each drone is primarily assigned to serve N = n m customers. We use Φ j to denote the primary serving customer set for drone j. According to the new number of customers, we set the primary serving customer set for drone j to be Φ j = {c (j−1)N+1 , · · · , c jN } for 1 ≤ j ≤ m. Each primary serving customer set contains exactly N customers, and each customer is contained in exactly one primary serving customer set. While the primary customer assignment is simple and intuitive, it fails to coordinate drones to improve the collaborative effectiveness in dealing with random daily demands. This motivates us to consider extended customers assignment as the next step.

Extended Customers Assignment
To improve the collaborative effectiveness of multiple drones to meet the random daily demands, we consider the k-overlapped serving customers of drones, where k is the parameter. k can take any integer value from 0 to n. We define the extended serving customer set as the primary serving customers plus k additional customers. We use Ψ j to denote the extended serving customer set for drone j. We set extended serving customer set for drone j to be Ψ j = {c (j−1)N+1 , · · · , c min{jN+k,n} } for 1 ≤ j ≤ m − 1, and Ψ m = Φ m . For k = N, the extended serving customer set of drone j is the concatenation of the primary serving customers for adjacent drones (drone j and drone j + 1). Therefore, the extended serving customer set of drone j and j + 1 overlap at Φ j+1 , which leads to collaboratively assigning drones to meet the random daily demand for customer set Φ j+1 . In executing the delivery for customers at the next step, our strategy uses drone j to ensure that the demands of all customers in the primary serving set Ψ j are met while using the remaining capacity of drone j to meet the demands of customers in the additional segment of the extended serving customer set Φ j .

Finding Delivery Trips
At this step, we design our strategy by sequentially determining delivery trips for drones from 1 to m. We follow three basic principles while finding delivery trips for drones. First, each drone starts a delivery trip from depot c 0 with full capacity Q and ends a delivery trip by returning to depot c 0 . Second, the delivery trip of drone j must first satisfy the demand for all customers in its primary serving customer set Φ j that have not been satisfied by the drone j − 1, and skip customers with zero demand. Third, if there is a remaining capacity of drone j after the demands for all customers in its primary serving customer set Φ j are met, drone j will continue to use its remaining capacity to satisfy the demands in the additional segment of its extended serving customer set Ψ j .
We define some notations used in designing our algorithm for finding delivery trips. For each drone V j , c s j represents the first customer that it serves and c e j represents the last customer that it serves; w j represents the amount of demand in its primary serving customer set Φ j ; e j is the remaining capacity after satisfying the demands in its primary serving customer set Φ j ;ẽ j denotes the portion of the remaining capacity e j that will be used to satisfy the demands in the additional segment of its extended serving customer set Ψ j . We illustrate our algorithm of finding delivery trips in Algorithm 2. Generally, we divide Algorithm 2 into three modules. In the first module, we determine the portation of the remaining capacity e j that will be used to satisfy the demands in the additional segment of its extended serving customer set Ψ j and the total carried demand of the drone. In the second module called the finding served customer module, for each drone V j , it finds the first served customer c s j using capacityẽ j−1 in primary serving customer set Φ j and the last served customer c e j using capacityẽ j in extended serving customer set Ψ j . In the third module called the determining individual delivery trip module as shown in Algorithm 2, the drone with the total carried demand serves the customers starting from the first served customer c s j in its first delivery trip, and it skips the customers with zero demand. It will stop serving customers in its first delivery trip in the three situations. First, when the total demands of served customers have achieved the total carried demand of the drone, the drone finishes its first delivery trip. Second, when the total demands of served customers have achieved their capacity, the drone finishes its first delivery trip. If there are remaining customers to be served, the drone can start a new delivery trip after returning to the depot c 0 . Third, when all customers, from the first customer to the last customer, have been served, the drone finishes its first delivery trip.

Algorithm 2 Algorithm for Determining Delivery Trips.
Input: Primary serving customer set Φ 1 · · · Φ m , extended serving customer set Ψ 1 · · · Ψ m , drone capacity Q, observed daily demandD 1 · · ·D n . Ouput: Delivery trip for each drone. Calculate the total demandD Φ 1D Φ m of each primary serving customer set; Calculate the total demand of the additional segment of the extended serving customer set; Set the total filled demand of drone V j F j = 0 else Setẽ j = min(Q * w j /Q − w j ,D Ψ j ) Set F j =w j +ẽ j Finding the first served customer c s j using capacityẽ j−1 in primary serving customer set Φ j ; Finding the last served customer c e j using capacityẽ j in extended serving customer set Ψ j ; Call Algorithm 3 with the instance (c s j , c e j , F j ) as input to obtain the delivery trips. Output the delivery trip for each drone.

Theoretical Analysis for MORE
In this subsection, we theoretically analyze the performance of our proposed algorithm for MORE. Lemma 1. We have the following: where Z * denotes the travel cost of reoptimization.

Lemma 2.
We have the following: where L * denotes the optimal TSP distance traveling through all customers.
Theorem 1. Suppose the relabelling customer sequence [c 0 , c 1 , · · · , c n , c 0 ] forms an α-optimal TSP tour. Let Z and r avg denote the expected total delivery distance and the average number of delivery trips per drone under our MORE algorithm for the fixed k. Then, as n and m are scaled to infinity while fixing N, we can obtain Proof. Suppose β 1 , β 2 , β 3 , and β 4 are universal constants. Under our MORE algorithm, let r j be the expected number of delivery trips of drone j. Then the expected distance traveled by drone j can be approximately where u(c i , c 0 ) is the distance from customer c i to the depot c 0 , l j is the total distance of the customer sequence in the delivery trip of drone j. According to the Markov chain convergence theorem, r j is approximately equal to each other, we can get Therefore, the expected total delivery distance Z of our MORE algorithm is approximately Applying the lower bound from Lemma 1, we can obtain As m goes to inifinity, we have Applying the upper bound from Lemma 1, we can obtain = lim m→+∞ Qmr avg Nµ = Qr avg Nµ .
Algorithm 3 Algorithm for splitting customers into the individual delivery trip.

Input:
First served customer c s j , last served customer c e j , filled demand F j , drone capacity Q, observed daily demandD 1 , · · · ,D n . Ouput: A sequence of serving customers for the individual delivery trip. Initialize the dictionary d dic for all customers with zero; Break this loop; If the total demands of served customers achieve the drone capacity End the current delivery trip, start a new delivery trip; If the current customer cur needs to be served Serve the current customer cur with 1 demand; Set d dic [cur]+ = 1 Set i+ = 1 Output the individual delivery trip.

MORE with Time Window Constraint for Customers
In this section, we consider a special case of MORE problem (MORE-TW for short) that consider the time window constraint for each customer. We first formulate the MORE-TW problem. Then, we propose an algorithm for MORE-TW.

Problem Definition for MORE-TW
The special case of MORE problem, MORE-TW, is to determine the delivery trip for at most m drones towards minimizing the total expected delivery time of all the drones under the constraint that each customer should be served after the starting time of the time window. Formally, MORE-TW can be defined as follows.
where V is the speed of the drone, and a j is the starting time of the time window for customer c j .

Solution for MORE-TW
We illustrate our algorithm in Algorithm 4. First, we only consider that a drone starts from the depot c 0 , passes through all customers in customer set C, and finally returns to the depot c 0 with the goal of minimizing the total traveling time under the constraint that each customer should be served between its time window, which can be transformed into the vehicle routing problem with a time window. Therefore, we solve the classical vehicle routing problem with a time window to obtain a VRP tour. We pick an arbitrary customer and its orientation in the VRP tour as the starting customer and orientation. We relabel all customers according to their order in the VRP tour to obtain new customer sequences. Second, we divide the new customer sequence into m disjoint groups to obtain the primary customer set Φ 1 · · · Φ m . Third, we assign k-overlapped serving customers of drones to obtain the extended customer set Ψ 1 · · · Ψ m . Fourth, we use Algorithm 2 on the primary customer set and extended customer set to obtain the individual delivery trip for each drone.

Algorithm 4 General Algorithm for MORE-TW.
Input: Customer C, time window t i = [a i , b i ] for each customer c i , drone capacity Q, observed daily demandD 1 , · · · ,D n . Ouput: A sequence of serving customers for the individual delivery trip. Solve vehicle routing problem with time windows with the instance (C, t) as input to obtain new customer sequence; Divide new customer sequence into m disjoint groups to obtain primary customer set; Assign extended customers to obtain extended customer set; Call Algorithm 2 with the instance ( Φ 1 · · · Φ m , Ψ 1 · · · Ψ m , Q,D 1 · · ·D n ) as input to obtain the delivery trip; Output the individual delivery trip.

Simulation Results
In this section, we conduct a simulation experiment to verify the performance of our proposed MORE algorithm in therm of drone capacity Q, number of customers n, number of overlapped served customers k, and number of drones m.

Simulation Setup
In our simulation, unless otherwise stated, we set Q = 100, n = 200, m = 20, k = 10, D j ∼ U(0, 8), (j = 1, · · · , i = n) and our area size is 100 km × 100 km. We consider the overall traveling cost as a metric for the evaluation of MORE. Moreover, each experimental result is obtained by averaging the results for 100 topologies of customers.

Baseline Setup
We develop three other algorithms for comparison with MORE, i.e., Selected Customers Dedicated Trips (SCDT), where the primary customer sets are determined via our proposed algorithm for the TSP problem, the extended customer sets are the same as the primary customer sets; Grouping Customers Overlapped Trips (GCOT), where the primary customer sets are determined by grouping the customers of the neighborhood, the extended customer sets are obtained using our extended customer assignment method; and Grouping Customers Dedicated Trips (GCDT). We set n = 20, m = 4, k = 5 and Q = 20. We can obtain the delivery paths of our proposed MORE algorithm and SCDT algorithm as shown in Figures 8 and 9. As shown in Figure 8, there are three drones to be set out to meet all the demands of the customers while there are three drones to be set out as shown in Figure 9.   Figure 10 shows that the overall traveling cost for four algorithms first decreases with drone capacity Q when Q is small because the larger the drone capacity, the more customers can be met in one delivery trip by drone. Then, the overall traveling cost keeps relatively stable when Q exceeds 90 because the capacity is large enough for customers to meet almost all its assigned customers.

Impact of Number of Customers n
Our simulation results show that, on average, MORE outperforms SCDT, GCOT, and GCDT by 31.15%, 68.55%, and 72.51%, respectively, in terms of n. Figure 11 shows that the overall traveling cost for four algorithms grows almost linearly as the number of customers increases from 20 to 320. This is because when the customers are uniformly distributed, the delivery trip patterns for these four algorithms will not change too much as the number of customers increases. Our simulation results show that, on average, MORE outperforms SCDT, GCOT, and GCDT by 23.37%, 76.70%, and 77.29%, respectively, in terms of k. Figure 12 shows the overall traveling cost for four algorithms almost stays stable with k. This is because the algorithm of SCDT and GCDT is not affected by k. Our algorithm first keeps stable when k does not exceed 10, and reaches the minimum cost when the number of overlapped served customers is 10. While k exceeds 10, Our algorithm keeps stable, which states that the overall cost is minimum when we select k as the number of primary serving customers. Our simulation results show that, on average, MORE outperforms SCDT, GCOT, and GCDT by 19.02%, 78.90%, and 79.60%, respectively, in terms of m. Figure 13 shows the overall traveling cost for four algorithms stays almost stable with m. Because when the number of drones is small, the drone can fly multiple delivery trips to meet the demands of customers.

Limitation
The key novelty of this paper is taking stochastic demand into consideration. However, we ignore the impact of weight on the maximum speed, maximum flight time, and maximum distance of drones. In the future, we will consider more characteristics of drones in designing our delivery routing. The main idea of our proposed algorithm for stochastic demands can be used in designing the new algorithm by taking more characteristics of drones into consideration. For example, the first three steps considering stochastic demands can be used due to the robust consideration for demands.

Conclusions
The key novelty of this paper is first taking into consideration the stochastic demand for customers and capacity constraints for each drone for designing drone delivery trips.
The key contributions of this paper are formulating the problem of determining drone delivery trips with stochastic demands for customers, proposing an effective algorithm, and conducting simulations to verify our theoretical findings. The key technical depth of this paper is in converting the original NP-hard problem to the realm of the classical vehicle routing problem, decoupling the interrelated relationship of stochastic demands and capacity constraint. Our simulation experimental result shows that our proposed algorithm achieves good performance and can outperform the other comparison algorithm by 79.60 percent.

Data Availability Statement:
The data presented in this study are available on request from the first author.