Planning delivery-by-drone micro-fulfilment centres

Delivery drones are a disruptive technology that is spurring logistics system change, such as the adoption of urban micro-fulfilment centres (MFCs). In this paper, we develop and implement a two-stage continuum approximation (CA) model of this disruptive system in a geographic information system. The model includes common CA techniques at a local level to minimise cost, and then these local solutions are used in a second stage regional location-allocation multiple knapsack problem. We then compare the drone MFC system to a traditional delivery-by-van system and investigate potential cost or emissions savings by adjusting time-window demand, logistical sprawl, electric van alternatives, and MFC emissions. Furthermore, we conduct a sensitivity analysis to show that uncertainty in demand and effective storage density both significantly influence the number of MFCs selected and benchmark our model against commercial solvers. This methodology may also be further developed and applied to other new delivery vehicle modes.


Introduction
Urban logistics, supply chains, and freight transportation worldwide are being disrupted by innovative technology and changing customer expectations.One such technology is delivery drones.Researchers estimate that current drone technology is advanced enough that drone delivery systems could economically service about 30% of the world's urbanised population (Aurambout, Gkoumas, and Ciuffo 2019), and which investment industries estimate is valued in 2020 between one and two billion USD globally.In North America, at least four large retail firms that have e-commerce presence, Amazon, Walmart, London Drugs, and CVS, are pursuing drone delivery as a delivery option for their customers.In this introduction, we first cover the state of the industry, then we will outline the objective of this paper and our contributions.

State of the industry
Over the last several years, many online retail companies have stated their interest in drone delivery.Amazon was an early proponent, which, according to their own promotional material online, has been designing and developing their own drone and delivery system to complement their already growing logistic service capabilities since 2013.Furthermore, Alphabet (the parent company of Google) has wide news coverage of their trial drone delivery service, called Wing, which has been in operation in the USA, Australia, and Finland, since 2019.In 2020, a pharmacy chain (CVS) paired with a major logistics service provider (UPS) to deliver prescription medicine in Florida via drone delivery.Other major supermarket chains such as Walmart and London Drugs, have declared similar business pairings with new drone delivery service companies, such as Flytrex and Indro Robotics, respectively.This phenomenon is occurring worldwide as, for example, JD.com showcases its drone delivery in China since 2020, and Swiss Post is currently delivering mail from three locations in Switzerland.
In parallel with this drone delivery mode adoption is the rise of e-commerce; because of e-commerce, planners and managers can no longer assume only a traditional supply chain from factory to warehouse to store, but rather a fast then direct, warehouse-tohome delivery option is often now expected in developed economies (Perboli et al. 2021).Moreover, the seller increasingly performs delivery themselves, and not a third-party logistics supplier, as the seller seeks greater control over their customer experience and cost savings.
Drones and other autonomous delivery vehicles are attractive to these Seller and Service Providers (SSPs) due to the public's perceived eco-friendliness of drones (battery-electric powered), and due to drone's ability to deliver in short time windows.
For instance, according to the Canada Post E-commerce Survey in 2020, consumers are increasingly influenced and aware of the environmental impacts of their shopping behaviour; reportedly, over half of Canadians are frustrated with excessive packaging, and over 40% of consumers report shopping with retailers who support an eco-friendly agenda.Drones may have further environmental sustainability benefits when they are paired with advances in micro-fulfilment centre (MFC) warehousing technology that promises to address both the transportation emissions issue of urban freight and reduce packaging requirements.In a TED talk, the CEO of Attabotics, a company that designs and manufactures MFC inventory storage systems, described how MFC technology can reduce packaging.Given that drones carry individual packages, then these items may require less packaging protection than if they were packed tightly among many other parcels in a traditional delivery van.In addition, due to the COVID-19 pandemic that struck the world in 2019, contactless drop-offs are increasingly desirable, making robotic drone delivery more attractive than traditional human delivery.
SSPs, such as Amazon, Walmart, and JD, enter the logistics market with different motivations than traditional Logistics Service Providers (LSPs) which make the above benefits of drones even more appealing.SSPs are committed to less logistics-only infrastructure than LSPs and so these new entrants to the market are more interested and able to adopt emerging technologies, such as drones.This new perspective changes costs and supply chain structures.Therefore, new models are required to understand this changing business environment and to incorporate drone delivery.

Study objectives and contributions
Our objective is to develop a model that can solve a multi-facility location-allocation with inventory problem.We desire that the tool can be easily understood by local stakeholders and used as a planning and discussion tool by local area planners and elected decision makers.Furthermore, we wish the model to be easily extended to incorporate other autonomous delivery modes, and other research questions, by future work.
Assuming the proposed benefits of a delivery-by-drone system appeal to an SSP, we consider how the SSP may implement such a system to serve a set of urban residential communities that regularly demand a variety of drone-transportable goods.In our scenario, the SSP wishes to deliver both standard and expedited packages.Expedited packages could be clothing, footwear, grocery items, electronics, books, health products, office supplies, jewellery, or other items that the SSP decides.The SSP controls which products are offered for this expedited delivery; Chen, Hu, and Solak (2021) given heuristics for choosing which products may be profitable for expedited delivery.We do not consider emergency medical items, such as defibrillators or blood.We assume that the last-mile delivery is performed by the item seller who also operates the MFCs.Previously, delivery was dominated light trucks (commercial vans) and human delivery drivers but, given advancements in robotics, we consider a near-future scenario where it is possible, safe, and legal for these items to be delivered autonomously via drones.
Our contributions are: 1.No previous study has considered an integrated multi-facility cost-optimised locationallocation for an urban delivery-by-drone system.2. We expand the field of continuum approximation (CA) modelling.The estimated costs of the two-stage method used to construct the location-allocation solution are for the first time compared to a GIS estimate of the solution's costs; the comparison validates the two-stage method for this case study and shows the method performs similarly to conventional CA techniques.a.By incorporating CA techniques into an integrated strategic model, we allow for future adaptations and extensions of the work, such as in van-drone hybrids and sidewalk autonomous delivery robots.
In summation, we address research gaps, as identified by Boysen, Fedtke, and Schwerdfeger (2020), of adapting existing strategic decision-making models for use with emerging technologies.

Literature review
In the following thematic literature review, we first outline previous studies and literature reviews that investigated drone delivery systems, micro-fulfilment centres, and recent twostage CA modelling.We highlight gaps that our work addresses throughout the review and summarise these in a research gaps section.

