Modeling optimal drone fleet size considering stochastic demand

Hadas, Yuval; Figliozzi, Miguel A. Article Modeling optimal drone fleet size considering stochastic demand EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Hadas, Yuval; Figliozzi, Miguel A. (2024) : Modeling optimal drone fleet size considering stochastic demand, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 13, Iss. 1, pp. 1-12, https://doi.org/10.1016/j.ejtl.2024.100127 This Version is available at: https://hdl.handle.net/10419/325205 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc-nd/4.0/ EURO Journal on Transportation and Logistics 13 (2024) 100127 Available online 24 March 2024 2192-4376/© 2024 Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Modeling optimal drone fleet size considering stochastic demand Yuval Hadas a , Miguel A. Figliozzi b , * a Department of Management, Bar-Ilan University, Ramat Gan, 5290002, Israel b Transportation Technology and People (TTP) Lab, Department of Civil and Environmental Engineering, Portland State University, Portland, OR, 97201, USA ARTICLE INFO Keywords: Drone Fleet sizing Stochastic demand Optimization ABSTRACT The last mile delivery is particularly challenging for stochastic deliveries with narrow time windows. This topic is timely due to the rise of e-commerce and courier type services and the impacts of fleet size and vehicle type on delivery costs. A novel contribution of this research is to provide an optimization approach, extending the newsvendor model, to provide an optimal drone fleet sizing solution with stochastic demand in terms of number of deliveries and deliveries weight or payload from one central depot. The solutions obtained are robust, as shown in a comprehensive sensitivity analysis. 1. Introduction In recent years, there have been a lot of exciting announcements and pilot studies of drone, or Unmanned Aerial Vehicles (UAVs), deployments in freight transportation and logistics. UAVs have been featured frequently in the media following announcements made by large corporations such as Amazon (Vincent and Gartenberg, 2019). In the US e-commerce grew at a 30% rate in 2020 (eMarketer, 2020) and drone deliveries are expected to become a 7 billion US dollar market by 2027 (Insights, 2020). Drones are increasingly being utilized to deliver medical supplies and in courier type services. The COVID-19 pandemic has accelerated this trend. Drones arrive quickly by taking more direct paths and avoiding ground-based obstructions or congested roads. Most of the drone literature has focused on routing, path optimization, or scheduling of truck-drone teams. This research studies a novel drone fleet sizing problem with stochastic demand in terms of demand volume and size of the delivery. Unlike most studies, in this research, the objective function is profit maximization when there are costs associated with the aircraft size and unmet customer demands. This paper is organized as follows: a literature review is presented in the next section, followed by a formulation of the drone fleet size optimization problem with a stochastic number of demands and demand weights or payload. Payload is important because for drones, unlike ground vehicles, payload is a major constraint that severely reduces drone range and impacts on its costs. A comprehensive simulation and sensitivity analysis of a fleet sizing scenario is later presented. The paper ends with conclusions. 2. Literature review This section provides a brief background regarding drone-truck applications, fleet sizing, and drone emissions. The idea of utilizing both UAV and trucks (Murray and Chu, 2015) to improve overall delivery efficiency has been analyzed by many authors focusing on the actual design of routes and logistics systems in the last five years. There has been an explosion in the number of papers related to drone routing optimization. Several reviews present an overview of modeling efforts. Since the focus of this research is on drone fleet sizing, not on routing or scheduling, the reader is referred to the following surveys for drone applications, vehicle routing, and scheduling. A survey of applications of drones in civil problems is presented by Otto et al. (2018). A more recent survey by Macrina et al. (2020) focuses mostly on routing problems with drones and the review by of Chung et al. (2020) on drone-truck combined operations. In the drone literature, it is possible to find drone a small number of research efforts and models that incorporate stochastic elements. For example, Baloch and Gzara (2020) study parcel deliveries with drones under competition using multinomial logit models. In this work, there is competition and uncertainty regarding drone market share but also regarding overall market size. The work of Chen et al. (2021) use a Markov decision process to develop closed analytical solutions that are useful to get insights regarding drone fee structures and delivery capacity for drone delivery operations with random demands, different product categories, and multiple service zones. Drone fleet sizing has received scant attention in relation to drone-truck routing. The work of Lee (2017) utilizes modularity and simulation in drone design to * Corresponding author. E-mail addresses: [email protected] (Y. Hadas), [email protected] (M.A. Figliozzi). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.sciencedirect.com/journal/euro-journal-on-transportation-and-logistics https://doi.org/10.1016/j.ejtl.2024.100127 Received 13 July 2022; Received in revised form 18 August 2023; Accepted 18 January 2024 EURO Journal on Transportation and Logistics 13 (2024) 100127 2 estimate fleet size and optimize their operation. Rabta et al. (2018) studies a drone fleet model for last mile deliveries in humanitarian logistics using a MIP model where the objective minimizes the total traveling distance or cost taking into account drone payload and energy constraints and the installation of recharging stations. Chauhan et al. (2019) present a MIP formulation and heuristics to maximize drone fleet coverage, also taking into account drone payload and energy constraints and providing facility location, allocating drones to facilities and drones to customers. This work is later extended to provide robust solutions when accounting for uncertainty in initial energy battery availability and energy consumption (Chauhan et al., 2020). Torabbeigi et al. (2020) formulate a set covering problem and a MIP model, considering payload and battery consumption, to locate facilities and minimize the number of drones in a parcel delivery system. Point-to-point operations that extend the range of drone have been also studied (Pinto and Lagorio, 2022). Fleet sizing under demand uncertainty has been studied for specific time-sensitive applications. For example, Boutilier et al. (2017) takes into account spatial demand uncertainty and utilizes optimization and queuing to determine drone fleet size for the emergency delivery of automated external defibrillators. This work was extended by Glick et al. (2021) to include not only demand uncertainty, but also the impact of weather conditions on drone fleet sizing. An integrated location-queuing model was later developed using realistic response times and analyzing the impact of congestion (Boutilier and Chan, 2022). Other practical applications include work on cold chain network design for vaccine distribution with time limits (Enayati et al., 2023). The models and results presented in this research are novel in several ways: (a) there is uncertainty regarding the number of requests per period, (b) there is uncertainty in terms of demand characteristics such as the parcel weight, and (c) it is a profit maximization problem where companies or drone operators face tradeoffs in terms of fleet size, type of drone, revenue, operating costs, and lost sales. The profit maximization approach is more appropriate when private operators are assumed and given the revenue/cost tradeoffs not all customers could be served. To the best of the authors’ knowledge, there is no similar published research related to drone fleet size optimization. 3. Drone fleet size problem formulation The formulation of the drone fleet size model utilizes ideas related to the newsvendor model for fleet sizing. The newsvendor model is a basic problem in stochastic inventory control and from the 1950s has been widely studied and applied to operations research and supply chain management problems (Wooldridge, 2015). The newsvendor approach has been used in the transportation literature, for example, to determine optimal fleet size in public transportation settings where agency costs are comprised by vehicle size, empty seats, and lost sales and user costs are related to waiting and overcrowding (Herbon and Hadas, 2015) or transportation adaptation of a supply chain model (Hadas and Shnaiderman, 2012). However, unlike the above-mentioned models, which are minimum cost models, the proposed model is based on profit maximization and applied to drones, which requires a different formulation incorporating payload constraints. 3.1. Assumptions The main assumptions in this research are the following. 1. Delivery times are short like in a typical courier service. Without loss of generality in this research, it is assumed that the planning period is ½ hour and that drone deliveries are made in this ½ hour period. Hence, for fleet size planning purposes, a ½ peak demand period is assumed to determine drone fleet size. 2. The distribution of demand is known, but the actual demand to be served in each period is only known at the start of the period. 3. One trip for one customer per drone can be completed in the ½ hour demand period. This is a reasonable assumption, taking into account flying time, takeoff and drop-off time, cargo preparation and battery swapping time. 4. Drone purchase cost is a function of payload capability. Drone purchase cost and operational cost is an increasing function of drone size. 5. Drone service area and range are constant and independent of drone size. Drone costs increase as a function of payload and reflect larger batteries to serve the same service area. 6. A customer can be served by a drone or a truck, but drones travel from the depot directly to the customer, i.e. the drone is not using the truck as a base for take-off or landing. The first two assumptions are similar to assumptions found in the traditional newsvendor problem, where the decision (fleet size in this case) must be made before actual demand realization is known. In particular, the ½ hour time window in the case of drone deliveries can be considered as planning for the peak hour demand where customers are willing to pay an extra fee or premium for a fast delivery service as well discussed in Baloch and Gzara (2020). The following three assumptions are needed for drones and not found in the newsvendor literature. As a result of these three assumptions, it is necessary to add a decision variable: drone payload (size) that greatly increases the complexity of the problem. These assumptions are realistic and not very restrictive. The main limitation of the newsvendor model is that demands cannot be carried over onto the next period, i.e. demands that are not served in a given period are lost. It is not possible to postpone servicing a demand or to allow for demand backlogs in inventory control terminology. Another limitation of the model is that in the real-world, fleet size is an integer variable and there is a finite set of drones (payloads) to choose from. The results of the model provide the best solution, assuming continuous variables, but because the objective function is fairly “flat” around the optimal, the rounded values provide a very good approximation. Alternatively, given that the number of potential payloads is limited, the payload value can be set and the model can be run optimizing only fleet size. The model is meant to provide managerial insights to guide drone fleet sizing decisions. In practice, there could be many other application specific constraints, like financing and cash flows that are outside of the scope of this paper. The model considers two decision variables, the fleet size (N) and drone’s capacity or payload (V). However, to simplify the presentation of the model, and without loss of generality, N is also referred to as the maximum number of planned trips per period because it is assumed that each drone can only serve a customer per period. Having N as the number of planned trips, can further relax the non-integer results, which is acceptable as N is the average number of trips per time window, similar to the way public transport headways are being used (Hadas and Shnaiderman, 2012). 3.2. Random variables and distributions The decision variables are related to two distributions or random variables: the demand distribution or number of trips per period and the distribution of customer payloads. In the proposed model, the two random variables should follow three requirements: 1) non-negative intervals of finite length, 2) twice differentiable, and 3) fit or estimate common distributions. The first requirement stems from the nature of the decision variables, i.e. both are non-negative, with the payload bounded by minimal and maximal size related to economical and physical properties, while the demand can be estimated within a certain confidence interval. The second requirement reflects the need for a closed-form derivation of the optimal solution. The last requirement satisfies the need to estimate different distributions representing demand and weight patterns. The four parameter Beta distribution (McDonald and Xu, 1995), with Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 3 the shape parameters α and β in the positive interval [l,u]satisfies all three requirements. It is bounded by parameters l and u and differentiable per equation (1). f(x) = (x−l) α −1(u−x)β−1 (u−x) α +β−1Γ( α )Γ(β) Γ( α +β) (1) Given that the Gamma function is defined as Γ(n) = (n−1)!equation (1) is also polynomial, which will be useful later for obtaining the global optimal solution. The shape parameters α and β can be modified to define both symmetric ( α =β=1)and unimodal skewed ( α >1,β>1, α ∕= β)distributions, with relatively straightforward calculation of the first two moments (equation (2) and equation (3)). E(x) = l+ α α +β(u−l)(2) Var(x) = α β ( α +β)2( α +β+1)(u−l)2(3) Furthermore, the most common distributions (normal, logistic, etc.) can be approximated using a Beta distribution. Furthermore, both uniform and exponential distributions are special cases when the parameters are ( α =β=1)and lim n→∞nBeta(1,n) = Exp(1)respectively. Fig. 1 provides four examples of the four parameters Beta distribution for different parameters α , β , l and u, in this research the demand attributes are number of customers per period and payload per customer (both stochastic variables). 3.3. Model formulation The model maximizes a profit function which is composed of four components: revenues (Tr), delivery costs (TCe), operation costs (TCo), and lost sales when customer requests cannot be satisfied due to fleet size and payload limits (TCl). Unlike the traditional newsvendor model, which only has one decision variable, the optimization model utilized in this research has two decision variables. In the traditional newsvendor model, there is a onedimensional probability space with two regions, shortage and surplus. In the proposed model, there is a joint probability space with two dimensions that are independent of each other based on the following properties. 1) The demand and weight distributions are independent. The former defines the number of deliveries, while the latter denotes the individual weight of each delivery. 2) The fleet size and payload decision variables are each uniquely associated with the demand and weight, respectively. 3) Following 1) and by definition, the decision variables are independent. Given the complexity of the proposed newsvendor model with two independent variables, a case where demand number and payload are positively or negatively correlated is left as a future research effort. Based on the previous three assumptions the two-dimensional probability space can be divided into four regions, as illustrated in Fig. 2. The bounds of the demand interval are denoted [ld,ud]and the bounds of the payload interval are denoted [lw,uw]. The two axes, x and y, correspond to the payload and fleet size, respectively. Each axis is bounded by the variable’s minimum and maximum values. The vertical and horizontal lines are projecting a given values of the decision variables. Region 1 refers to the so-called surplus region, as these actual demand and weight can be satisfied. Region 4 reflects the shortage related to both payload and fleet size. Region 2 refers to the situation in which unsatisfied demand is related to payload only (payload shortage). Region 3 refers to the situation in which unsatisfied demand is related to insufficient fleet size (fleet size shortage). Unlike previous regions, region 3 can be divided into two sub-regions, 3a and 3b, in which the latter can be satisfied by unused trips associated with region 2. Hence, the formulation of the problem should include the not a priori obvious case where the demand per period exceeds the number of drones available but some of the shortage is captured already in region 2 and should not be double counted. Furthermore, each region is associated with one or more of the revenue or cost components. The realized revenue and delivery costs are associated with region 1, as they are a function of the actual demand bounded by the selected fleet size and payload. Lost sales due to payload limit occur for all requests with weight over the payload limit (region 2). Lost sales due to fleet size limit occur for fleet size limits (region 3), and lost sales due to both variables are associated with region 4. The exact formulation of each component is herein defined, with fd(x)and fw(y)as the demand and weight probability distribution functions, respectively. Since both probability distributions are independent the joint probability fdw(x,y)can be written as fd(x)fw(y). Following that, the integrals over fd(x)and fw(y)can be replaced by the probability functions: Fig. 1. Examples for the four parameters Beta distribution. Fig. 2. probability space’s zones. Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 4 Fd(V) = ∫ N ld fd(x)dx (4) Fw(V) = ∫ V lw fw(y)dy (5) In turn, equations (6)–(9) define the aforementioned regions’ expectancies. Ed(N,V) = ∫ ud ld Min(x,N)fd(x)dx Fw(V) = ⎛ ⎝∫ N ld x fd(x)dx +N(1−Fd(N))⎞ ⎠ Fw(V) (6) Ew(N,V) = ∫ ud ld w(N,V,x)fd(x)dx =∫ ud ld Min(x,N)fd(x)dx (1−Fw(V)) =⎛ ⎝∫ N ld x fd(x)dx +N(1−Fd(N))⎞ ⎠(1−Fw(V)) (7) Es(N,V) = ∫ ud N s(N,V,x)fd(x)dx =⎛ ⎝∫ ud N x fd(x)dx −N(1−Fd(N))⎞ ⎠Fw(V)(8) Esw(N,V) = ⎛ ⎝∫ ud N x fd(x)dx −N(1−Fd(N))⎞ ⎠(1−Fw(V)) (9) The additional deliveries based on the unused trips resulting from the payload limits can be formulated as the expectancy of the minimum between regions 2 and 3. Equation (10) defines the additional successful deliveries resulted with the unused trips due to payload limits. i.e., as each planned trip has a probability Fw(V)of not being served it can potentially satisfy part of the unmet deliveries. Ed ′ (N,V) = ∫ ud N Min (s(N,V,x),w(N,V,x))dx =∫ ud N Min ((x−N)Fw(V),N(1−Fw(V)))fd(x)dx (10) Equation (11) defines the total revenues’ expectancy, where R is the revenue or price paid per delivery. Equation (12) and equation (13) define the operation costs’ expectancies related to the fleet size and drone size (maximal payload). The former are the costs incurred by the actual deliveries, with the coefficient Ce, which is the average cost to complete a trip related to energy and battery consumption. The latter is the fixed cost required to secure the availability of the drones regardless of the actual deliveries (purchasing, labor, facility). Two coefficients are associated with the fixed costs, Cv and Cf, the former with the drone size and the latter with the fleet size. Finally, equations (14)–(16) define the penalties’ expectancies associated with unserved customers due to fleet size and payload limits. Both penalty components are associated with the same penalty per unserved customer Cl. Tr(N,V) = R(Ed(N,V))+Ed ′ (N,V)(11) TCe(N,V) = CeV(Ed(N,V) + Ed ′ (N,V)) (12) TCo(N,V) = N(Cf+CvV)(13) TCs(N,V) = ClEs(N,V)(14) TCw(N,V) = ClEw(N,V)(15) TCsw(N,V) = ClEsw(N,V)(16) The objective function is then to maximize the profits (for clarity, the decision variables were removed from the functions). max Tp =Tr −TCo −TCe −TCs −TCw −TCsw ld ≤N≤ud lw ≤V≤uw However, part of the penalties should be eliminated. This fact follows from equation (10) where additional deliveries are carried out. Hence, TC sw, which is the absolute difference between w(N,V,x)and s(N,V,x), replaces TCs and TCw as follows: max Tp =Tr −TCo −TCe −TC sw −TCsw. ld ≤N≤ud lw ≤V≤uw (17) where Esw (N,V) = ∫ N ld x(1−Fw(V)) fd(x)dx +∫ ud N |N(1−Fw(V)) − (x−N)Fw(V)| fd(x)dx (18) TC sw (N,V) = ClE sw (N,V)(19) The model can be rewritten as follows. Equation (10) can be simplified given that regions 3a and 3b are divided at the line defined by equation (20). N ′ (V) = N Fw(V)(20) as a result, equation (10) can be rewritten as: Ed ′ (N,V) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∫ N Fw(V) N (x−N)fd(x)dx ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ Fw(V) + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∫ N Fw(V) ud N fd(x)dx ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (1−Fw(V)) (21) With N ′ (V) ≤ ud, eliminating the use of the minimum operator (equation (10)) by dividing the integral from N to ud into two separate integrals from N to N Fw(V)and from N Fw(V)to ud. For clarity, the full development of equation (21) is provided in Appendix A. Based on equation (20), equation (18) can be rewritten as: Esw (N,V) = ∫ N ld x(1−Fw(V)) fd(x)dx +∫ N Fw(V) N (N−x−Fw(V))fd(x)dx +∫ N Fw(V) ud (x Fw(V) − N)fd(x)dx (22) Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 5 Once again, for clarity, the full development of equation (22) is provided in Appendix A. As a result, the revised objective function is max Tp =Tr −TCo −TCe −TC sw −TCsw. ld ≤N≤ud lw ≤V≤uw s.t.N ′ (V) ≤ ud (23) 4. Global optimal search algorithm The objective function (23) is non-convex, as depicted in Fig. 3, hence traditional optimization algorithms does not guarantee finding the global optima. However, the model has several properties that enable the construction of an efficient global optimal search algorithm. The objective function components are twice continuously differentiable functions, as such, existing algorithms that use convex relaxation of non-linear twice continuously differentiable functions can obtain the global optima. Hence the α BB algorithm (Androulakis et al., 1995) or the QBB algorithm (Zhu and Kuno, 2005) are possible algorithms that can be used. Furthermore, as fd(x)and fw(y)are polynomial functions it follows that Tp is a polynomial function as well, as Tp’s components are integrals over polynomial functions. This in turns, simplified the search for the global optima, even though Tp is non-convex. Specifically, the search procedure can be constructed as follows. 4.1. Step A: model relaxation 1) Objective function (23) is solved unconstrained. 2) Find all roots of the unconstrained function (23): ∂ Tp ∂ N=0, ∂ Tp ∂ V=0(24) 3) Iterate over all roots and identify the global optima. 4) If the global optimum is feasible, given the constraint, stop. Otherwise continue to step B. Solving (24) is a relatively simple task, as the problem of finding all roots of a multivariate polynomial is a well-studied problem (Geil, 2015), specifically when solving over a bounded interval, such as the problem at hand. Then, commercial packages can be used that implement traditional search algorithms such as the Newton-Raphson. 4.2. Step B: Lagrange multiplier If the global optimum is outside the feasible region, then solve the objective function with equality constraint, as the solution resides along the constraint (equation (25)). max Tp ′ =Tr −TCo −TCe −TC sw −TCsw −λ(N ′ (V) − ud) ld ≤N≤ud lw ≤V≤uw (25) 1) Find all roots of equation (25) based the Lagrange multipliers method. ∂ Tp ′ ∂ N=0, ∂ Tp ′ ∂ V=0, ∂ Tp ′ ∂ λ=0(26) 2) Iterate over all roots and identify the global optimum. 5. Drone fleet sizing results The implementation of the formulation and search algorithm is applied in this section to a base case study. To further study the stability of the solution, a sensitivity analysis and noise in the parameters are also studied in this section. 6. Case study The model was implemented with Maple 2020 package (Maplesoft, 2020). Solutions were obtained instantly, as the model is differentiable with only range constraints for the decision variables. Assuming a ½ hour period, with 1 trip per period per drone, the following parameters were selected for a case study. The revenue (R)per delivery is based on rush courier delivery rates (Breakaway, 2021) after considering the service area of the drone and the number of deliveries per hour. For the drone costs, two key assumptions are labor costs per trip and the useful life of the aircraft (Figliozzi, 2018) and a base cost for unmet delivery is the cost of an alternative ground delivery like Uber Eats (Gridwise, 2020). The appeals of drones are both potential delivery time savings and cost (Vincent and Gartenberg, 2019). It is important to note that the profit for a drone delivery assumes a premium due to fast and reliable delivery and the cost of unmet demand is related to using a ground service to cover the unmet drone demand and/or the loss of a customer to a rival service provider. Drone services are still in its infancy and it is not trivial to estimate the value of lost customer (Hogan et al., 2003), due to the difficulties to identify precise or narrow ranges for most of the parameters an extensive sensitivity analysis exercise is performed in the following section. R =$12.5/delivery Cl =$ 5.0/unmet delivery Cf =$ 1.5/delivery Ce =$ 0.2/delivery-kg payload Cv =$ 0.1/delivery-kg payload For this study, we assume that a multi-copter can carry up to 2.5 kg of payload with an effective range of 10 km. As a reference, Amazon is studying drones to deliver up to 5 pounds (2.27 kg) in 1/2 h or less, which roughly agrees with the assumptions made in this paper (Amazon, 2021). In additions, up to 5 pounds accounts for 80–90% of Amazon packages (Manjoo, 2016). As a reference, an effective range of 5 km covers downtown Manhattan in New York City. For service times, it is assumed that an average drone delivery requires 30 min. Note that for a total round trip of 10 km, flying at an average cruise speed of 18 m s the flying time is approximately 9–10 min, and the rest of time is consumed by other activities which include: take-off preparation, delivery time, Fig. 3. Objective function (Tp) as a function of payload V and fleet size N. Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 6 landing, battery swap, and cargo loading and securing (Glick et al., 2021). The values of Ce and Cv are determined taking into account that for delivery drones the payload is approximately 1/5 of the drone total takeoff weight and that drone energy consumption and costs are a function of drone size (Figliozzi, 2023a). It is implicitly assumed that many of the drone dimensions and cost data are approximately linear for small size drones with a tare (empty weight without payload and battery) of up to 10 kg as shown when analyzing real-world drone specifications (Figliozzi, 2018). The tare of the drone is assumed to be 8 kg and the battery 2.5 kg. The estimated cost related to energy and battery replacement are approximately $0.5 per delivery assuming: a payload of 2.5 kg, a battery energy density of 200 wh/kg, a cost of $275/Kwh, a battery life of 300 cycles, a 50% battery utilization to account for the safety factor and aging, a cost per Kwh of approximately $0.15/Kwh, and finally a drone energy consumption of 25 wh/km when flying. The estimated cost of drone as a function of size is estimated as $0.25 per delivery assuming: a payload of 2.5 kg, 3000 flight hours per drone, and an initial drone cost of $7500. A comprehensive sensitivity analysis is provided because drones and battery technologies are quickly evolving (Figliozzi, 2023b) and therefore it is recommended that assumptions regarding drone costs and capabilities in future research efforts should be adequate for the specific application and year of analysis. Demand and weight probability distributions parameters were set to fd(x) ∼ Beta( α =3,β=3,l=0,u=100)and fw(y) ∼ Beta( α =3,β=3, l=0,u=2.5)respectively. The weight distribution is in kilograms. With these parameters the expected demand and variance are 50 and 357 respectively and the expected weight and variance of 1.25 and 0.22 respectively. The optimal solution value is $458 obtained for N =75 drones and V =2.38 kg (5.25 lbs) of maximum payload per drone. The solution space is presented in Fig. 3, along with the break-even plane. The shape of each of the components and values associated with the optimal solution are presented in Appendix B. All the six figures are aligned with the axes of Fig. 3. 6.1. Sensitivity analysis The sensitivity analysis constitutes of the following assessments: 1) distributions’ shape analysis, 2) coefficients sensitivity, 3) noise level analysis and 4) robustness analysis. 6.1.1. Distribution shapes First, the effect of the demand and weight distribution shapes was investigated. For that, the base solution, which has symmetrical distributions was compared to positive and negative skewed distributions. The results are presented in Table 1. When compared to the symmetrical distribution, the positive and negative skewed distributions have a drastic effect on the fleet size and drone size. Positive skew distributions decrease the optimal fleet size and drone size (demand and weight averages are smaller), while negative skew distributions increase them. Moreover, as fd(x)and fw(y) are independent, each one affects its corresponding decision variable. The slight change of the other decision variable can be attributed to TCo (equation (13)), which has a multiplication of N and V. Please note that N was rounded, the slight change (0.1–0.2 for each column) is not presented. 6.1.2. Coefficients sensitivity Additional sensitivity analysis was performed for each of the coefficients with selected coefficient values, while holding all other coefficients to their initial values. This analysis examines the effect on N and V resulting from increasing or decreasing each coefficient. The results are summarized in Table 2. The results provide insight into the change of N and V with regards to the specific coefficient. The coefficient Table 1 Sensitivity analysis of the demand and weight distributions. Demand distribution Positive skew E=29, V=16 2 α =2,β=5 Symmetrical E=50, V=19 2 α =3,β=3 Negative Skew E=71, V=16 2 α =5,β=2 Weight distribution Positive skew E=0.71, V=0.4 2 α =2,β=5 N 51 75 91 V 1.90 1.93 1.94 Tp 243 465 703 Symmetric E=1.25, V=0.47 2 α =3,β=3 N 51 75 91 V 2.37 2.38 2.39 Tp 239 548 685 Negative Skew E¼1.79, V¼0.4 2 α =5,β=2 N 51 75 91 V 2.49 2.50 2.50 Tp 238 457 693 Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 7 that is being examined has a gray background. The objective is to verify that increasing or decreasing each coefficient is logical: 1) Increasing R (from 12.5 to 30, 100, and 550) increases both N and V, as higher revenue reduces the effect of the costs. The same applies to decreasing R (from 12.5 to 10, 5, and 1). 2) Similar behavior can be observed from Cl. Higher shortage costs must reduce then shortage, which is realized by increasing N and V. 3) Cv and Cf are the coefficients related to the so called fixed costs (equation (13)), so higher fixed costs (increasing the coefficients) decreases both N and V. However, as Cf is only related to N, the effect on V is minimal. 5) Ce is related to the operational costs (equation (12)), and have similar effect to Cv and Cf. However, Cv is multiplied by V (aside from the demand), hence it affects only V when the change is closer to the initial coefficient value. However, due to the power effect, the larger the coefficient, the higher the change of both N and V. The results verify that the behavior of all the coefficients are logical. 6.1.3. Noise level analysis Each of the objective function coefficients was injected with a random, uniformly distributed noise, in which c’ is a randomized coefficient c with a noise level in a given % range pr. The noise was injected for each coefficient before obtaining the optimal solution. A total of 1000 independent simulations were performed for each scenario with the optimal decision variables results recorded. c ′ =c+rand (1−pr…1 +pr (27) Fig. 4 illustrates the distribution of the decision variables optimal values in a range determined by pr = ±20% for each of the coefficients (Cl,Ce,Cf,Cv,R). The black dot represents the original optimal solution. All the coefficients have linear effect on N and V. Clearly, based on the slopes, Cf affect mostly N, while Ce,Cv has more effect on V. The parameters Cl,R have similar impact on N and V. In order to investigate extreme values, each of the coefficients was randomized within the following range: Cl−0.1000, Ce−0.5.1, Cf – 0.10, Cf – 0.5.1, and R – 1.500. Fig. 5 illustrates that for a wider range of noise levels, not all the coefficients are linear, especially the effect of R, Cf and Ce on N and V. Table 2 Sensitivity Analysis of the Model’s Coefficients. Fig. 4. Decision variables distribution by coefficient with pr ±20%. Fig. 5. Decision variables distribution by coefficient with extreme values. Y. Hadas and M.A. Figliozzi EURO Journal on Transportation and Logistics 13 (2024) 100127 8 6.1.4. Robustness analysis The robustness of the solution was investigated by introducing noise to all coefficients, excluding R. This is due to the uncertainty with drones’ operations, while R is associated with the delivery cost and is more predictable. Three scenarios were constructed for three noise levels: 5%, 10%, and 20%. The decision variables distribution is illustrated in Fig. 6, as a function of optimal fleet size (N) and payload (V). The black dot is the optimal solution without noise. The observations with less noise (5%) are close to the black dot and mostly covered by the observations with higher noise levels. To facilitate the observation of the distribution, there are three boxes bounded by the two-dimensional 5 and 95 percentiles of N and V and by noise level. It is evident that the solution is robust, for the noisiest level (20%) the solution range is +/−~1.5 drones and +/−~0.012 kg. when considering the 5 and 95 percentiles bounds. The impact of noise on profits is observed in Fig. 7 utilizing histograms. As expected, the histogram range is roughly proportional to the noise level. The black vertical bar is the optimal solution without noise. It is also observed that the solution is also robust in terms of profit with a range ±5% with the noisiest level (20%) with respect to the profit without noise. 7. Conclusions The last mile delivery is particularly challenging for stochastic deliveries with narrow time windows. Due to its characteristics, drones are suited for situations where fast and reliable deliveries are needed. This research developed a novel optimization approach for drone fleet sizing, extending the newsvendor model. The model provides an optimal drone fleet sizing solution with stochastic demand in terms of two decision variables: 1) number of deliveries and 2) deliveries weight or payload. Unlike other studies, in this research, the objective function is profit maximization and there are costs associated with the aircraft size and unmet customer demands. An efficient algorithm guarantees that the optimal solution is found, which can be used in large-scale delivery scenarios. The solutions obtained are robust, as shown in the comprehensive sensitivity analysis. The sensitivity analysis showed the higher importance of fleet size in relation to drone size, though as drone technology is rapidly evolving it is important to consider that this finding may change in the future or other specific drone applications. This research focused only on drone fleet sizing, however, having a fleet of drones and a fleet of trucks is a more resilient approach since drones may not be able to operate effectively with adverse weather conditions and similarly ground vehicles may be hindered sometimes by congestion or ground network disruptions. However, having two different vehicle types is also likely to increase costs not only in terms of labor but also in terms of facilities and capital costs, and modeling these costs and tradeoffs requires a major future research effort. Further research is also necessary to study both drone and truck fleet sizes when considering profit and/or sustainability goals. Other research direction can extend the model with additional decision variable, the range (or coverage). The larger the range, the higher the potential demand. However, the number trips per drone will decrease due to the larger round trip. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research was supported by a grant from the Freight Monitoring Research Institute (FMRI), a USDOT Transportation Center. Any errors or omissions are the sole responsibility of the authors. APPENDIX A In this appendix, equation (21) and equation (22) are fully developed and proofed. For both equations two propositions are introduced, followed by each equation development. Fig. 6. Decision variables for 5%, 10%, and 20% noise levels. Fig. 7. Profit distribution for 5%, 10%, and 20% noise levels. Y. Hadas and M.A. Figliozzi