Drone delivery service
The study of drone delivery has gained increasing attention in literature as the technology has shown promise in industry.Moshref-Javadi and Winkenbach (2021) provide a broad literature review on drone delivery models.They find that most studies use either algorithm or mathematical programming methods, which aim to minimise system travel distance or financial system cost.They conclude that there is a severe lack of integrated problem studies, where topics such as inventory management or costs are considered.The authors also identify that few multi-facility problems have been addressed for drone delivery and they recommend that these multiple depot problems should be investigated further as, they argue, they are more likely to be the logistical system that is implemented in the real world than single facility systems.
At an operational level, there are studies that investigate the optimisation of individual drone flight paths; Coutinho, Battarra, and Fliege (2018) provide a review of these studies.Moreover, hybrid truck and drone systems (known as 'horsefly systems') at a tactical and operational level have been extensively studied, by authors such as Salama and Srinivas (2020).These, and other works, are covered by literature reviews of Rojas Viloria et al. (2020) and Macrina et al. (2020).Such horsefly systems are out of the scope of our current paper.Our paper covers an integrated multi-facility location-allocation problem for urban delivery-by-drone direct from MFCs.
At a strategic level, drones require infrastructure (e.g.launch locations) to support their effective deployment, and so many related facility location studies exist in literature.
Some studies consider delivery and response for disaster relief.Chowdhury et al. ( 2017) utilised a CA and GIS model to estimate optimal inventory levels and emergency response logistics locations for three counties in the southern United States.MacKle et al. ( 2020) investigated a similar medical response drone system, one specifically designed for longterm facilities that serve cardiac arrest patients across Northern Ireland.This study also used GIS, but it complemented GIS with a genetic algorithm to evaluate the financial costs and estimated lifesaving benefits of the proposed drone system.
Other studies have investigated location and drone network designs for non-emergency medical scenarios.Kim et al. (2017) considered the drone delivery of regular prescription medication to patients' homes in a rural area using a bi-level integer programming model.Arenzana, Javier, and Macias (2020) also used a programming method, but they looked at the inter-hospital case of blood delivery in the city of London, United Kingdom, where congestion makes drone delivery both a faster and cheaper solution than the current road ambulance inter-hospital delivery network.Both studies also included a fleet size estimate for the drones, which added to the modelling complexity and limited the problem size to under ten locations.
Chen, Hu, and Solak (2021) investigated cost-optimal drone delivery fleet size for only one MFC location using heuristic solving methods.These authors found that SSPs can profit from having a large fleet of drones to deliver high-value items with short wait times and that consolidating even just two packages per drone per trip can increase profits by over 50%.However, Moon, Salhi, and Feng (2020) investigated a similar location-routing problem for multi-compartment last-mile delivery which could be applied to a multi-capacity drone system; as their case study accounted for an added cost to operate the larger and more complex multi-compartment vehicles, they found that this added capacity was often not justified and led to increased system cost.Hong, Kuby, and Murray (2017) showed how the drone battery swap station location problem could be analysed in GIS and how this method allowed for modelling of obstacles and non-Euclidean pathing.Asadi and Pinkley (2021) conducted a detailed analysis of these replacement and recharging strategies at such battery stations to find that a policy of always swapping for fully charged batteries when replenishing the battery stations inventory was the best strategy, for a majority of cases.
Yet other studies have examined package delivery drone systems with various logistical structures.In a case study of Seattle, US, Shavarani et al. (2018) proposed a heuristic solution to determine an optimal arrangement of MFC locations and supporting drone recharge stations to completely replace Amazon's current van delivery system.Cokyasar et al. (2021) also developed and applied a heuristic programming method in a Chicago case study; the method located battery-swapping stations and allocated discrete demand points to the stations and allowed for a parallel truck and drone system, where some areas of the city were served by trucks and other areas were served by drones to achieve a low cost system -we also investigate such parallel systems in our case study.
In a statedly unique study, Baloch and Gzara (2020) investigated a drone delivery versus local grocery store delivery problem expressed as a multi-nominal logit market share model that they solve using mixed integer algorithmic solver that interchanges between solving a master problem (MFC location) and sub-problem (customer allocation).These authors find that, in their case study of New York City, that current regulations on the ratio of drone remote operators to drones controlled and the technological capability of drones to deliver to dense urban areas are barriers to drone delivery profitability.Baloch and Gzara also varied customer sensitivity to price, time, and an 'inherent attractiveness' of drones to find several instances where drone delivery may be profitable in a range or urban and sub-urban scenarios.Although these authors consider a profit-maximising model, their cost function only includes a piece-wise fixed cost per facility and a fixed cost per delivery whereas our work considers these costs in continuous space plus inventory costs.
Most studies only investigate financial costs as their objective metric, but some investigate carbon emissions.Goodchild and Toy (2018) used GIS to compare the emission intensity of truck and drone delivery systems in Seattle; they adjusted the energy efficiency of the drones in a sensitivity analysis to understand what drone characteristics were required if lower emissions were desired.Stolaroff et al. (2018) included emissions of both vehicles and warehouses in a comparison of drone delivery systems that used MFCs and systems that used traditional van and electric van delivery.Furthermore, Figliozzi (2017) investigated lifecycle emissions of replacing a van fleet with a drone system using only the CA method and determined that a drone system would result in fewer emissions.Figliozzi (2020) also compared van and drone delivery system emissions under different operational requirements of a system, such as logistical sprawl distance and time windows.Figliozzi (2020) showed that delivery time windows did not affect energy consumption per customer of a delivery drone system, but delivery time windows did affect vans, with shorter time windows decreasing van efficiency.These studies found several scenarios in which drones were the most emissions efficient option, but none found that drones were always best, thus showing that the strategic decisions made before the operational ones significantly affected the environmental impact of delivery drones.However, these authors did not consider the total financial costs of a system or inventory factors, so they could not state the cost changes incurred in the reduced emissions scenarios and thus could not state if profit-seeking companies would pursue these options without external incentive or regulation.
There are criticisms of using drones for delivery that may prove a barrier for deployment, such as privacy, security, safety, environmental, social, and employment implications.Interested readers are encouraged to see Chung, Sah, and Lee (2020) for further discussion on these topics.For our case study, we assume that these barriers are overcome.

Micro-fulfilment centres (MFC)
Urban logistics centres are a topic of growing interest in the urban freight literature and industry practice.Terminology is still not well defined nor agreed upon in the literature (Gunes and Goodchild 2021), but generally, and in this paper, the term 'micro-fulfilment centre' (MFC) is used to refer to locations that have deliberate short-to-medium term inventory storage and are owned and operated by one company.'Micro-hubs' may also be used, but this term may carry connotations of shared use.MFCs are like urban consolidation centres (UCCs), and the two types of centres share many transferrable insights, which have been explored in the literature.However, UCCs generally have no deliberate storage and are intended to reduce delivery vehicle volume into city centres, which can include municipally organised, owned, and even operated facilities, so work that examines inventory considerations is required to understand MFCs specifically (Urban Freight Lab 2020).
A primary benefit that is shared by MFCs, UCCs, and micro-hubs are the benefits as a distribution centre, where deliveries can be de-consolidated from large vehicles appropriate for inter-city travel to smaller vehicles more suitable for urban travel.These facilities are becoming more feasible given developments in innovative last-mile delivery business practises and technology.For example, Ballare and Lin (2020) show how a disruptive business practice, 'crowd-shipping', can reduce number of trucks, truck kilometres travelled, total operating costs, and total fuel consumption when compared to a classical hub-andspoke distribution system.The authors also highlight how locating the facilities is essential for the cost-competitiveness and sustainability of the system.Other authors, such as Sheth et al. (2019), show how an innovative vehicle types, such as electric assist cargo-bikes, can be used to economically replace even relatively lengthy urban delivery tours when paired with these urban distribution centres.
However, it is important that both business practises and environmental impacts of the innovative vehicles are examined.Lin, Chen, and Kawamura (2016) used CA methods to optimise the vehicle type, vehicle number, and route choices to minimise the cost of a UCC system in a Chicago case study then estimated the emissions output using a US Department of Energy model.The authors found that, counter to intentions, the expected outcome of some municipal policies would be cost savings for the system owner at the expense of increased vehicle kilometres travelled and increased emissions.It is clear from this study that well-intended local municipal policies should be modelled before implementation to ensure that the desired scenarios, leading to both system cost savings and emissions savings, are incentivised.
Most academic studies on urban delivery drones have so far not considered financial costs, including inventory costs, nor have they included the wider strategic decisionmaking implications of MFCs (deliberate multi-day storage).Rather, many studies focus on UCC models (simple cross-docking and temporary storage) and although some findings are transferrable to MFCs, some findings are not.For instance, Stolaroff et al. (2018) do not evaluate the financial cost of the studied urban drone delivery system, and so the cost of adopting such a system, and whether an SSP may be motivated to implement such a system, are unknown.Although Shavarani et al. (2018) consider the unit transportation cost and facility opening cost, they do not consider the variable sizing costs of each facility, inventory costs, or upstream logistic costs that an SSP must consider.Lemardelé et al. (2021) evaluate launching drones from UCCs, but they do not consider stored inventory nor variable sizing of facility costs.Figliozzi (2020) evaluates only emissions; furthermore, the author does not consider the impacts of multiple facilities in an urban area nor inventory at the facilities.There is a research gap in the literature, and therefore, there is a lack of understanding of the holistic costs of an urban drone delivery system using MFC locations as proposed by industry.

Continuum approximation location-allocation modelling
Turning now to CA facility location problems, Newell (1973) developed the convex optimisation method and demonstrated how to minimise analytically the combined transportation and warehouse set-up costs.Additional works over the decades are covered in literature reviews by Langevin, Mbaraga, and Campbell (1996) and Ansari et al. (2018).Notably, Erlenkotter (1989) accounted for economies of scale in the objective cost functions and introduced three new distance metrics that have various benefits.Moreover, Rutten, van Laarhoven, and Vos (2001) then added further cost terms, including inventory stock, and tailored the cost function for a specific case involving trucks in a Manhattan grid.Ultimately, these analytical approaches provide the first stage of understanding a locationallocation problem and give estimates of the number of facilities, size of facilities, and the magnitude of the catchment areas and are not intended as a final design; rather CA methods are used to narrow the field of possibilities and to quickly inform high-level decision making and public discussions.
The two-stage CA approach that built on this early CA work, which we extend, was developed by Tsao et al. (2012).The authors examined the whole of the US in a numerical example to demonstrate their approach.They divided the country into many small, equally sized areas that could each be assumed to have a uniform demand.Then, they applied the optimisation methods shown by Newell (1973) to determine optimal numbers of warehouses within these areas.In this instance, the sum of the facilities across these areas then formed the national solution, and they used the case study to investigate iterations of inventory considerations.Tsao et al. (2012) found an almost 12% reduction in total system cost using this two-stage divisive method compared with assigning the whole country a uniform average demand density.The authors, however, did not validate their final solution as they did not compare their estimated values for transportation costs in the optimisation process with the estimated transportation costs arising from their final solution.Furthermore, they did not further allocate demand points to the facilities; rather, they left this step for a future stage.We address both points in the current work.Chowdhury et al. (2017) used a similar two-stage method in a drone delivery system for a disaster relief study of the south coast of the US.The authors described a 'grid-couple-cover' approach that was implemented using a 'trial-and-error' method.This method was used to create a raster-like grid over the study area.There is not enough information on the trialand-error method as described to replicate the study.Furthermore, the described approach may lose data precision; since every grid square must be of the same geographical size, the data must be gathered for known grid sizes or interpolated for these exact sized areas, which limits the accurate use of past data.By allowing for arbitrarily sized and irregularly shaped areas to form the solution, the method we develop allows for existing data to be more easily used to design systems.Furthermore, similar to Tsao et al. (2012), Chowdhury et al. (2017) did not compare their estimates used to construct their location allocation solution with true costs likely resulting from implementing that solution, and thus, they did not validate their method.Moreover, neither set of authors compare their two-stage CA methods with a discrete method, such as that produced by an in-built solver to a GIS software, which we do.In summary, our review of the literature agrees with the results of Ansari et al. (2018), that there are few studies which utilise a CA method to model location-allocation with inventory problems in recent years, which we attribute to the success of discrete and heuristic methods.However, there has been little attention to the specific area of CA integrated models motivates us to investigate the approach as a research gap, as suggested by Ansari et al. (2018).

Research gaps
To address the research gaps identified in the literature review, we present a two-stage CA and algorithm allocation method, as summarised in Table 1.
We then apply this model to a case study to validate the method and explore deliveryby-drone MFCs as an urban logistics solution.We estimate emissions, costs, and inventory factors as previous studies have done and, by considering the strategic problem, are able to make novel insights into the issue.

Methodology
In this section, we formulate the problem, including notation and assumptions, and introduce the CA modelling method used to solve it.We present the objective cost function for the drone system and how it is optimised with respect to the number of MFCs, and finally we show the algorithm for allocating communities in the location-allocation second stage.

Problem formulation
In this section, we present the system boundaries and define the system variables.
Consider an SSP that desires to implement a multi-commodity delivery-by-drone system in a city.The SSP wishes to offer customers fast, same-day delivery for some items (expedited items) within a time-window that the SSP sets (such as 'two-hour delivery').Consider also that the local and national regulations allow these flight operations, and that cost minimisation while serving the entire city is the aim of the SSP; what logistical and engineering factors must they include when they determine their infrastructure network?
The SSP aims to offer their service to every community (subscript C) within the municipal region (subscript R) by using MFCs (subscript U).The SSP wishes to know the best number of MFCs (N), and the location-allocation arrangement of the MFCs to minimise their total annual operating costs (C Drone ).Each community area (A C ) within the regional area (A R ) must be allocated to an MFC and thus be within one MFC catchment area (A U ).We assume transhipments are not allowed between MFCs.
The city's adult population demands a typical number of packages per unit area per unit time, which is the mean regional demand density (μ R ).The SSP meets this demand with regional cycle stock (W α,R ) held across all MFCs in the region [subscript ordering is where variable subscripts are first, followed by applicable area subscripting second, e.g.W α,R is cycle stock (α) in the region (R)].Although the SSP has good demand prediction models, future observed demand is uncertain, and so the SSP also holds a regional safety stock (W β,R ) related to the standard deviation of regional demand density (σ R ).A portion of both stocks are stored locally at each MFC for use within each service area in an arborescent supply chain structure.Each MFC has capacity for the combined MFC cycle stock (W α,U ) and MFC safety stock (W β,U ).The former stock is held to meet catchment expected demand density (μ U ), and the later meets the variation in demand related to the catchment standard deviation in demand density (σ U ) between resupplies of the MFC.The sum of these two stocks prescribes the total capacity per MFC (W U ), which may be different between MFCs.The sum of total capacities across all MFCs in the region forms the regional capacity (W R ).
Considering inventory costs, the average wholesale value (u) of the various commodities per delivery to be stored in the MFC must be paid by the SSP before they receive deliveries to the MFCs; they receive these deliveries at a fixed resupply frequency (E) from a regional warehouse outside of the city boundaries.These resupply deliveries are conducted by an appropriately sized vehicle determined by the rate of resupply and the catchment cycle stock.The distance from this regional warehouse to a community is the community logistical sprawl (d c ).While in storage, packages can be stored at an effective storage density (m), which includes the need for walkways, sorting areas, item sizes, and the mix of standard and expedited packages; empirically, effective storage density is the observed number of packages sold in the resupply time frame divided by the total area of the storage facility.After being stored for a time, when a package is demanded, it is loaded onto a drone and flown to the customer at an average speed and following a path that can be approximated according to the product of a configuration factor (ϕ) and the Euclidean distance from the MFC.The drone then returns, having not travelled further than its maximum flight range (MFR).Only one customer is serviced per dispatch of a drone meaning that changing time-windows does not affect the total last-mile distance travelled given that we assume an adequate drone fleet size to always serve demand.
These outlined operations come at a cost.The cost of each MFC is separated into a fixed annual cost per facility (C f ) and a linear annual cost per square meter of storage space per facility (C s ); this cost must be either leased or its purchase financed at a yearly rate.To fill the MFC, the SSP must pay for each delivery (C d ) to each MFC, either paid to a third-party supplier or managed themselves.Typically, while packages are in storage, their wholesale cost has been paid by the SSP, but the packages have not yet been bought by a customer.This unmet cost is covered by the SSP, which, incurs an inventory holding cost rate (C h ).Other financing methods, such as a seller providing a platform on which others sell goods, and they only facilitate the transaction, (known as 'drop shipping') can avoid this cost, but we do not consider these options.Once bought, outbound delivery begins.The drone operation has costs of electricity, equipment, and remote pilot wages that can be expressed as a last mile cost per kilometre (C l ).
For the comparison to the traditional van delivery system, we consider the van travel cost (C t ), which includes the cost of fuel, capital expenditure, driver wages, maintenance, and all other operating expenses expressed as a cost per kilometre travelled.We assume that the drones only carry one package per trip whereas vans have a van package capacity (C), and the number of deliveries they can make per time-window is linearly affected by the number of time-windows per day (T) set by the SSP.
All above operations constitute the total operating cost of the urban consolidation centre delivery-by-drone system that the SSP wishes to consider.This total cost forms the objective function to be minimised.We assume revenue and demand are independent of MFC operations and structure, and so minimising cost is equal to maximising profit.The costs of this operation will be compared to the costs of a traditional delivery-by-van system, formulated following CA methods outlined by Daganzo (2005).We consider these costs for comparison as the significant distinguishing costs between the van and drone systems.We assume that costs further 'upstream' than the regional warehouse are similar between the systems.

Solution method
Our solution method is that we first solve CA cost optimisation model for each community in the region and determine the optimising design parameters for each community (N * C , W * C [we denote optimisers with an asterisk]).These costs and associated optimisers form input parameters for a second stage allocation where the communities are collected into service areas for the final location-allocation system design.This method leverages the benefits of CA models while using a framework to introduce an explicit spatial term that can implement the method and provide solutions in a discrete allocation environment, as recommended by Ansari et al. (2018).
The two-stage method is useful because it can estimate optimal solutions for areas of high spatially varying demand, areas not suitable for conventional CA methods, while also accounting for many costs not easily optimised holistically in discrete methods simultaneously or in commercial GIS software.The two-stage method has not yet been applied to an urban goods delivery problem.

Delivery-by-drone system objective function
Equation (1) shows the objective cost function that is to be optimised for each community in the first stage of the model.
The first term captures the cost of regular deliveries to the MFC; it is the number of resupply trips multiplied by the unit cost of these trips.The regional warehouse where these resupply trips originate from is the furthest 'up' the supply chain these costs will go; it is of interest to compare the drone delivery system to an alternative, traditional, delivery-by-van system, and distribution from the regional warehouse is the first substantial difference we model.The second term captures the expected transportation distance of the drones, the return trip from an MFC to and from an average delivery point within the MFC catchment area.From Newell (1973), we know that this distance increases with respect to the root of the catchment area.Increasing the number of MFCs in the community thus decreases the last-mile distance by the same relationship.Increasing community demand for an item (the product of community area and demand density) will linearly increase last-mile distance as too will changing the configuration factor (ϕ), which may be taken from Erlenkotter (1989) to represent a given distance metric and catchment shape.
The third term captures the sum of fixed costs across all MFCs in the community and is linear given the definition of fixed cost per facility.The fourth and fifth terms capture sizing costs; the former determines the storage space required for the cycle stock, and the latter determines the storage space for the safety stock.The sixth and seventh terms capture the holding costs of inventory, and, like storage costs, the former determines costs resulting from cycle stock, and the latter determines costs resulting from safety stock.For both storage space and holding costs, the cycle stock grows, as may be expected, with increases in demand and is not dependent on number of MFCs.The safety stock however does grow with the square root of the number of MFCs, due to a deconsolidation effect (Geoffrion 1979) which we show in Figure 1.Safety stock also increases, as may be expected, with increasing uncertainty in the community demand (σ C ) and sensitivity to stock-outs (β).
Once the optimisation is conducted, the resulting number of MFCs can be input back into this equation to determine the operating cost of the system in the community.See the supplementary information for further details regarding the terms for transportation and inventory costs.

Traditional Van system objective function
To compare the proposed drone system to a traditional delivery system, we use an estimate of transportation distance covered by a van system.Van delivery distance is estimated following a method by Daganzo (2005) and adjusted to include the number of time windows per day (T) in which the SSP offers expedited packages to be delivered.We assume that the vans deliver proportionally fewer packages in a shorter time window; for instance, if a ten-hour day is partitioned into two-time windows, then half of the packages delivered in Other estimates of touring distance are available, such as Figliozzi ( 2008), but we do not consider these because the parameters that require sample data for regression were not available for the study area.Similarly, we only examine financial cost to the SSP and do not consider societal costs of delivery van vehicle class (Holguín-Veras, Torres Cruz, and Jeff Ban 2013).
Equation (3) shows the objective cost function used to compare a traditional delivery by van and a drone delivery system.The first term is the cost of travelling the transportation distance.
The remaining three terms are the costs associated with operating and maintaining a larger regional warehouse than would be needed with the MFC drone system.These costs are mathematically similar to those in (1) but with one fixed warehouse location.We add to the size and cost of the regional warehouse equal to one MFC, which we can see in the comparison of storage space in Figure 1.Given this formulation, our model implies that the delivery-by-drone system always incurs more warehousing and inventory costs than the traditional van system does.

Stage one: community optimisation
To minimise the delivery-by-drone system cost, (1), we determine the value for N C that sets the first derivative of the function to zero, as in (4).However, this function does not have a closed-form solution (it is not analytically solvable), and so numerical, graphical, or algorithm analysis must be used.These techniques will also confirm that the solution for N C is a minimiser.Later, in the case study, we choose to discretise (4) by evaluating the function over many values of N C (implemented in Excel, version 2110, Build 16.0.14527.20234, 64-bit) to determine an approximate solution to a resolution of ±0.005 MFC units per community.
With the number of MFCs determined after solving (4), other outputs of the model can be determined.No optimisation is conducted on the van system as all the considered parameters monotonically affect the determinate total system cost.

Stage two: the allocation and location problem
Although stage one has determined estimates for the number of MFCs required and other parameters, we have not yet used any spatial knowledge of the region.In stage two, we use the optimal number of MFC per community as an input weight parameter for a spatiallyadapted multiple knapsack optimisation problem (Church and Murray 2008) and then solve a series of single facility location problems.This stage is in place of the clustering step in Tsao et al. (2012) and Chowdhury et al. (2017) as our method preserves more spatially-varying demand data.We highlight that this stage could be solved in several different ways, most interestingly is that it could be solved near manually in a stakeholder-led discussion.Although we present an algorithm process and mathematically exact objective, by separating the problem into two stages then this second stage can be easily presented as a puzzle for discussion; the objectives of compact catchments and one MFC per catchment (as suggested by the sum of N C values) are simple enough to observe on a map and calculate respectively that many stakeholders could understand and attempt their own solutions.We also highlight that we do not use an in-built GIS solver for the determination of the allocation solution, but we do use the 'Mean Center' tool in ArcMAP for the series of single facility location problems.We suggest that this stakeholder method would complement the following algorithmic method in a real-world decision situation.

The multiple knapsack allocation problem
Although the multiple knapsack problem is a conventionally NP-hard problem, modern algorithms can solve instances of thousands of items (communities) and over a hundred knapsacks (MFCs) in a reasonable amount of time (Dell'Amico et al. 2019).Joint locationallocation models, however, typically solve smaller instances of less than one thousand demand points and less than ten facilities because they also must compute transportation routing costs for each considered solution combination (Daskin and Tucker 2018).In other words, there exist more efficient optimisation algorithms for multiple knapsack problems than there exist optimisation algorithms for multi-facility location-allocation problems.This is our motivation to investigate the two-stage method; as a better solution than the classical CA approach can be achieved in less time than the alternative commercial location-allocation algorithms, especially when the problem is large.The typically a-spatial knapsack problem requires a spatial constraint or objective to model a location-allocation problem.This spatial constraint can be formulated in several different ways to enforce varying degrees of strictness to meet the desired catchment characteristics (Church and Murray 2008).In this paper, we use two objectives, catchment perimeter and closeness to MFC unity, in a bi-objective method (objectives expressed in ( 5) and ( 6) respectively).Furthermore, we enforce that every community must be allocated to a catchment area.
Table 2 shows the notation that we use in the multiple knapsack allocation problem.We determine the perimeter of a catchment area (P U,j ) by calculating the sum of the perimeters of the communities allocated to the catchment (P C,i ) and subtracting the internal neighbouring edges (P C,i,h ).An internal neighbouring edge is the amount of a community's perimeter that neighbours another community's when both communities are allocated to the same catchment (g i,h = 1).

Minimise
(5) Equation ( 5) allows us to evaluate the total perimeter of the catchments for a given MFC allocation.We then aim to minimise this total perimeter as this objective makes the catchment areas as compact as possible and, thus, transportation efficient as possible.
The second objective of the bi-objective model ( 6) is to partition the region into MFC catchment zones (subscript label j) such that the sum of the optimal MFCs in all allocated communities (6c) is as close to an even split as possible of the total non-integer number of MFCs required (6b).The number of MFCs for the region is the sum of the MFCs required for each community, first as a real number as in (6a), and second, rounded to the nearest integer as in (6b).The optimal number of MFCs in a catchment is the sum of the optimal number of MFCs in the communities that are allocated to the catchment.We do this by using a binary categorical variable (U i,j ) as in (6c).

Minimise
We also add constraints so that all variables are restricted to non-negative real values.We implemented and solved this stage in a desktop GIS using a greedy algorithm method supplemented by local interchanges, like Hong, Kuby, and Murray ( 2017), so that we preserve the spatial contiguity of the catchment areas.See the supplementary information for further information on the algorithm followed.The next step in the allocation-location process is to determine the optimal MFC locations within the catchment areas.

The series of MFC location problems
With the MFC catchment areas allocated, the location problem becomes a series of independent location problems.The objective is to minimise the transportation distance to and from all demand points in the region.In addition, the assumption of uniform demand within each community may even be relaxed to find a locally optimal solution for MFC location within a catchment.The location may be optimised using any number of solution methods to find the optimal location, but ultimately the decision may be quite restricted by the availability of appropriate space within the catchment area.With the aid of heuristics, the optimal discrete choice may be manually identifiable (especially when given only one MFC catchment area and a finite choice of locations) or the choice may require subjective input of intangible factors of the area.Several discrete optimisation techniques are also available, well known, and available in commercial GIS packages to solve single facility location problems (Daskin 2013;Church and Murray 2008).We use the 'Mean Center' tool in ArcMAP in the case study.

Case study
Now that the method has been described, a future drone delivery system, as potentially deployable by an SSP, is investigated in the city of Calgary.We assume that an Unmanned Aerial System Traffic Management system has been developed and implemented so that the intended operations are legal and safe.The resulting delivery-by-drone system is then compared to a traditional delivery-by-van system in terms of both cost and expected operational emissions.We also benchmark the solution against a commercially available solver.Finally, a sensitivity analysis is conducted on some estimated parameters, and their significance is discussed.The following results are displayed using Esri ArcMAP version 10.7.1 using 2016 geographical and census data made available by Statistics Canada.

Summary of parameters
Table 3 shows a summary of the input variables for the numerical example.See the supplementary information for further supporting information about how we determined the parameter values for the case study.Population density, community area, and logistical sprawl distance vary per community.Additionally, demand uncertainty at the community level is considered unknown.This level of detail in the data may be available to the retailer, as is assumed in Tsao et al. (2012), but as this data is not available in our study area, we use an approximation, (7).This relationship between regional and community level demands, demand variations, and areas is required to maintain the assumption of arborescent network design and maintenance of the stockout factor (Schwarz 1981).This formulation also preserves the expected linear scaling of the optimal MFC number in a community area.

Initial classical analysis
For a baseline understanding of the study area, a classical analysis is first conducted as if the whole region is one community that has slowly varying parameters.This preliminary analysis is also used to set the default resupply frequency.As shown in (1), this rate will have a minimiser for the objective cost.Table 4 shows a cost comparison of the van and drone systems; as shown, the traditional system is almost three times cheaper.On a per item basis, costs are $2.30 and $0.79 for the drone and van systems, respectively.
Consider a few different scenarios from the baseline.First, logistical sprawl in Calgary is relatively short at an average of only 23 km.For this study, we assume that the regional warehouse is in Balzac, about 5 km north of the Calgary city limits, because it is the location of similar existing e-retailer facilities.If we consider a regional warehouse in the same Balzac location was to support a delivery-by-drone system in the city of Red Deer (a smaller city about 130 km north of Calgary) then the sprawl distance would increase, and the cost estimate for the van system overtakes the drone equivalent: $2.36 for the van versus $2.30 for the drone.The expected cost of the drone delivery does not change in the Red Deer scenario as we assume a conservatively high value for the cost to resupply the MFCs which is independent of the logistical sprawl.
Next, we consider changes to the time windows parameter for the city of Calgary.Imposing time-windows to the deliveries only affects the van system (see section 3.4) as we assume the drone deliveries are fast enough, and the standard and expedited deliveries are planned well enough, to always meet the given time-windows.We experimented with the model using different time-windows between one and ten and find that if the van system offers five two-hour time windows in a day, then the cost per item estimate for the van increases to $2.42, which is above the estimate for the drone system ($2.30).This preliminary analysis agrees with intuition and previous research, that these two situations, satellite cities that have large logistical sprawl from the regional hub and in the case of rapid delivery demands, are where the drone MFC system will be cost competitive with traditional van delivery.Furthermore, in section 4.4, we show how the two-stage method can assist in further understanding the time window case.

Two-stage method location-allocation
There are 202 residential communities in the city of Calgary as of the 2016 census.We excluded one community from the case study due to it being non-contiguous with the other 201.The final row, Region, is the sum of the community values and represents the regional value result of the two-stage method.For instance, the optimal number of MFCs for the region is the sum of the optimal number in the communities (17.47) as per (6a).Rounded to an integer this is 17 as per (6b), one lower than recommended by the classical method, which is 18 (rounded down from 18.3).Furthermore, total system cost for the region is calculated similarly, by summation of the community cost estimates.The cost of the baseline drone system is predicted to be $11,013,855 CAD and the van system to be $3,689,137 CAD, which is within 3% and 5% of the classical method, respectively.The total cost is later investigated in a sensitivity analysis by adjusting the input parameters.
Table 5 shows a sample of the results from applying the community optimisation stage of the two-stage method to the relevant community level data to determine the optimal number of MFCs; the table also includes a comparison of the optimal costs for the drone and van systems.
The first column lists the community number, from 1 to the total number, 201.The second column shows the community area, which we determined as the geodesic area of the community using ArcMap Desktop 10.7.1.We also used ArcMap Desktop 10.7.1 to determine the perimeters of the communities and their neighbours using the 'Polygon Neighbours' tool.We obtained population data for each community from the 2016 census (Statistics Canada 2016).Demand density is the population density, calculated from population number and area, multiplied by the mean demand parameter, which is packages per year per person.We calculated the optimal number of MFCs by solving (4) for each community by evaluating the function over many values of N C (implemented in Excel, version 2110, Build 16.0.14527.20234, 64-bit) to determine an approximate solution to a resolution of ±0.005 MFC units per community.This numerical method confirmed that the optimising value of number of MFCs was a cost minimiser.Note that the second stage of the method, where communities are grouped into catchment areas, is most easily conducted when the optimal number of MFCs is lower than one, indeed when it is lower than one-half on average, for each community.The total cost for the drone and van systems is then calculated using (1) and (3), respectively.
The final row, Region, is the sum of the community values and represents the regional value result of the two-stage method.For instance, the optimal number of MFCs for the region is the sum of the optimal number in the communities (17.47) as per (6a).Rounded to an integer this is 17 as per (6b), one lower than recommended by the classical method, which is 18 (rounded down from 18.3).Furthermore, total system cost for the region is calculated similarly, by summation of the community cost estimates.The cost of the baseline drone system is predicted to be $11,013,855 CAD and the van system to be $3,689,137 CAD, which is within 3% and 5% of the classical method, respectively.The total cost is later investigated in a sensitivity analysis by adjusting the input parameters.
We then solve the second stage of the two-stage method, the multiple knapsack problem by implementing the allocation algorithm from section 3.6.1 in ArcMap Desktop 10.7.1.using an Intel i7-9750H CPU @ 2.60 GHz, 64-bit, 16GB RAM computer.We also use ArcMap to evaluate the objective functions during the solution process.Figure 2 shows the resulting location-allocation solution, with MFC catchment areas outlined in thick lines and neighbouring catchments textured, in four different patterns, to visually distinguish them.
We determined the optimal locations of each MFC within the catchments (circles with a centre dot) using the 'Mean Centre' tool in ArcMap, with the demand of each community as a weight.We then obtained the drone transportation distances in ArcMap by first estimating an average travel distance from each community to its allocated MFC using the 'Point Distance' tool from the centre of each community to the MFC; multiplying by the annual package demand of the community; multiplying by two for a round trip.We also calculated the maximum distance any drone in the region would be required to fly by using the 'Construct Points' and then the 'Point Distance' tools in ArcMap to evaluate the maximum distance from the MFCs to their respective boundaries.We found that the maximum distance from an MFC to the perimeter of the allocated catchment found was 11.15, 22.30 km for a return journey, which is less than the maximum flight range (24 km) of the drones and so the location-allocation solution is feasible.
The allocation also allowed us to determine the cycle and safety stock required at each MFC, and the related MFC size and cost.For van travel distance, we selected the location of the regional warehouse (a hollow cross) as the current real-world location of an existing Amazon fulfilment centre.We then determined the logistical sprawl distance from this regional warehouse to each community along road network by using the ArcMap 'Make OD Cost Matrix Layer' tool.The location-allocation solution shown in Figure 2 has objective function values of 1459 km and 0.015 from the functions ( 5) and ( 6), respectively.The transportation distance estimate of 20,470,629 km compared to 20,880,871 km estimated by ArcMap is approximately 2% lower, which is within a typically accepted range when using a CA method.Similarly, the estimated safety stock of 1,969 packages is about 6% higher than the GIS implemented result of 1,856 packages.

Benchmarking
It is important to consider the performance of our two-stage method when compared with commercially available methods.Hence, we performed a location-allocation optimisation in ArcMAP using the Spatial Analyst extension.
For the 510 potential MFC locations, we used the geometric centres of every industrial and commercial zone in the city, sourced from publicly available data online provided by the City of Calgary.The 201 demand points were represented by the geometric centres of every community, weighted by their respective demand determined as in the two-stage method.We set up the drone travel network by creating a layer using the 'XY to Line' tool which connected every demand point to every facility point.The ArcMAP commercial solver cannot account for financial cost, so we set the objective to minimise the sum of the weighted distances from each selected MFC to the allocated communities.This gives a location-allocation solution optimised for minimal last-mile transportation distance.We used this solution as a benchmark for model comparison.
Figure 3 shows the resulting location-allocation solution of the commercial solver for the case of seventeen MFCs.
Table 6 shows the last mile distance and the estimated annual cost of operating the delivery-by-drone system for the respective location-allocation solutions.As commercial solvers typically require the number of facilities to be an input parameter for the analysis, we also performed the commercial solver analysis for eighteen MFCs, as might have been done if the classical method were used to determine the input number of MFCs.
From these results, we can see that the two-stage method better estimates the benchmarks from the commercial solver, the last-mile transportation distance and the total annual cost of the system, than the classical method does.Taking the commercial solver solutions as benchmarks, we can see that the classical analysis overestimates the last mile distance by 9% in the initial estimate, and overestimates by 11% when the solution is implemented into a location-allocation.Whereas the two-stage method does still overestimate  the last-mile distance of the benchmark, it does so by only 3% and 5% for the a-spatial estimate and the implemented estimate, respectively.The result is similar for the total system annual costs.
Although the two-stage method does not produce a lower total cost solution than the commercial solver for our dataset, this is due to the specific cost parameters of this scenario as the last-mile transportation distance represents over half of the total system cost in every scenario.We expect that the two-stage method would produce good solutions across a wider variety of cost variables than the commercial solver, given that the solver minimises only transportation cost whereas the two-stage method accounts for inventory costs.Indeed, we expect there may be scenarios with high holding costs where the twostage method produces a lower-cost solution than a commercial solver.We expect future work to examine the effectiveness of the two-stage method in a wider variety of scenarios.

Operating emissions
Some companies are concerned about the environmental impact of delivery-by-drone systems, and therefore, it is important to evaluate the systems for carbon dioxide emissions.Table 7 shows the carbon dioxide emissions of the last-mile transportation distances for the delivery-by-drone two-stage method solution and the van delivery solution.The vehicle energy per kilometre parameter in Table 7 is from Figliozzi (2017Figliozzi ( , 2020)), and this parameter is also within the range of parameters as studied by Goodchild and Toy (2018) for vans such as a Dodge RAM van (diesel) and a Renault Kangoo EV (electric van).The carbon dioxide emissions per energy parameter uses emissions per fuel estimates published by the Government of Canada in 2017 adjusted for the portion of fuel currently used in the Alberta energy grid, as published by the Alberta Electric Systems Operator in their 2019 report.We assume that diesel and electric vans cost the same and conduct routes in the same manner (same routes, capacity, transportation distance) and we assume that electric vans will only be adopted when they perform similar in cost and routing to the incumbent diesel vehicles; this assumption could be altered and investigated by further studies.We are dominantly interested in comparisons between the drone system and the two van types.
The results support the common assumption that drones are less carbon intensive compared to diesel vans, even in an urban setting.However, the significantly higher total fleet travel distance of the drones caused by their one-package payload capacity ultimately makes the drone system more emissions intensive than a potential future of electrically powered vans using the same energy grid makeup.Use of multi-capacity drones may be able to reduce the total distance travelled (Chen, Hu, and Solak 2021).Our findings are consistent with Figliozzi (2020) and Goodchild and Toy (2018).
Absent from Table 7, however, are the energy requirements of the MFCs in the drone system, as discussed by Stolaroff et al. (2018).We estimate MFC energy requirements for  Alberta using information from a commercial building survey from 2013 published by Natural Resources Canada.The square meterage required for the storage of goods is determined by the two-stage model.The total drone, diesel van, and electric van system emissions in tons of carbon dioxide, including warehouse emissions, is 392.1, 702.7, and 275.7, respectively.These results show that warehousing emissions are a significant factor to consider as they erode the relative benefits of the drone MFC system.

Mixed drone and van system by time windows
Another advantage of the developed two-stage model is that a dual system (vans and MFCs with drones working simultaneously but serving different parts of the same city) can be easily visualised and estimates of optimal joint systems can be created quickly.The classical analysis suggests that with five two-hour time windows, the drone-only system is more cost-effective than the van equivalent.However, this classical analysis is restricted to the scale of the entire city whereas the two-stage method is a partitioning process and can help answer the following questions: can part of the city be served by drone and others by van; and how may the build order of this eventual system be rolled out considering gradually improving time-window offerings to customers?Table 8, Figures 4 and 5 show the two-stage method's answers to these questions.For different delivery window lengths (fractions of a ten-hour delivery day), the optimal cost of drone and van delivery systems can be estimated and compared.W recall that this analysis assumes that the transportation effectiveness of the drone system is unaffected by the time windows due to the effective mixing of standard and expedited packages, whereas the van delivery system is affected as the van delivers fewer parcels per vehicle per tour in shorter time-windows.We add an MFC from one time-window scenario to the next when the cost estimate of the catchment is lower for the drone system than the van system.The decisions to change a catchment area from van service to MFC drone service are independent and so that, under a given set of circumstances, some areas of the city are served by drone and others by van, but no single catchment is served by both van and drone.
Figure 3 shows a mixed mode solution for the 100-minute time window problem.The same catchment areas that were regionally optimal for a full delivery-by-drone system are shown in bold black lines again.These catchment areas are served by one MFC each, which are numbered and represented by white circles with centre dots.Cross-hatched catchments are served by drones, and the non-hatched catchments are served by vans.Moreover, the colour-coding of the communities (outlined with thin grey lines) reflects the difference in estimated cost between the two modes per expected package demanded by that community, according to the legend in the figure.Figure 3 also shows that the e-retailer regional hub is located to the north of the city, above and between catchments 1 and 14.The logistical sprawl distance from the hub to the south of the city is the most significant factor in raising the cost of expedited deliveries as the vans must make this haul multiple times per day, increasing with increasing numbers of time windows.

Sensitivity analysis
We considered three logistical parameters of interest to examine as part of a sensitivity analysis: effective storage density (m), uncertainty in demand (σ ), and resupply frequency (E).These parameters show the three main relationships that parameters have with the decision variable in (4) and the system cost in (1). Figure 6 shows a linear relationship, Figure 8 is a reciprocal function, and Figure 7 is a convex function that is the union of the linear and reciprocal functions.The figures show the last mile cost (last mile drone distance cost), the MFC costs (sum of fixed, resupply, and shelf costs), the holding cost (sum of cycle and safety holding costs), the sum of these costs for total system cost and the related optimal number of MFCs, both rounded ( N * R ) and as the determined CA number (N * R ).We separate the shelf cost from MFC costs in Figure 6 for added clarity to the uncertainty in demand analysis.These three figures are representative of how any single variable may affect the number of MFCs in the two-stage method and are expected when examining the system cost and derivative equations as in (1) and (4).
Of these three parameters, demand uncertainty is the parameter that the seller has the least influence over as it is a function of customer behaviour.This factor is shown to near linearly relate to higher optimal drone MFC system cost as less knowledge of customer demand leads to larger safety stock inventory, which leads to increases in both storage space and inventory holding costs.Alternatively, the seller may accept more frequent stock outs and accept the associated indirect costs to branding and customer loyalty.To avoid either increased stock costs or indirect customer costs, our analysis and Figure 6 suggests that decreasing the number of MFCs to utilise the consolidation effect is appropriate (Geoffrion 1979).Although decreasing the number of MFCs will increase last-mile transportation costs, this action will also create larger catchment areas, which will aggregate more demand.A greater aggregation of demand reduces the increase in the uncertainty in demand for the expediated products relative to not changing the number of MFCs.This can be seen analytically in (7), in section 4.1.This lower increase in uncertainty in demand per MFC catchment, created by expanding the catchment areas, helps contain the increase  in total safety stock required when summated across the region, thereby preventing superlinear increases in inventory holding costs.Our results show that, in this case study, the increased holding costs avoided by decreasing the number of MFCs is greater than the trade-off increase in last-mile transportation costs for decreasing the number of MFCs, therefore making decreasing MFCs the correct choice in response to higher demand uncertainty.Average package value (u) also affects parts of the objective function linearly, as understood by examining (4), and is a similarly significant factor to consider which affects optimal number of MFCs in the same manner.
Figure 6 can also inform how much effort and resources should be put towards reducing demand uncertainty.For example, the seller could issue a customer survey with the expectation that this survey will reduce demand uncertainty (the standard deviation in predicted user demand) by 0.2 packages per person per year.Figure 6 suggests that this reduced uncertainty could save approximately $500,000 CAD per year in system costs.This cost is saved because the reduced uncertainty, as provided by the survey, allows for the opening of two more MFCs and a more than commensurate reduction in transportation costs.Thus, such a survey may be a positive investment if the survey and the resulting non-recurring implementation work associated with it could be conducted for less than the expected amount saved.
Figure 7 shows the resupply frequency, which also affects order quantity.Within realistic operational ranges, the system cost is relatively flat, which suggests that other tactical or operational constraints should be considered to determine the resupply frequency between deliveries once weekly and once every three days, as within these ranges the system cost does not increase or decrease by more than 3%.However, over this same range, the number of MFCs to facilitate this system cost changes significantly, from 19 to 15. Consequently.this analysis suggests that an effective strategy could be to establish 19 MFCs which have fewer deliveries (once weekly) at the start of a multi-year plan, and then the seller should expect to increase the resupply frequency over time to accommodate demand growth.This strategy allows for longer-term growth while still achieving a near optimal estimated cost in the short-term.The resupply frequency may also be constrained by the maximum size and thus minimum frequency of resupply trucks needed to service the product of the size and demand for the products sold.The effective storage density in the MFC is likely to be lower than the storage density in a truck or van during transit.
Figure 8 shows that changes in effective storage density above 15 packages per square meter do not result in any changes in the number of MFCs.The system cost savings approach a limit of about 8% beyond this effective storage density.Moreover, effective storage densities above some value will become unreasonable given minimum commercial building lot size, demand, and resupply frequency.This 'cut-off' value, 15 packages per square meter in this case study, however, is affected by the per area shelf space cost, and it is clear by examining the cost function in (1) that a higher shelf space cost leads to a higher storage density cut-off value.Consequently, we define the cut-off value as when decreases in shelving cost become insignificant.Since the relationship is linear, there is also an effective shelf cost per item (C s /m) cut-off value.Our case study and analysis suggest that delivery-by-drone systems with an inventory management strategy and rental market combination that can support an effective shelf cost above 0.09 CAD per item are insensitive to further decreases in rental cost or increases in effective storage density.
Effective storage density is achieved by: accurate prediction of demand at a time scale relative to resupply frequency, measured in days; small item sizes; efficient rack placement; improved warehousing technology; and good inventory management.If these densities cannot be achieved, then the effective storage density factor is significant because it is the only parameter that results in MFC numbers below ten and even as low as one or two MFCs, an unfeasible MFC arrangement for current drone technology maximum flight range.Thus, confidence in selling goods stocked in MFCs and achieving a high effective storage density are initial, but not continuing, barriers for companies using a delivery-by-drone system.Furthermore, an effective storage density above a certain range (above 15 in this case) is effectively wasted as the optimal system, as estimated by the first stage in the CA, is below the minimal lot size of commercial space to rent or buy observed in our study, and thus cannot be realised.The SSP may instead increase the variety of items offered for expedited delivery until this threshold of effective storage density is reached as this will garner more sales, and more profit, from their customers.Understanding the relationship shown by Figure 8 is important to not over or under offer the range of expedited delivery packages.

Conclusion
Urban goods delivery has been dominated by fossil fuel powered, human driven vans and customer-attracting physical retail stores of many sizes for the history of traditional retail.With the advent of disruptive technologies, such as drones and e-commerce, and the adoption of disruptive business practices, such as seller and service providers and micro-fulfilment centres, the traditional structure of urban retail logistics will change.Multiple recent literature reviews on urban logistics facilities and urban delivery-by-drone have shown the research gaps that this paper addresses.
This paper has shown that, for the example seller and service provider studied, a microfulfilment centre delivery-by-drone system for typical packages is not yet justifiable on a pure cost effectiveness basis but may be in the future given advances in technology, regulations, and/or implementation of short time-window delivery.We have also shown that the mix of expedited and standard packages the seller will offer is an essential factor for success of a delivery-by-drone and MFC system, and conversely this mix of goods can be a critical barrier if managed incorrectly.This mix of goods affects both drone utilisation efficiency outside of MFCs and effective storage density within each MFC, a factor that determines the size of the MFCs needed.If these challenges are overcome, a cost-optimal drone system can be a lower emissions alternative to a traditional diesel van system but not lower than a future electric van system.
Furthermore, decision makers in this field need adaptable tools, such as CA models, which will aid in the understanding of these new logistical structures.Within a city, however, the common first assumption of uniform demand density that many CA models make is often invalid due to spatially varying socio-economic factors.We have shown that CA methods can be used to understand this coming urban disruption while acknowledging non-uniform demand density across space; in addition, CA methods can be used to build allocation maps in a methodological way with commercially available software.Therefore, we propose that this two-stage method can be a less 'black-box' methodology for stakeholders to understand than alternative methods whilst still being an effective planning estimate tool.This paper has shown that this two-stage method estimates transportation distance and inventory to a better accuracy than classical single uniform demand approximation methods, and thus the insights of CA methods at a local level will approximately hold for cost-optimal regional level solutions.

Future work
We assumed drones work constantly over their daily routine, meaning that they were delivering standard parcels in-between time-sensitive deliveries.This intermixing of expedited and standard deliveries requires effective operational algorithms which have yet to be developed and are of interest for future study.
Another algorithm development could be improved and automated allocation procedures for the multi-knapsack allocation problem, leveraging the Python programming capabilities of ArcMap or comparable languages for other GIS software.
A strength of the two-stage method is that the objective cost function is independent from the allocation stage and future work that changes the objective function, to either account for different delivery modes or include emissions directly in the objective, are next of interest.

Figure 1 .
Figure 1.Conceptual diagrams of Van System (Top) versus MFC and Drone System (Bottom), Highlighting storage space needs.

Figure 4 .
Figure 4. Drone, Van, and mixed system costs by number of time windows.

Figure 6 .
Figure 6.Regional uncertainty in user demand.

Table 1 .
Model features and paper contributions.By incorporating CA techniques into an integrated strategic model, we allow for future adaptations and extensions of the work, such as in truck-drone hybrids and sidewalk autonomous delivery robots.2a.1.It can be extended to account for parallel and hybrid delivery systems where vans, drones, sidewalk robots, cargo-bikes, and other modes are deployed from one location.

Table 3 .
Case study default parameters.

Table 4 .
Regional uniform demand results.

Table 5 .
Example community optimal MFC number results.

Table 6 .
Benchmarking of methods.

Table 7 .
Energy usage by vehicle type considering two-stage model estimate solution.

Table 8 .
MFCs by time windows.