WrEXfuZv

1Contents lists available at ScienceDirect
2Int. J. Production Economics
3journal homepage: www.elsevier.com/locate/ijpe
4Drones for disaster response and relief operations: A continuous
5approximation model
6Sudipta Chowdhurya, Adindu Emelogua, Mohammad Marufuzzamana,⁎, Sarah G. Nurreb,
7Linkan Biana
8a Department of Industrial and Systems Engineering, Mississippi State University, Starkville, MS 39759-9542, USA
9b Department of Industrial Engineering, University of Arkansas, AR 72701, USA
10A R T I C L E I N F O
11Keywords:
12Drone routing
13Humanitarian logistics
14Continuous approximation
15Facility location
16Inventory management
17A B S T R A C T
18This paper proposes a Continuous Approximation (CA) model that designs the potentiality of drones as a mode
19of transportation to supply emergency commodities in a disaster-affected region. The model determines the
20optimal distribution center locations and their corresponding service regions and ordering quantities to
21minimize the overall distribution cost for the disaster-relief operation. We propose a two-phase CA approach
22that solves the model efficiently. We conduct extensive sensitivity analysis to reveal insights into how system
23design varies with key drone design parameters. We use three disaster-prone coastal counties of Mississippi as a
24testbed to visualize and validate the modeling results.
251. Introduction
26Between 1995 and 2015, Emergency Events Database (EM-DAT)
27of The Center for Research on the Epidemiology of Disasters (CRED)
28recorded a total of 6457 weather-related disasters, which claimed over
29606,000 lives and left 4.1 billion people injured, homeless, and/or in
30need of emergency assistance. More specifically, between 2005 and
312014, EM-DAT recorded an average of 335 weather-related extreme
32events which is approximately 14% higher than the number recorded
33between 1995 and 2004 and twice the number between 1985 and 1995.
34Scientists predict that we will witness a continued upward trend in
35weather-related extreme events in the decades ahead (UNISDR, 2015).
36Among all the countries in the world hit by the highest number of
37weather-related disasters, United States (US) comes first, with 472
38disasters occurring over the past 20 years (Shaw, 2015). In recent
39years, natural disasters such as flooding, wildfire, hurricane, tornadoes,
40and many others hit frequently in different parts of the US and cause
41significant damages to both economy and human life. The most
42devastating natural disaster that shocked the entire country in recent
43years was Hurricane Katrina which rampaged the Gulf Coast of the US
44in the morning of August 29, 2005. It had a Category 3 rating on the
45Saffir-Simpson Hurricane Scale and it brought sustained winds of 100–
46140 miles per hour.1 The storm itself did a great deal of damage, but its
47aftermath was even more catastrophic. Many people charge that the
48federal government reacted slowly to meet the needs of the disaster
49affected people. The main reason behind this is that the government,
50particularly the federal government, seemed unprepared for the
51disaster. It took the Federal Emergency Management Agency (FEMA)
52days to establish operations in the disaster-affected regions, and even
53then did not seem to have a sound plan of action. People suffered
54significantly because of the scarcity of clean water, food, and medical
55aid. A proper plan regarding facility location and inventory policy of
56emergency products as well as deciding on which mode of transportation
57to use would have been much more beneficial during the recovery
58phase of the disaster and helped people in the affected region
59significantly.
60Logistics activities in response to a disaster are commonly known as
61humanitarian logistics. Humanitarian logistics can be defined as the
62process of planning, implementing, and controlling the efficient, costeffective
63flow and storage of goods and materials, as well as related
64information, from a point of origin to a point of consumption for the
65purpose of meeting the end beneficiaries requirements (Thomas and
66Mizushima, 2005). The aim is to deliver the right commodity at the
67right time, to the right place, and to the right people. Although the basic
68methodology behind humanitarian logistics is the same as classical
69logistics operation, the effect of decision making holds more significance
70for humanitarian logistics. Note that while planning a robust
71logistics network, decisions regarding where to place emergency
72http://dx.doi.org/10.1016/j.ijpe.2017.03.024
73Received 1 October 2016; Received in revised form 20 February 2017; Accepted 30 March 2017
74⁎ Corresponding author.
75E-mail addresses: sc2603@msstate.edu (S. Chowdhury), aae39@msstate.edu (A. Emelogu), maruf@ise.msstate.edu (M. Marufuzzaman), snurre@uark.edu (S.G. Nurre),
76bian@ise.msstate.edu (L. Bian).
771 Available from: http://www.cnn.com/2005/WEATHER/12/21/katrina/.
78International Journal of Production Economics 188 (2017) 167–184
79Available online 12 April 2017
800925-5273/ © 2017 Elsevier B.V. All rights reserved.
81MARK
82facilities from where all the commodities will be transported to the
83demand sites is of significant importance as it has a direct impact on
84how fast disaster-affected points can be served. Simultaneously,
85inventory policy decisions are equally important since they determine
86whether enough commodities are available to serve all the disaster
87affected regions.
88In this study, we propose an integrated facility location-inventory
89allocation model for a disaster affected region where drones can be
90considered as a potential mode of transportation to transport emergency
91supplies to the demand points. Drones can become very useful as
92a mode of transportation in humanitarian logistics since they do not
93need any preexisting path to fly. Thus, if a natural disaster strikes and
94roads are blocked, drones can easily be used to serve a disaster affected
95region. Nowadays drones are used in a diverse range of civilian
96activities, primarily humanitarian aid and rescue actions in various
97natural disasters which include earthquakes, hurricanes, flooding,
98volcanic eruption, tsunamis, and many others. The crisis management
99of these natural disasters poses very critical decisions related to the
100well-being of the affected people. Integrating drones into emergency
101and disaster response protocols can be considered as an effective tool to
102provide immediate relief benefits needed by civilians, communities,
103and first responders. The disaster relief life cycle can be split into four
104major stages: prevention, preparation, response, and recovery (Center
105for Disaster, 2016). Drones have the potential to play a major role in all
106four stages, though currently are used overwhelmingly in the response
107stage. Unfortunately, even though they have the potential to respond
108quickly in a critical situation, drones are surprisingly under-utilized.
109When a disaster strikes, drones are able to provide support with
110risk assessment, mapping, and planning for the affected region
111(Measure-Red cross, 2015). Al-Tahir et al. (2011) present an overview
112and assessment of the technology relevant to low cost cameras and
113platforms to acquire aerial photographs of a disaster-affected
114Caribbean region. The authors suggest that drones can be employed
115further to establish temporary communication structures, creating upto-
116date maps of the affected region, and searching for hot spots where
117the rescue teams may have more chances of finding victims. There are a
118number of cases where drones have already been used in humanitarian
119settings. For instance, in 2013, a US start-up company called Matternet
120announced that it had tested humanitarian aid drone prototypes in
121Haiti and the Dominican Republic (Cohen, 2014). Drones were further
122used to create up-to-date maps of the devastated areas following
123Typhoon Haiyan in 2013 where the mapping efforts were coordinated
124between multiple aid and drone organizations (Greenwood, 2015).
125Similarly, when the earthquake and Tsunami hit the Fukushima area of
126Japan, the nearby nuclear power plant named Fukushima Daiichi
127plant was highly affected. A drone was launched to fly around the
128Fukushima plant to gauge the radiation levels. It took off around six
129kilometers away from the damaged power plant in Namie City. It had a
130flying time of 30 min and it was able to collect radiation levels in realtime
131and went back to the scientists safely (Pamintuan-Lamorena,
1322014). Table 1 summarizes the application of drones in response to
133natural disasters in recent years between 2005 and 2015 (Murphy,
1342014).
135In situations, where transportation network is highly impacted by
136natural disasters, obvious mode of transportation e.g. trucks often can't
137be used in response operations. Despite high expectations of the role of
138drones to serve in humanitarian logistics, there is not yet a comprehensive
139design framework on how to economically deploy a drone that
140is capable of serving a disaster affected region considering the restrictions
141provided by Federal Aviation Administration (FAA)2 (Federal
142Aviation, 2016) and the technological limitations (e.g., limited battery
143and weight carrying capacity) that the existing drones hold. This paper
144proposes a continuous approximation (CA) model for determining the
145optimal configuration of a humanitarian logistics system, including the
146location of distribution centers and the corresponding emergency
147supply inventories at these centers under stochastic demand and drone
148flying range constraints. This optimal configuration shall minimize the
149overall system cost including the costs of opening distribution centers,
150ordering and holding costs of emergency supplies, and the expected
151transportation costs of using trucks and drones to serve a disaster
152affected region. We approximate the nonlinear drone transportation
153cost by considering a number of routing specific factors such as
154climbing, hovering, descending, turning, acceleration and deceleration,
155rotational, and constant speed cost into account. We then propose a
156two-phase continuous approximation approach to solve this nonlinear
157programming problem efficiently. Finally, we apply this model to three
158disaster prone coastal counties of Mississippi i.e., Hancock, Harrison,
159and Jackson counties and draw interesting managerial insights into the
160optimal system design and the total system cost (Table 2).
161Table 1
162Drones used in response to natural disasters between 2005 and 2015 (Murphy, 2014).
163Year Natural disaster Name of drone Application
164A B C D
1652005 Hurricane Katrina Response
166(USA)
167AeroViroment Raven ✓ ✓
168Evolution ✓ ✓
169iSENSYS T-Rex ✓ ✓
170Silver Fox ✓ ✓
1712005 Hurricane Katrina Recovery
172(USA)
173iSENSYS IP3 ✓
1742005 Hurricane Wilma (USA) iSENSYS T-Rex ✓ ✓
1752007 Berkman Plaza II iSENSYS IP3 ✓
1762009 Laquila Earthquake (Italy) Custom ✓ ✓
1772009 Typhoon Morakot (Taiwan) Unknown ✓
1782010 Haiti Earthquake (Haiti) Elbit Skylark ✓
1792011 Christchurch Earthquake
180(NZ)
181Parrot AR. Drone ✓
1822011 Tohoku Earthquake (Japan) Pelican ✓
1832011 Fukushima Nuclear
184Emergency (Japan)
185Custom ✓
186Honeywell T-Hawk ✓ ✓
1872011 Evangelos Florakis
188Explosion (Cyprus)
189AscTec Falcon ✓ ✓
190AscTec Hummingbird ✓ ✓
1912011 Thailand Floods (Thailand) FIBO UAV-1 ✓
192FIBO UAV Glider ✓
193SIAM UAV ✓
1942012 Finale Emilia Earthquake
195(Italy)
196NIFTi ✓
1972013 Typhoon Haiyan
198(Philippines)
199unknown ✓
2002013 Lushan Earthquake (China) HW18 (Ewatt
201HoverWings)
202✓ ✓
2032013 Boulder Colorado floods
204(USA)
205Falcon Fixed ✓
2062014 SR350 Mudslides Response
207(USA)
208DJI Phantom ✓
209AirRobot 100 ✓
210Precision Hawk ✓
2112014 SR350 Mudslides Recovery
212(USA)
213AirRobot 180 ✓
214Precision Hawk ✓
2152014 Balkans flooding (Serbia,
216Bosnia-Herzegovina)
217ICARUS custom ✓ ✓
2182014 Collbran landside (USA) Falcon Fixed ✓ ✓
219Falcon Hover ✓ ✓
2202014 Yunnan China Earthquake
221(China)
222Parrot AR Type 2 ✓
223A=Search, B=Reconnaissance and mapping, C=structural inspection, and D=estimation
224of debris.
2252 In the United States, Federal Aviation Administration (FAA) requires UAVs to be
226operated under a ceiling of 400ft which will severely limit the effective range of utilizing
227drones for disaster relief operation (Federal Aviation, 2016)
228S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
229168
230In summary, the main contributions of this paper are as follows:
231• As far as the author's know, this study is the first work that
232addresses design of a humanitarian logistics system by considering
233drone as a viable mode to transport emergency supplies to the
234disaster affected regions.
235• The nonlinear drone transportation cost along with the inventory
236cost decisions are effectively modeled. We then employ a two-phase
237continuous approximation approach to solve the nonlinear distribution
238problem.
239• Provided a real life application to visualize and validate the modeling
240results and deduced a number of managerial insights about different
241key drone parameters (e.g., drone flight height, speed), road access
242availability following a natural catastrophe, and how inventory
243holding and reordering cost impact the humanitarian logistics
244system performance.
245The exposition of this paper is as follows. In Section 2 we review the
246literatures that are most pertinent to this research. In Section 3, we
247present a continuous approximation model to supply emergency
248commodities in a disaster affected region. In Section 4, we outline
249the developed two-phase continuous approximation approach to solve
250our optimization model. We then present the model input parameters
251in Section 5 and a series of experimental results in Section 6. Lastly, we
252conclude and present avenues for future research in Section 7.
2532. Literature review
254The proposed research is built upon existing studies on facility
255location and inventory allocation problems which researchers have
256attempted to solve over the years. Nozick and Turnquist (1998)
257perform regression analysis to approximate the safety stock cost of a
258distribution center and incorporate the inventory decisions in a fixedcharge
259facility location model. Nozick and Turnquist (2001) extend
260their prior work by developing a model that simultaneously integrates
261facility costs, inventory costs, transportation costs, and service responsiveness
262by borrowing ideas from queuing theory, discrete choice
263location analysis, and multi-objective decision-making. Erlebacher
264and Meller (2000) propose a two-stage heuristic procedure to identify
265the interactions that exist between location and inventory in designing
266distribution system. Shen et al. (2003) propose a formulation to
267integrate facility location and inventory allocation decisions under
268the same decision making framework to minimize the overall system
269cost. Note that the formulation incorporates a nonlinear inventory and
270safety stock cost. Miranda and Garrido (2004) propose a non-linear
271mixed-integer programming model that incorporates inventory control
272decisions, i.e., economic order quantity and safety stock decisions into
273typical facility location models. The authors employ a Lagrangian
274relaxation approach to solve their proposed nonlinear optimization
275problem. Teo and Shu (2004) study a distribution network design
276problem integrating transportation and an infinite horizon multiechelon
277inventory cost function. The authors consider the trade-off
278that exists between inventory, direct shipment, and facility location
279cost for such a system. Romeijn et al. (2007) propose a framework for
280the two-echelon supply chain design problem that incorporates location
281decisions as well as location-specific transportation and inventory
282decisions. Fairly recently, Tsao et al. (2012) have developed a detailed
283non-linear cost model where the authors propose a solution technique
284that preserves the interrelations between facility location and inventory
285management decisions.
286Disaster preparedness and response has been studied extensively by
287many researchers in recent years. Ozdamar and Ertem (2015) have
288reviewed different logistics models developed for the response and
289recovery planning phases considering their modeling and formulation
290characteristics. They also investigate the technological advances that
291facilitate the execution of proposed models and solutions on various
292types of information systems. This in turn leads to the conclusion in
293order to achieve unobstructed international humanitarian cooperation,
294standard software that integrates different planning phases is highly
295needed. Luis et al. (2012) propose an analysis of the use of operations
296research models in transportation of relief goods from the perspective
297of both practitioners and academics. Authors conclude in this study
298that, if models are flexible enough to address the high level of
299uncertainty in modeling disasters, the same framework can also be
300carried over into other areas with similar challenges. A novel framework
301is proposed by Sahebjamnia et al. (2015) that is developed for
302integrated business continuity and disaster recovery planning for
303efficient and effective resuming and recovering of critical operations
304after being disrupted. In this study, at the strategic level, the context of
305the organization is first explored and primary features of the organiza-
306Table 2
307Summary of key parameters and variables used in this study.
308Notes unit
309Set
310I Set of clusters (i.e. i = 1, 2, …, 12)
311Symbol
312Nri
313Total number of distribution centers –
314Decision variables
315Ari
316Influence area mile2
317Qri
318Ordering quantity –
319Parameters
320fr Cost of opening a refueling station for trucks and battery
321charging station for drones
322$
323gr Constant scalar –
324γi Discrete demand points in each cluster Ci expressed as a
325spatial density slow varying function
326–
327θi Customer demand for each cluster Ci expressed as a spatial
328density slow varying function
329–
330δ Planning horizon days
331p Probability of using truck to transport emergency supplies –
332q Probability of using drone to transport emergency supplies –
333αp Transportation cost per mile per item using truck $
334αq Transportation cost per mile per item using drone $
335Pcl Required power of the motor while climbing watt
336v Speed of drone while climbing miles/h
337k Altitude of drone mile
338Cc Cost per unit of energy consumed $
339β Minimum power needed to hover over the ground watt
340λ Motor speed multiplier –
341th Hovering time s
342Pdsv′ Required power of the motor while descending watt
343Speed of drone while descending miles/h
344Pacc Required power of drone during acceleration watt
345Pvc Required power of drone while maintaining constant speed watt
346Pdec Required power of drone during deceleration watt
347T1 Time spent during acceleration phase s
348T2 Time spent during constant speed phase s
349T3 Time spent during deceleration phase s
350m Total payload of a drone including its own weight kg
351d Distance covered with acceleration acc mile
352T Time taken to cover distance d covered with acceleration acc s
353w Percentage of total distance covered during acceleration phase –
354b Percentage of total distance covered during constant speed
355phase
356–
357y Percentage of total distance covered during deceleration phase –
358Pturn Required power of drone during turning watt
359ηtΔθ number of turns a drone takes in a given time period –
360Angle covered during turning radian
361wturn Angular rotation speed rad/s
362Rr Reorder cost $
363μr Mean lead time days
364σr
3652 Variance of lead time days
366αr Service level at distribution center –
367hr Holding cost $
368Rd Range of drone miles
369S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
370169
371tional resilience are recognized. Afterwards, a new multi-objective
372mixed integer linear programming model is formulated to allocate
373internal and external resources to both resuming and recovery plans
374simultaneously. Finally, hypothetical disruptive events are investigated
375to evaluate the applicability of the plans at the operational level. A
376stochastic linear mixed-integer programming model with integrated
377decision considering facility and stock pre-positioning, evacuation
378planning, and relief vehicle planning is proposed by Manopiniwes
379and Irohara (2016). In this work, authors suggest a compromise
380between the cost and equity of relief victims. Ransikarbum and
381Mason (2016) propose a goal programming-based multiple-objective
382integrated response and recovery model to investigate strategic supply
383distribution and early-stage network restoration decisions. Given
384limited capacity, budget, and available resources, this model prescribes
385equity- or fairness-based compromise solutions for user-desired goals.
386Numerical experiments in this study shows that goal programming
387based model provides a compromise solution when no solution exists
388that satisfies hard or real constraints. Tofighi et al. (2016) propose a
389two-stage scenario-based possibilistic-stochastic programming approach
390for designing a relief network in Tehran in preparation for
391potential earthquakes to cope with the main logistical problems before
392and after disaster strikes. During the first stage, the locations for
393central warehouses and local distribution centers are determined along
394with the pre-positioned inventory levels for the emergency supplies. In
395the second stage, a relief distribution plan is developed based on
396various disaster scenarios with the objective of minimizing total
397distribution time, the maximum weighted distribution time for the
398critical items, total cost of unused inventories, and weighted shortage
399cost of unmet demands. A tailored differential evolution algorithm is
400also developed to find feasible solutions within a reasonable CPU time.
401Jia et al. (2007) create models and solution approaches for determining
402the facility locations of medical supplies in response to large-scale
403emergencies. The authors include the uncertainty in demand and
404insufficient medical supply by ensuring each demand point receives
405services from multiple facilities located at different distances from the
406demand points. Using a maximal covering problem, they incorporate
407both the quantity-of-coverage and quality-of-coverage. They develop
408three different heuristic approaches to determine the best facility
409locations. Considering various factors such as demand locations, the
410transportation network condition after a disaster, the locations of
411existing hospitals, and the number of available temporary facility
412locations, Chen et al. (2016) propose an integer programming model
413to improve the effectiveness of emergency medical service after a
414disaster by determine the locations for EMS facilities. They use
415Lagrangian relaxation to enable the ability to solve large, realistic sized
416problems. A two-stage stochastic mixed integer program (SMIP) is
417presented by Rawls and Turnquist (2010) for an emergency response
418pre-positioning strategy for different disasters. The authors investigate
419effect of uncertainty in the inventory level and transportation network
420availability after a disaster strikes in a mathematical model. This model
421is then solved using a customized Lagrangian L-shaped method. Noyan
422(2012) has developed a risk-averse two-stage stochastic programming
423model for disaster management. Here, the author has specified the
424conditional-value-at-risk as the risk measure and has developed two
425decomposition algorithms based on the generic Benders-decomposition
426approach to solve such problems. Both the concepts of the value of
427perfect information and the value of the stochastic solution are utilized
428for the proposed framework.
429Transportation of emergency supplies to a disaster affected region
430has been an area of interest among researchers for many years.
431Different researchers have explored this area considering a vast array
432of factors. Ben-Tal et al. (2011) propose a robust optimization model
433for dynamically assigning emergency response and evacuation traffic
434flow problems with time dependent demand uncertainty. Hamedi et al.
435(2012) propose a genetic algorithm-based heuristic to solve a humanitarian
436response planning problem for a fleet of vehicles with
437reliability considerations. Yuan and Wang (2009) present two mathematical
438models where the first model is a single-objective path selection
439model that minimizes the total travel time along a path. The travel
440speed on each arc is modeled as a continuous-time decreasing function.
441The second model extended the first model to consider a case where
442chaos, panic, and congestion in the time of a disaster are taken into
443consideration. The authors propose an ant colony algorithm to solve
444their proposed optimization models. Yan and Shih (2009) developed a
445multi-objective, mixed-integer, multiple-commodity network flow problem
446that minimizes the length of time required for both emergency
447roadway repair and relief distribution. Oran et al. (2012) present novel
448formulations of the facility location problem and vehicle routing
449problem with time windows with priority considerations. This means
450that the proposed models ensure that higher priority locations are
451considered first, and then lower priority ones are considered. This goes
452for both facility and routing decisions. The facility location problem is
453solved using an MIP solver, while a tabu search based metaheuristic is
454developed for the solution of the vehicle routing problem with time
455windows. Yi and Özdamar (2007) propose an integrated locationrouting
456model for coordinating logistics support and evacuation
457operations after emergencies and natural disasters. The objective of
458this study is to maximize response service level by enabling quicker
459transportation to affected areas and locating temporary emergency
460units in potential locations. The resultant model generates queues in a
461minimized way through the full utilization of facility capacities
462achieved by the interaction of the routing problem with service rate
463equilibrium and location. Ji and Zhu (2012) present an optimization
464framework that quantifies the risk of disasters as well as the support
465mechanisms for the disaster relief logistics. Interested readers can
466review the studies of Abounacer et al. (2014), Sheu et al. (2005), Sheu
467(2007), Zhang et al. (2012), Balcik and Beamon (2008), Lin et al.
468(2009), and Rath and Gutjahr (2014) for a detailed discussion about
469the different variants of location-routing problems applied to a disaster
470relief operation.
471Another stream of research focuses on locating facilities in different
472disaster response problems. Ahmadi et al. (2015) propose a multidepot
473location-routing model to solve a last mile delivery problem for
474serving a disaster affected region. A novel Variable Neighborhood
475Search (VNS) algorithm is developed to solve large instances of the
476proposed model. Khayal et al. (2015) propose a location-allocation
477model that transfers excess resources between temporary facilities
478operating in different time periods with the goal of minimizing
479deprived victims affected by disasters. One uniqueness of this model
480is that it allows delayed satisfaction of demand when resources in a
481planning period are not enough to fulfill the customer demand. Salman
482and Yücel (2015) study the problem of locating emergency response
483facilities in such a way that majority of the demand can be satisfied in a
484reasonable amount of time. The authors developed a tabu search
485heuristic to solve the proposed optimization model. Afshar and
486Haghani (2012) propose a mathematical model in response to natural
487disasters that controls the flow of several relief commodities from
488multiple sources to demand points. The model further considers
489various level of details such as vehicle routing schedules for relief
490commodities and the optimal locations of temporary facilities. Kılcı
491et al. (2015) propose a mixed-integer linear programming model that
492simultaneously controls the utilization of the shelters and assigns
493population to operating shelter areas. This model has been primarily
494designed to assist the planning phase for a possible disaster recovery
495process. Rennemo et al. (2014) propose a three-stage stochastic
496programming model that considers the opening of local distribution
497facilities, initial allocation of supplies, and last mile distribution of aid.
498The model further captures multi-commodity flows and transportation
499modes while explicitly considers different vehicle types. The authors
500use utility as a measure to achieve a fair distribution of aid to serve a
501disaster affected region.
502The majority of the studies discussed above use discrete models to
503S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
504170
505formulate the integrated network design problem. These models often
506provide unrealistic solutions when the availability of data is limited.
507Further, the increased volume of data will add computational complexity
508for the discrete models. A Continuous Approximation (CA)
509approach can be employed to overcome this problem since the method
510requires less data to approximate solutions. This approach defines
511decision variables in terms of continuous functions which in turn
512reduces the complexity of the model. Newell (1973) has demonstrated
513the idea of applying continuous techniques to finite-dimensional
514operational research problems. Daganzo (1996) has addressed a
515number of different transportation-related problems using a CA
516approach. Blumenfeld and Beckmann (1985) have developed an
517analytical framework for estimating the cost of distributing freight
518from one origin to many destinations. The analysis uses a continuous
519space modeling approach, which requires only the spatial density of
520destinations and the average and variance of demand. Langevin et al.
521(1996) present an overview of CA models that have been developed for
522freight distribution problems. Geoffrion (1976) has developed a continuous
523model to serve warehouse demands which are distributed
524uniformly over a plane. Erlenkotter (1989) has developed a model to
525determine the optimal area served by a single production unit where
526demand is assumed to be distributed uniformly over a market area.
527Rutten et al. (2001) have extended the model proposed by Erlenkotter
528(1989) to develop a more precise model by approximating accurately
529the fixed location, inventory, and transportation costs. The aim is to
530determine the optimal number of depots serving a set of uniformly
531distributed customers in a given area.
532An analytic method is developed by Burns et al. (1985) that uses the
533spatial density of customers to minimize the freight transportation and
534inventory costs. This study analyzes and compares two distribution
535strategies: direct shipping (i.e., shipping separate loads to each
536customer) and peddling (i.e., dispatching trucks that deliver items to
537more than one customer per load). Dasci and Verter (2001) present a
538CA framework to represent spatial distributions of cost and customer
539demand. The proposed CA model provides a number of insights about
540the impact of problem parameters on facility design decisions.
541Wangand and Lu (2006) have studied spatial modeling and proposed
542smoothing techniques for non-homogenous processes by considering
543details at different levels of the distribution network. Further, this
544study has proposed fine refinements to the approximation model based
545on the level of details captured by the data. A CA framework is
546developed by Murat et al. (2010) where the market demand is modeled
547as a continuous density function. This methodology prioritizes the
548allocation decisions rather than the location decisions and a local
549search heuristic (steepest-descent algorithm) is employed to solve the
550optimization model. Pujari et al. (2008) have implemented a CA
551approach to determine the optimal number of shipments along with
552their sizes while location, production, inventory, and transportation
553decisions are taken into consideration. Different from other works, our
554approach provides a first comprehensive design framework on how to
555economically deploy drones to serve a disaster affected region. The
556model provides numerous managerial insights which allow the decision
557makers to understand how different key drone specific parameters
558(e.g., drone flight height, speed), road access availability following a
559natural catastrophe, and inventory holding and reordering cost impact
560the overall system design performance.
5613. Problem description and model formulation
5623.1. Problem description
563In this section, we develop an integrated mathematical model of a
564logistics network for a disaster affected region that employs drones and
565trucks as the potential modes of transportation to transport emergency
566supplies from multiple source points (e.g., distribution centers) to
567destination points (e.g., customer demand points). The objective is to
568minimize the expected overall cost which includes the cost of locating
569facilities(i.e. distribution centers), allocating inventories, and transporting
570emergency supplies using trucks and drones. The proposed
571model and the solution approach have been designed for the disaster
572preparedness phase. Fig. 1 illustrates the service area of the proposed
573logistics network for a distribution center.
5743.2. Model assumptions
575Before we present our logistics network design problem, let us first
576introduce a number of assumptions that we made related to the
577network structure, demand pattern for emergency commodities, availability
578of modes of transportation, and inventory replenishment
579policies at the distribution centers. Note that developing a humanitarian
580logistics network is a challenging job in itself, and incorporating
581drones as a potential mode of transportation further complicates the
582modeling efforts. These assumptions help to simplify the logistics
583model without sacrificing our understanding of the problem. Most of
584the assumptions made in this study are influenced by the works of
585Ganeshan (1999), Dasci and Verter (2001), Teo and Shu (2004), and
586Tsao et al. (2012). The assumptions are stated as follows: (i) we
587consider an “arborescent network” where each distribution center can
588serve multiple customer demand points but the opposite is not allowed.
589This assumption is common for any transportation network (e.g.,
590see Melo et al. (2006)). The assumption is even more significant for
591the case of a disaster where a fixed distribution center is assumed to
592serve emergency supplies to multiple customer demand points. (ii)
593demand at each distribution center follows a Poisson process as it is
594generated by the demand originating from customer demand points in
595its influence area (service area); (iii) demand per unit time for customer
596demand points in cluster Ci is an independent and identically
597distributed Poisson process with rate θi (cluster Ci groups customers
598with similar demand). Assumptions (ii) and (iii) are complementary;
599thus, we have explained them together. It is prevalent in literature to
600model demand by Poisson process (see Melo et al. (2011), Araman and
601Caldentey (2011), McCarthy et al. (2008)). Consequently, if demands at
602customer demand points are generated via a Poisson process, demand
603at each distribution center will also be generated via a Poisson process
604as distribution centers will be used to serve the customer demand
605points. (iv) Lateral shipment of products among distribution centers
606and customer demand points is not allowed in the model. All the
607shipments from the distribution centers to the customer demand
608points are via direct shipment. Lateral shipments, through extensively
609used in the supply chain literature, are typically ignored in humanitarian
610logistics literatures due to their impracticability in real life disasterspecific
611applications. Therefore, in this study we assume no lateral
612shipments between the distribution centers and customer demand
613Fig. 1. Service area of the proposed logistics network.
614S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
615171
616points. (v) The Euclidean distance measure is used to calculate the
617distance between a distribution center and a customer demand point.
618For the target area that has a sparse transportation network, we
619understand that Euclidean distance is not a good choice. Therefore,
620we have multiplied the Euclidean distance with a factor (e.g., 1.2) to
621approximate realistic travel distances. (vi) Dasci and Verter (2001) have
622shown that the irregular shapes of the service region (i.e., circular,
623hexagon, square, diamond, triangle) has little effect on the optimal
624solution. Based on this result, we assume that the influence area of each
625distribution center to be almost circular. Moreover, we assume that
626each distribution center is located at the center of the influence areas.
627Note that, in our studyg A r ri
628represents the average distance from a
629distribution center to the demand points in the service region, provided
630that the facility is located at the center of the region. This is a wellknown
631result used by many authors in the CA literature (e.g., Newell,
6321973; Erlenkotter, 1989). These prior studies have experimentally
633validated that the optimal solution has little effect based on the type of
634shape (e.g., circular, hexagon, square, diamond, triangle) considered
635for the service region. (vii) the capacity limitation of the facilities are
636ignored at all network levels; (viii) a sufficient number of drones and
637trucks are available at all times to serve the demand points; (ix)
638customers within an influential area of a distribution center will only be
639served by that specific center and the probability of facility failure is not
640considered; and (x) drones always carry loads with full capacity. The
641last four assumptions are in alignment with a number of such
642assumptions made in prior disaster-specific application studies (e.g.,
643Yushimito et al., 2012; Eksioglu et al., 2009; Moshref-Javadi and Lee,
6442016).
6453.3. Model
646To formulate the cost functions, the entire logistics network is
647expressed in terms of smooth continuous functions. In a two dimensional
648space 2, let the integrated network be represented by a
649continuous service area denoted by R where R ⊆ 2. All the demand
650points occupy a discrete position in this service area R and for each
651cluster Ci, these discrete demand points can be expressed as a spatial
652density slow varying function γ(x) i , where x≔{x , x , …, x } ∈ R 1 2 N .
653Similarly, the customer demand for each cluster Ci can be expressed
654as a spatial density slow varying function θ (x); x ∈ R i . Thus, the
655customer demand at each point x ∈ R can be expressed as a product
656of the customer density and the customer demand, and is given by
657γ(x)θ (x); x ∈ R i i . Let A (x) ri
658be the influence area associated with each
659distribution center in cluster Ci. If we cover the entire area R with
660influence areas of size A (x) ri
661, then the total number of distribution
662centers N (x) ri
663needed at each cluster Ci is given by3
664N x
665R
666A x
667( ) =
668( )
669r
670r
671i
672i (1)
673In our model, three cost functions are considered; namely, total
674facility cost, total transportation cost, and average inventory holding
675cost for the distribution centers. Note the fixed cost of purchasing
676trucks and drones are not considered while calculating the total
677transportation cost. A brief discussion on how to calculate these costs
678are provided as follows:
679• Total facility cost, TF(x): A fixed cost is incurred for establishing and
680running a distribution center. The total facility cost can be obtained
681by multiplying the fixed cost of opening a distribution center fr with
682the total number of distribution centers N (x) ri
683. Note that fr consists
684of opening a refueling station for trucks and a battery charging
685station for drones. Thus, the total facility cost can be expressed as
686follows:
687TF(x) = f N (x) r ri (2)
688• Total transportation cost, TOT(x): It has already been discussed
689that both drones and trucks will be made available to transport
690emergency supplies from multiple distribution centers to demand
691points. Let δ be the length of the planning horizon and let αp and αq
692be the unit transportation cost ($/mile/unit) of transporting emergency
693supplies via truck and drone, respectively. Taking N as the
694total demand points in the disaster affected region, the euclidean
695distance between a distribution center j and a demand point at x ∈ R
696is taken as ∥ x − x ∥ j where j = 1, 2, …, N. Since the impact of
697disaster can not be predicted accurately in advance, let p(x|x) be the
698probability that a demand point x ∈ R can still be served via truck
699after a natural catastrophe (e.g., hurricane, tornado). This implies
700that some roads are still accessible (with a probability p(x|x)) for
701trucks to transport emergency supplies to the customers (e.g.,
702disaster affected zones). In case of the unavailability of road access,
703the demand points x ∈ R can be served by drones with a probability
704q(x|x) where p(x|x) + q(x|x) = 1. Note that, the average travel distance
705from the distribution center, located at the center of the influence
706area, to the customers is proportional to A (x) ri
707. This average distance
708can be denoted by g A (x) r ri
709where gr is a constant scalar. The value
710of gr is approximated by simple geometry given by Li and Ouyang
711(2010) as follows:
712g
713π
714≈
7152
7163 r (3)
717Let [D (x)] ri
718be the expected demand per unit time experienced by
719each distribution center in cluster Ci and
720[D (x)] = γ(x)θ (x)A (x); x ∈ R ri i i ri (4)
721Thus, the total transportation cost, TOT(x) can be expressed as follows:
722????????????? ?????????????
723TOT(x) = p(x|x)α g A δ[D (x)]N (x) + q(x|x)α g A δ[D (x)]N (x) p r r r r q r r r r
724truck transportation cost drone transportation cost
725i i i i i i
726(5)
727Note that the total customer demand during the planning horizon δ
728in whole service area R is given by
729
730∫ γ(x)θ (x)dx i i . Since γ(x)θ (x) i i is a slow
731varying function of x ∈ R, it can be argued that
732
733∫ γ(x)θ (x)dx = δγ(x)θ (x)R i i i i . Therefore, the total transportation cost
734can be rewritten as follows:
735????????????? ???????????
736TOT(x) = p(x|x)α g A δγ(x)θ (x)R + q(x|x)α g A δγ(x)θ (x)R p r r i i qr r i i
737truck transportation cost drone transportation cost
738i i
739(6)
740where R can also be expressed as R = A (x)N (x) ri ri
741.
742Approximation of αp: Determining unit transportation cost
743($/mile/unit) for trucks has been studied in detail in literature.
744Interested readers are encouraged to read works done by Gonzales
745et al. (2013), Barnes and Langworthy (2003), and Marufuzzaman et al.
746(2015).
747Approximation of αq: Drone transportation cost needs to be
748estimated in such a way that resembles drone behavior in a real world
749scenario. In this study, drone transportation cost is estimated by
750considering different stages of a drone flight such as climbing,
751hovering, descending, turning, acceleration and deceleration, rotational,
752and constant speed. Moreover, total flight time of drone is also
753considered in this study. For each of these stages, the amount of battery
754charge depleted and its associated cost is estimated. For calculation of
755service time, we incorporate this component into the hovering time
756calculation. Considering all these factors into equations is nontrivial
7573 Eq. (1) is valid under the assumption that the customer demand is a slow varying
758function of x and thus the influence area of each distribution center can be approximated
759by a circular region as described by Daganzo (1996)
760S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
761172
762and so a non-linear cost function has been developed to represent them
763in a realistic manner.
764As stated, calculating the unit transportation cost, αq ($/mile/unit)
765using drones is not straightforward. A number of factors (e.g., weight,
766carrying capacity, types of power supply, propellers, flight time) govern
767this cost; thus, it becomes very difficult to derive and predict an exact
768energy consumption model. In this study, we represent the energy
769consumption of a drone as a function of its speed and operating
770condition(see Franco and Buttazzo (2015)).
771There are a number of significant factors that contribute to the unit
772cost calculation of drones such as climbing, hovering, descending,
773turning, acceleration and deceleration, rotational, and constant speed
774cost. A pictorial representation of all the costs associated with drones
775are illustrated in Fig. 2. To represent these costs, let us assume that
776each drone can be located in a three-dimensional space (i, j, k) where i
777and j are the component on the X-Y plane and k represents it's altitude.
778Let β be the minimum power needed to hover over the ground and χ be
779a motor speed multiplier. Note that both β and χ are constants that
780depend on the weight of the drone as well as on certain characteristics
781of it's motor/propeller. We further denote Pcl to be the power of the
782motor, v to be the speed of drone while climbing, and th to be the
783hovering time. If Cc is the cost per unit of energy consumed, then
784climbing cost (CC) of drone can be approximated as follows:
785⎛
786⎝ ⎜
787⎞
788⎠ ⎟
789CC P
790k
791v
792= C cl c (7)
793where P ( ) cl
794k
795v
796is the power consumption of the drone that is needed to
797climb at altitude k with speed v. Subsequently, the hovering cost (HC)
798of a drone is given by
799HC = (β + χk)t C h c (8)
800Pugliese et al. (2016) presume that this estimate of the hovering
801cost is almost the same as in real case. Now, it is natural to assume that
802descending speed as well as motor power needed while descending will
803always be lower than the climbing speed and power consumption while
804climbing. Thus, define v′ to be the descending speed where v′ < v and
805Pds to be the power consumption of the drone while descending from
806altitude k with speed v′ where P < P ds cl , the descending cost (DC) can
807now be approximated as follows:
808⎛
809⎝ ⎜
810⎞
811⎠ ⎟
812DC P
813k
814v
815= C
816′
817ds c (9)
818We now calculate the transportation cost of drone due to moving
819from the distribution center to the destination points. We assume that
820the trajectories consist of three phases, i.e., an acceleration phase
821where the drone gains speed from zero, an intermediate phase in which
822the drone has a constant speed, and a final deceleration phase when the
823drone reaches close to its destination point (shown in Fig. 2). Note that
824from the distribution center, which is the center of each Ari
825, the total
826distance covered after these three phases will be on average equal to
827g A (x) r ri
828. Let us divide the overall time to cover g A (x) r ri
829into three
830time periods, i.e., t , t 1 2, and t3. Acceleration occurs to reach speed vc
831from zero during time interval [0, t ] 1 . This constant speed vc is
832maintained by the drone during time interval [t , t ] 1 2 . Finally, the drone
833decelerates to zero again in time interval [t , t ] 2 3 from constant speed vc.
834These variable speed costs(VC) can be approximated as follows:
835⎛
836⎝ ⎜
837⎞
838⎠ ⎟
839VC ∫ P dt ∫ P dt ∫ P dt C
840P T P T P T C
841= + +
842= ( + + )
843t
844acc t
845t
846vc t
847t
848dec c
849acc vc dec c
8500
8511 2 3
8521
8531
8542
8552
8563
857(10)
858where Pacc, Pvc, and Pdec denote the power consumption during
859acceleration, maintaining constant speed, and deceleration, respectively
860and T = t − 0 1 1 , T = t − t 2 2 1 , and T = t − t 3 3 2. Now, the power
861required by the drones can be computed as follows:
862⎛
863⎝ ⎜
864⎞
865⎠ ⎟
866⎛
867⎝ ⎜
868⎞
869⎠ ⎟
870P
871T
872m a d
873T
874=
875Work done by the drone
876=
877* * cc
878(11)
879where m is the total payload a drone can carry including its own weight,
880acc is the acceleration in which drone covers a distance, d is the
881covered distance whose value can be taken as g A (x) r ri
882, and T is the
883total time it takes to cover g A (x) r ri
884with full payload. Using Eq. (11),
885acc can be calculated as follows:
886⎛
887⎝ ⎜
888⎞
889⎠ ⎟a
890P T
891m d
892=
893*
894*
895cc (12)
896As explained before, a drone spends time T1 to reach from zero
897velocity to vc, then flies for time T2 with constant speed vc, and at the
898end spends time T3 during deceleration where it reduces down its
899speed from vc to zero (shown in Fig. 3). From the law of motion we
900know that
901v = u + a *T c cc (13)
902where u is the initial velocity, acc is the acceleration, vc is the final
903velocity, and T is the total time spent by the drone. In this case, the
904initial velocity can be taken as zero i.e., u=0. Thus, Eq. (13) becomes
905v = a *T c cc (14)
906Let us assume that a certain percentage of the drones covered
907region is attributed to accelerating, maintaining constant speed, and
908decelerating the drones. In other words, a drone covers w% of the total
909distance during acceleration phase, b% during constant speed, and the
910remaining y% during deceleration phase (shown in Fig. 3). Thus, using
911Fig. 2. Cost components of drone transportation.
912Fig. 3. Approximate distribution of flight times with distance in drone transportation.
913S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
914173
915Eqs. (11)–(14), we can calculate the time spent at different stages of
916the drone's flight as follows:
917T
918v m w g A x
919P
920=
921* * %* ( ) c r r
922acc
9231
924i
925(15)
926T
927b g A x
928v
929=
930%* ( ) *2 r r
931c
9322
933i
934(16)
935T
936v m y g A x
937P
938=
939* * %* ( ) c r r
940dec
9413
942i
943(17)
944Finally, while traveling from one point to another, motor power is
945also required for providing the rotational movement of the drones. Let
946wturn be the angular rotational speed, pturn be the power needed per
947turn, and nt be the number of turns a drone takes in a given time
948period.With this, the rotational cost (RC) incurred to cover an angle Δθ
949can be computed as follows:
950⎛
951⎝ ⎜
952⎞
953⎠ ⎟
954RC p
955n θ
956w
957= C
958Δ
959turn
960t
961turn
962c (18)
963We are now ready to approximate the unit cost incurred by drone
964αq ($/mile/unit) for the entire traveling distance as follows:
965? ? ?
966???????
967?????
968????????????????? ???????
969? ?
970⎛
971⎝ ⎜
972⎞
973⎠ ⎟
974⎛
975⎝ ⎜
976⎞
977⎠ ⎟
978⎛
979⎝ ⎜
980⎞
981⎠ ⎟
982⎛
983⎝ ⎜
984⎞
985⎠ ⎟
986∫ ∫ ∫
987α CC HC DC VC RC
988P
989k
990v
991C β χktC P
992k
993v
994C
995P dt P dt P dt C p
996n θ
997w
998C
999= + + + +
1000= +( + ) +
1001′
1002+ + + +
1003Δ
1004q
1005cl c h c ds c
1006t
1007acc t
1008t
1009vc t
1010t
1011dec c turn
1012t
1013turn
1014c
1015climbing cost
1016hovering cost
1017descending cost
10180
1019variable speed cost rotational cost
10201
10211
10222
10232
10243
1025(19)
1026Replacing αq in Eq. (6), we obtain the following total transportation
1027cost TOT(x) to serve disaster affected zones using trucks and drones:
1028⎪
1029⎪
1030⎧⎨⎩
1031⎛
1032⎝ ⎜
1033⎞
1034⎠ ⎟
1035⎛
1036⎝ ⎜
1037⎞
1038⎠ ⎟
1039⎛
1040⎝ ⎜
1041⎞
1042⎠ ⎟
1043⎛
1044⎝ ⎜
1045⎞
1046⎠ ⎟
1047⎫⎬⎭
1048∫ ∫ ∫
1049TOT x p x α g A δγ x θ x R C P
1050k
1051v
1052β χkt P
1053k
1054v
1055P dt P dt P dt
1056p
1057n θ
1058w
1059q x g A δγ x θ x R
1060x
1061x
1062( ) = ( | ) ( ) ( ) + + ( + ) +
1063′
1064+ + +
1065+
1066Δ
1067( | ) ( ) ( )
1068p r r i i c cl h ds
1069t
1070acc
1071t
1072t
1073vc
1074t
1075t
1076dec
1077turn
1078t
1079turn
1080r r i i
10810
1082i
1083i
10841
10851
10862
10872
10883
1089(20)
1090For notation convenience, from now on we will use p and q to
1091represent p(x|x) and q(x|x), respectively.
1092• Average inventory cost for distribution center, TOI(x): Each
1093distribution center orders emergency supplies in batches of Q (x) ri
1094and with each batch ordered incurs a reorder cost Rr(x). Therefore,
1095the total reorder cost, TR(x), for all the distribution centers can be
1096expressed as follows:
1097⎛ 
1098⎝
1099⎜⎜
1100⎞
1101⎠
1102⎟⎟
1103TR x N x R x
1104δ D x
1105Q x
1106( ) = ( ) ( )
1107[ ( )]
1108( )
1109r r
1110r
1111r
1112i
1113i
1114i (21)
1115Now, the average inventory at each distribution center is given as
1116the sum of cycle inventory, Q (x)
11172
1118r and safety stock (Graves and Willems,
11192003; Schmidt et al., 2012), z Var[D ] αr ri,LT
1120, where
1121Var[D ] = μVar[D ] + σ [D ] r ,LT r r r r
11222 2
1123i i i
1124and Var[D ] = γ(x)θ (x)A (x) ri i i ri
1125. Here, μr
1126and σr
11272 are the mean and variance of the lead time, LT, respectively and
1128αr is defined as the service level at each distribution center.
1129To calculate the average inventory cost, we further need inventory
1130holding cost at the distribution centers. Let hr be the inventory holding
1131cost of a distribution center for each item over the planning horizon, δ.
1132Therefore, the total average inventory holding cost, TOI(x), can be
1133expressed as follows:
1134??????????????? ???????????
1135⎛ 
1136⎝ ⎜
1137⎞
1138⎠ ⎟
1139⎛
1140⎝
1141⎜⎜
1142⎞
1143⎠
1144⎟⎟
1145TOI x h N x
1146Q
1147z VarD N xR x
1148δ D x
1149Q x
1150( ) = ( )
11512
1152+ [ ] + ( ) ( )
1153[ ( )]
1154( )
1155r r
1156r x
1157α rLT r r
1158r
1159r
1160( )
1161,
1162inventory holding cost reorder cost
1163i
1164i
1165r i i
1166i
1167i
1168(22)
1169With this, we are now ready to calculate the total location-routing
1170cost to serve a disaster affected zone using trucks and drones. Note that
1171all the cost terms derived in this section are in terms of each point x in
1172the service region R. The problem of minimizing the total cost for the
1173entire region R can be expressed as follows:
1174Minimize ∫ (TF(x) + TOT(x) + TOI (x))
1175R (23)
1176subject to
1177N (x)A (x) = R ri ri (24)
1178g A ≤ R r ri d (25)
1179Q (x) ∈ Z r
1180+
1181i (26)
1182A (x) ∈ R r
1183+
1184i (27)
1185where Qri(x)
1186and A (x) ri
1187are the decision variables. Constraints (24) are
1188the area coverage constraints which ensure that the entire service
1189region is covered by the sum of the distribution centers influence area.
1190Constraints (25) indicate that the average distance between the
1191distribution center and customer demand points should never be more
1192than the range of drone Rd. Constraints (26) and (27) guarantee
1193integer and continuous values for Qri(x)
1194and A (x) ri
1195, respectively.
11964. Solution methodology
1197To minimize the total network cost, we have used a two-phase
1198continuous approximation technique developed by Tsao et al. (2012),
1199which in it's essence is a basic extension of the work conducted by
1200Daganzo (1996). The three coastal counties of Mississippi i.e.,
1201Harrison, Hancock, and Jackson counties are used to visualize and
1202validate the modeling results. This is shown in Fig. 4. However, from
1203Fig. 4 it is clear that the customer demand density does not follow a
1204homogeneous Poisson process. As a result, the two-phase continuous
1205approximation method cannot be applied directly to this case since the
1206input function violates the slow varying property which is essential for
1207the analysis using the continuous approximation technique. However, a
1208more detailed analysis of the region considered in the study reveals that
1209there are smaller areas where these functions can be considered as slow
1210varying i.e., the functions are smooth. Thus, to apply the two-phase
1211approximation method in this study, we first divide the entire region
1212into smaller sub-regions over which the functions maintain the slow
1213varying property (phase-I). In phase-II of the solution methodology,
1214the problem is modeled over the sub-regions using the cost functions
1215described in Section 3 and later solved using a continuous approximation
1216approach.
12174.1. Phase-I approximation: a grid cover-couple approach:
1218In phase-I, a grid cover-couple approach is used to divide the entire
1219studied coastal region into sub-regions with slow varying functions. To
1220do this, a mesh of equal sized squares needs to be designed to cover the
1221entire region R. However, while designing this mesh, within each grid,
1222demand must be slow varying. A trial and error method is utilized to
1223determine the feasible shape of the grid. Based on this method, we
1224determine the grid size by taking the variability of demand in the grids
1225into account. Based on our experimentation, we set a 5×5 square mile
1226grid since we observed that within this grid size the demand is found to
1227S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1228174
1229be slow varying. This idea is illustrated in Fig. 4. In doing so, each grid
1230can be attributed a density level that can be easily measured.
1231Now, the grids with similar densities can be clustered together to
1232form areas over which the customer demand point function is slow
1233varying. To partition the region into clusters (C , C , …, C1 2 N), a tolerance
1234limit ϵ is specified that measures the similarity among grids. The
1235tolerance limit defines the amount of variability in the demand density
1236for each grid that is acceptable while treating them as similar. In other
1237words, grids with density at most ϵ apart are considered similar.
1238However, choice of the value of ϵ depends on the customer demand
1239point density pattern of the existing network. Fig. 5 illustrates this
1240process, where in Fig. 5a we divide the three coastal counties of MS into
1241seven clusters based on the ϵ value. After applying this process, it
1242becomes apparent that approximately 70–80% of the population in the
1243three studied counties are concentrated within four of these clusters.
1244Thus, we repeat the same clustering procedure inside these four
1245clusters, however we use a smaller ϵ value to obtain a larger number
1246of clusters. Throughout this process, a total of twelve clusters are made
1247available from these four initial clusters. We proceed by using these
1248twelve clusters and discard the rest. Note, almost all the clusters are
1249located close to the coastal region which is highly prone to natural
1250disasters, and a majority of the people in these counties live within
1251these clusters. Fig. 5b illustrates this process of coupling.
12524.2. Phase-II approximation: distribution center influence area using
1253continuous approximation (CA) approach:
1254Phase-I approximation divides the total service region, R, into a
1255number of clusters (e.g., C , C , …, C1 2 N where N ∈ Z+) with slow varying
1256demand. Within the boundary of each cluster, the Continuous
1257Approximation (CA) technique can be used to formulate the integrated
1258facility location and inventory allocation model for the disaster affected
1259scenarios. The optimization model developed in Section 3 is used for
1260modeling the total logistic costs in each cluster. After solving the
1261optimization model by CA, the solution will give the optimal service
1262region A (x) ri
1263and order quantity Q (x) ri
1264for each distribution center. Note
1265that for each cluster {C}i ∀ i∈N, the size of the influential area A (x) ri
1266is the
1267same, but varies for a different cluster. Moreover, the optimal value of
1268reorder quantity Q (x) ri
1269can also be approximated for each distribution
1270center. Using the influence area of each distribution center, we can
1271easily calculate the total number of distribution centers N (x) ri
1272needed in
1273each cluster {C}i ∀ i∈N. The total number of distribution centers needed
1274for the total affected region can be obtained by summing up all the
1275distribution centers in each and every cluster.
12764.3. Continuous approximation model
1277Each cluster in the service region, R, has a slow varying property.
1278Therefore, we can ignore the dependence of all continuous functions on
1279parameter x. Let Ari
1280be the size of the influence area and Qri
1281be the
1282optimal reorder quantity for each distribution center in cluster {C}i ∀ i∈N.
1283Let TNC(A , Q ) ri ri
1284be the total network cost and is given by the sum of
1285the facility, transportation, and inventory cost functions. The integrated
1286logistics model now becomes:
1287Minimize:
1288⎪
1289⎪
1290⎛
1291⎝
1292⎜⎜
1293⎞
1294⎠
1295⎟⎟
1296⎛
1297⎝ ⎜
1298⎞
1299⎠ ⎟
1300⎧⎨⎩
1301⎛
1302⎝ ⎜
1303⎞
1304⎠ ⎟
1305⎛
1306⎝ ⎜
1307⎞
1308⎠ ⎟
1309⎛
1310⎝ ⎜
1311⎞
1312⎠ ⎟
1313⎛
1314⎝ ⎜
1315⎞
1316⎠ ⎟
1317⎫⎬⎭
1318⎛
1319⎝ ⎜
1320⎞
1321⎠ ⎟
1322⎛
1323⎝
1324⎜⎜
1325⎞
1326⎠
1327⎟⎟
1328⎛
1329⎝
1330⎜⎜
1331⎞
1332⎠
1333⎟⎟
1334∫ ∫ ∫
1335Σ Σ Σ
1336Σ Σ
1337A Q f
1338C
1339A
1340pα g A δγθ C qC P
1341k
1342v
1343β χkt P
1344k
1345v
1346P dt P dt P dt p
1347n θ
1348w
1349g A δγθC
1350h
1351Q
1352z γθμA σ γθA
1353C
1354A
1355R
1356δγθ C
1357Q
1358TNC( , )= + + +( + ) +
1359′
1360+ + + ) +
1361Δ
1362+
13632
1364+ + ( ) +
1365ri ri
1366i
1367N
1368r
1369i
1370ri i
1371N
1372p r ri i i i
1373i
1374N
1375c cl h ds
1376t
1377acc
1378t
1379t
1380vc
1381t
1382t
1383dec turn
1384t
1385turn
1386r ri i i i
1387i
1388N
1389r
1390ri
1391αr i i r ri r i i ri
1392i
1393ri i
1394N
1395r
1396i i i
1397ri
1398=1 =1 =1
13990
14001
14011
14022
14032
14043
1405=1
14062 2
1407=1 (28)
1408subject to
1409g A ≤ R r ri d (29)
1410Q ∈ Z r
1411+
1412i (30)
1413A ∈ R r
1414+
1415i (31)
1416where Ci is the area of cluster i in service region  and
1417⎛
1418⎝ ⎜
1419⎞
1420⎠ ⎟
1421N = r
1422C
1423i A
1424i
1425ri
1426,
1427which denotes the number of distribution centers in cluster i. The
1428decision variables are Ari
1429and Qri
1430for i = 1, 2, …, N.
14314.4. Solution approach
1432As described previously, TNC(A , Q ) ri ri
1433contains two decision variables
1434for each cluster: the optimal influence area of each distribution
1435center, Ari
1436and the ordering quantity for each distribution center, Qri
1437.
1438The optimal values of Ari
1439and Qri
1440minimize the total network cost
1441TNC(A , Q ) ri ri
1442of the system. The objective function of the optimization
1443model, defined by (28)–(31), is a non-linear function containing 2N
1444variables (Ari
1445and Qri
1446for i = 1, 2, …, N). As a result, this model is very
1447challenging to solve. Furthermore, the direct evaluation of the convexity
1448property of the objective function is not straightforward and
1449hence requires a sound methodological approach. To solve this
1450problem, we use the following solution approach that provides optimal
1451values of the decision variables of the problem. The detailed procedure
1452to determine the optimal solution of TNC(A , Q ) ri ri
1453is outlined below.
1454Let us assume that Ari
1455is known. Given Ari
1456, the first and second
1457order derivative of TNC(A , Q ) ri ri
1458with respect to Qri
1459ared Q A
1460dQ
1461TNC( r | ) i ri
1462ri
1463and
1464d Q A
1465dQ
1466TNC( r | ) i ri
1467ri
14682
14692 . We know that if the second derivative of a function is
1470greater than zero, the function is convex. Simple algebraic calculation
1471indicates that
1472Fig. 4. Customer density in grid for Hancock, Harrison, and Jackson counties in MS.
1473S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1474175
1475d Q A
1476dQ
1477R δγθC
1478Q
1479i N
1480TNC( | )
1481=
14822
1483> 0, = 1, 2, …,
1484r r
1485r
1486r i i i
1487r
14882
14892 3
1490i i
1491i i (32)
1492This means TNC(Q |A ) ri ri
1493is a convex function of Qri
1494. The optimal Qri
1495of the model TNC(Q |A ) ri ri
1496can be obtained by equating the first order
1497derivative to zero and solving for Qri
1498. This means that = 0:
1499d Q A
1500dQ
1501TNC( r | ) i ri
1502ri
1503Q A
1504A R δγθ
1505h
1506( )=
15072
1508r r
1509r r i i
1510r
1511i i
1512i
1513(33)
1514Substituting the value of Qri
1515into Eq. (30) yields
1516⎪
1517⎪
1518⎛
1519⎝
1520⎜⎜
1521⎞
1522⎠
1523⎟⎟
1524⎛
1525⎝ ⎜
1526⎞
1527⎠ ⎟
1528⎧⎨⎩
1529⎛
1530⎝ ⎜
1531⎞
1532⎠ ⎟
1533⎛
1534⎝ ⎜
1535⎞
1536⎠ ⎟
1537⎛
1538⎝ ⎜
1539⎞
1540⎠ ⎟
1541⎛
1542⎝ ⎜
1543⎞
1544⎠ ⎟
1545⎫⎬⎭
1546⎛
1547⎝
1548⎜⎜⎜⎜⎜
1549⎞
1550⎠
1551⎟⎟⎟⎟⎟
1552⎛
1553⎝
1554⎜⎜
1555⎞
1556⎠
1557⎟⎟
1558⎛
1559⎝
1560⎜⎜⎜⎜
1561⎞
1562⎠
1563⎟⎟⎟⎟
1564∫ ∫ ∫
1565Σ Σ
1566Σ
1567Σ
1568Σ
1569A Q f
1570C
1571A
1572pα g A δγθC
1573qC P
1574k
1575v
1576β χkt
1577P
1578k
1579v
1580P dt P dt P dt
1581p
1582n θ
1583w
1584g A δγθC
1585h
1586A R δγθ
1587h
1588z γθμA σ γθA
1589C
1590A
1591R
1592δγθC
1593TNC( , )= +
1594+ +( + )
1595+
1596′
1597+ + + )
1598+
1599Δ
1600+
16012
16022
1603+ + ( )
1604+
1605r r
1606i
1607N
1608r
1609i
1610r i
1611N
1612p r r i i i
1613i
1614N
1615c cl h
1616ds
1617t
1618acc
1619t
1620t
1621vc
1622t
1623t
1624dec
1625turn
1626t
1627turn
1628r r i i i
1629i
1630N
1631r
1632r r i i
1633r
1634α i i r r r i i r
1635i
1636r
1637i
1638N
1639r
1640i i i
1641A R δγθ
1642h
1643=1 =1
1644=1
16450
1646=1
16472 2
1648=1
16492
1650i i
1651i
1652i
1653i
1654i
1655r i i
1656i
1657ri r i i
1658r
16591
16601
16612
16622
16633
1664(34)
1665Based on the above discussion, the following algorithm (Algorithm
16661) determines the optimal value of Qri
1667and Ari
1668.
1669Algorithm 1.
1670Step 1. Verify > 0
1671d A
1672A
1673TNC( r ) i
1674ri
16752
16762 to ensure the convexity of the
1677function.
1678Step 2. Solve for = 0
1679d A
1680A
1681TNC( r ) i
1682ri
1683and determine the local minimum
1684points.
1685Step 3. Based on the local minimum points calculate TNC(A ) ri
1686and pick the one with the smallest value.
1687Step 4. Determine the value of Qri
1688, Q (A )≔ ri ri
1689A R δγθ
1690h
16912 ri r i i
1692r
1693by Eq.
1694(33).
1695Step 5. Adjust Qri
1696for i = 1 ∼ N to get the nearest integer values.
1697Step 6. Calculate TNC(A , Q ) ri ri
1698with optimal values of Qri
1699and Ari
1700.
17015. Data description
1702In this section, we provide a detailed description of the data
1703collection procedure undertaken by this study. Historically, the coastal
1704regions of Mississippi are impacted by various types of natural
1705disasters such as hurricanes, tornadoes, and flooding. Among them,
1706hurricanes affect the coastal regions of Mississippi far worse than any
1707other type of disaster. Therefore, in this study we develop the logistics
1708network keeping this disaster type in mind. Three coastal counties of
1709Mississippi have been considered as potential test beds, i.e., Hancock,
1710Harrison, and Jackson counties. Note that whenever a hurricane hits a
1711region, the coastal areas are the one which is usually impacted the
1712most. For instance, when Hurricane Katrina hit in the morning of
1713August 29, 2005 and began its rampage for two days in Mississippi, the
1714most impacted regions were the three coastal counties, although
1715several other counties also experienced significant damage. Keeping
1716that in mind, we collected data from Hancock, Harrison, and Jackson
1717counties for use in this study. Fig. 6 illustrates the effect of Hurricane
1718Katrina in these three counties compared to hurricane Camille that hit
1719Mississippi in 1969. It is observed that Hurricane Katrina had a more
1720devastating effect, because Katrina's surge was higher and covered a
1721wider range as compared to Camille's. The major point here is that in
1722both cases the highly impacted counties in Mississippi were Hancock,
1723Harrison, and Jackson county.
1724Based on the above discussion, we are now ready to fix a set of
1725parameters that can be used as base parameters while determining the
1726optimal values for Ari
1727, Qri
1728, and TNC. We begin by setting the base case
1729values for p and q equal to 0.7 and 0.3, respectively. This implies that in
1730the aftermath of a disaster there is a 70% chance that trucks can be
1731used for relief efforts given the road conditions, while the remaining
1732percentage indicates that trucks are not capable of distributing goods to
1733certain regions and thus drones must be used to serve these regions.
1734We assume that only fuel cost is incurred during transportation via
1735trucks, and this cost increases linearly as more distance is covered. We
1736set αp=$0.40/mile/unit while transporting the emergency commodities
1737via truck (http://www.fedex.com/ratefinder/home). We further
1738set the fixed cost for opening a refueling station for both trucks and
1739drones fr=$140,000 (Rocky Mountain Institute, 2014), Hall (2016)
1740and inventory holding cost hr=$1.0/unit for our base case
1741experimentations.
1742In this study we considered an EHANG-184 drone (http://www.
1743ehang.com/ehang184/) to transport emergency commodities (e.g.,
1744medicines, first-aid, dry food) from multiple distribution centers to
1745the demand points. This drone has a wingspan of 3401 mm. This is one
1746of the largest and heaviest drones available in the market weighing
1747Fig. 5. Original and revised number of clusters for our studied region.
1748S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1749176
1750approximately 200,000 g. In addition to the weight of its built-in
1751components, EHANG-184 can carry up to 99,700 g for a period of
175223 min.4 Note that EHANG-184 can carry more weight, but because of
1753this flight time will reduce significantly. It has a maximum flying
1754altitude (k) of 2.17 miles (Specout, 2016) and contains up to eight
1755motors that put out 142 hp or 106 kW to eight propellers. Charging of
1756its battery pack takes up to four hours in trickle mode and two hours in
1757fast charge mode (Coxworthe, 2016). EHANG-184 can cover up to 10.
175852 miles round trip in each flight with a fully charged battery. This
1759means that the drone can cover up to 5.26 miles of distance one way
1760before returns back to the distribution center (Dronethusiast, 2015). To
1761calculate the drone transportation cost, we set the climbing (v),
1762descending (v′), and constant speed (vc) equal to v=55 m/s, v′=20 m/
1763s, and vc=25 m/s, respectively. The power required for climbing (Pcl),
1764acceleration (Pacc), constant speed (Pcc), deceleration (Pdec), and
1765descending (Pds) by a drone is set to be equal to Pcl=12,000 W,
1766Pacc=90,000 W, Pcc=85,000 W, Pdec=70,000 W, and Pds=5000 W,
1767and respectively. The power needed to hover close to the ground is
1768β=90,000 W while the motor speed multiplier χ is set to 1.0. To
1769calculate the rotating cost, we set the angular rotation speed as
1770wturn=2.1 rad/s, the power needed to turn as pturn=1000 W and
1771the covered angle during turning as Δθ=1.0 rad. Finally, while covering
1772a distance, it is assumed that drones cover 5% of the distance during
1773the acceleration phase (w), 90% of it during constant speed phase (x),
1774and the rest of it during deceleration phase (y). Table 3 summarizes the
1775key input parameters used in this study.
17766. Numerical experiments
1777This section presents the key lessons learned from solving our
1778proposed model TNC(A , Q ) ri ri
1779. We use the three disaster prone counties
1780in Mississippi i.e., Hancock, Harrison, and Jackson counties as a test
1781bed to visualize and validate our modeling results and provide
1782managerial insights under different realistic settings.
17836.1. Base case results
1784We first conduct numerical experiments of the base case scenario to
1785determine the optimal values of service area Ari
1786, ordering quantity Qri
1787,
1788and total network cost TNC. Table 4 presents the computed results
1789after applying the two-phase continuous approximation method proposed
1790in Section 4. Results show that for the base case scenario, we
1791obtain total facility cost TF=$6,716,592, total delivery cost TOT=
1792$10,003,838, total inventory cost TOI=$422,741, and total network
1793cost TNC=$17,143,168 to serve the disaster affected region. Table 4
1794further provides the optimal values of Ari
1795and Qri
1796for the clusters {C} i i∈N
1797depicted in Fig. 5. Based on our observation, clusters with higher
1798demand densities generate lower Ari
1799. Thus, we experimentally validate
1800the finding that the influence area of a distribution center Ari
1801is highly
1802correlated with the demand density of that region. Fig. 7 provides the
1803graphical illustration of TNC versus decision variables Ari
1804and Qri
1805for
1806the twelve clusters depicted in Fig. 5. It is clear from the figure that the
1807two-phase continuous approximation method is capable of determining
1808the optimal system cost TNC under different values of Ari
1809and Qri
1810(Table 4).
18116.2. Impact of drone flying altitude k on disaster relief operation:
1812We now analyze the impact of drone flying altitude k on disaster
1813relief operation. In the United States, Federal Aviation Administration
1814(FAA) requires UAVs to be operated under a ceiling of 400 ft which will
1815severely limit the effective range of utilizing drones for disaster relief
1816operation (Federal Aviation, 2016). However, increasing the flying
1817altitude k for drones will add additional cost in the relief operation.
1818Therefore, we now conduct sensitivity analysis by varying the drone
1819altitude k between 200 and 800 ft and observe its impact on disaster
1820relief operation. Table 5 reports how Ari
1821, Qri
1822, and TNC are impacted by
1823different values of k. Clearly, increasing the k value decreases the
1824average service area Ari
1825provided by drones during the disaster relief
1826operation and we observe the same trend in all the clusters {C} i i∈N
1827Fig. 6. The affected region caused by Hurricane Katrina and Camille (Masters, 2016).
1828Table 3
1829Key input parameters used in this study (Coxworthe, 2016; Specout, 2016).
1830Input Parameter Symbol Value
1831Unit truck transportation cost αp $0.40/mile/unit
1832Motor speed multiplier χ 1.0
1833Climbing power Pcl 12,000 W
1834Descending power Pds 5000 W
1835Acceleration power Pacc 90,000 W
1836Constant speed power Pcc 85,000 W
1837Deceleration power Pdec 70,000 W
1838Turning power pturn 1000 W
1839Hover power β 90,000 W
1840Maximum altitude k 2.17 miles
1841Constant speed vc 55 miles/h
1842Climbing speed v 55 miles/h
1843Descending speed v′ 45 miles/h
1844Angular rotation speed wturn 2.1 rad/s
1845Hovering time th 20 s
18464 This is the only drone available in the market that can carry approximately 100kg in
1847one trip while the closest one is XactSense Titan (http://drones.specout.com/l/144/
1848XactSense-Titan) which can carry only 26kg. This is one of the reasons to consider
1849EHANG type drones in this study.
1850S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1851177
1852depicted in Fig. 5. Since the average service area Ari
1853drops due to the
1854increase in k values, the average relief ordering quantity Qri
1855also drops
1856to serve those reduced areas. This in turn increases the overall system
1857cost TNC for the disaster relief operation Table 7. It is observed that on
1858average the total system cost TNC increases by 4.9% if the drone
1859altitude k increases from 200 to 800 ft.5 Fig. 8 provides a detailed view
1860as how different cost components (e.g., total ordering cost (TOI), total
1861transportation cost (TOT)) vary with the corresponding k values. We
1862see that in Fig. 8(a), as k increases, both TOI and TOT increases. We
1863further show the distribution of TOT cost between drones and trucks in
1864Fig. 8(b). It is observed that as the value of k increases, both the
1865contribution of costs in TOT for drones and trucks also increases. On
1866average, the total transportation cost via drone increases by 20.1% as
1867the value of k increases from 200 to 800 ft. Overall, the Ari
1868and Qri
1869decreases but TNC increases as k increases.
18706.3. Impact of drone speed vc on disaster relief operation:
1871We now study the impact of drone velocity vc on disaster relief
1872operation. Table 6 reports how Ari
1873, Qri
1874, and TNC are impacted by
1875different values of vc. Clearly, if the drone speed vc increases, the
1876average service area Ari
1877and average relief ordering quantity Qri
1878both
1879increase. These results are obvious since drones can now cover more
1880area due to their improved speed capabilities. On average, the drone
1881service area increases by 32.2% if the drone speed vc increases from
188235 miles/h to 75 miles/h. This further increases the average relief
1883ordering quantity Qri
1884by 13.6% since the distribution centers now
1885require additional supplies to support the increased areas of Ari
1886. Since
1887Ari
1888increases, the system now requires fewer number of distribution
1889centers Nri
1890to satisfy the customer demands. Thus, the overall system
1891cost TNC for the disaster relief operation drops as the value of vc
1892increases. Fig. 9 further depicts this relationship along with the cost
1893components TOI and TOT. Note that we observe a significant drop in
1894total transportation cost TOT (especially drone transportation cost) as
1895the value of vc increases from 35 miles/h to 75 miles/h. In summary,
1896we conclude that the disaster relief operation can be conducted in a
1897much cheaper way if the technology of drone speed vc improves in
1898coming years.
18996.4. Impact of road access probability p on disaster relief operation:
1900In this experiment, we analyze the impact of road access probability
1901p on disaster relief operation. Results reported in Table 7 shows that if
1902the road access drops after a disaster (i.e., the value of p decreases),
1903then the value of Ari
1904and Qri
1905decreases but the total network cost TNC
1906increases. The results are obvious due to the fact that drones are still
1907hindered by limited technological capabilities and thus can not carry
1908much weight or fly longer. In other words, the area Ari
1909covered by
1910drones are much smaller than trucks. An example of this can be shown
1911for cluster 1 (i.e., C1 in Table 7) where it is observed that the service
1912area Ar1
1913for drones and trucks are 5.99 and 9.32 miles, respectively.
1914Thus, if the road condition deteriorates after a natural catastrophe,
1915meaning that trucks can no longer serve the demand points, then
1916customer demands will primarily be served by drones at a higher cost.
1917Fig. 10 further illustrates this phenomenon and showed how sensitive
1918different cost components (e.g., TOT, TOI) are against different p
1919values. It is observed that when p=0.0, meaning that all the customer
1920demands will be served by drones, incurs a transportation cost of $14.9
1921million. On the other hand, when we assume that all the road
1922conditions are still accessible even after a massive disaster strikes in
1923Hancock, Harrison, and Jackson counties of Mississippi (i.e., p=1.0),
1924then the overall transportation cost to serve those counties becomes
1925$13.4 million. Therefore, compared to trucks, drone transportation
1926adds an additional $1.5 million to serve the customer demands based
1927on a condition when the road access completely ruins after a massive
1928disaster.
19296.5. Impact of emergency commodities holding cost hr and
1930reordering cost Rr on disaster relief operation
1931In this experiment, we analyze the impact of emergency commodities
1932holding cost hr and reordering cost Rr on disaster relief operation.
1933Holding emergency commodities during disaster costs more than the
1934regular operations. Thus, we first conduct sensitivity analysis on
1935holding cost hr and analyze its impact on disaster relief operation.
1936Table 8 reports how Ari
1937, Qri
1938, and TNC are impacted by different values
1939of hr. Clearly, hr significantly impacts the ordering quantities Qri
1940for all
1941the clusters {C} i i∈N depicted in Fig. 5. It is observed that if hr increases
1942from $1.0/unit to $20.0/unit, then on average Qri
1943drops by approximately
1944346.0% and thus adds an additional 45.1% system cost for the
1945disaster relief operation. In the next experiment, we investigate
1946whether increasing reorder cost Rr has any impact on decision
1947variables and total network cost. We have conducted experiments with
1948five experimental instances where Rr is varied from $3/order to $40/
1949order. Table 9 reports how Ari
1950, Qri
1951, and TNC are impacted by different
1952values of Rr. The results show that when Rr increases, the distribution
1953centers increase their ordering quantities Qri
1954to reduce ordering
1955frequencies. Note that in both instances (i.e., cases with hr and Rr),
1956Table 4
1957Optimal decision variables and cost components for different clusters.
1958Cluster Ari
1959Qri
1960TNC TF TOT TOI NDI
1961Ci (mile2) (units) ($) ($) ($) ($)
19621 7.38 1135 3,366,995 1,400,000 1,885,444 81,551 10
19632 6.45 1194 638,928 280,000 342,774 16,154 2
19643 7.04 1155 4,822,130 1,960,000 2,740,011 122,119 14
19654 6.13 1216 960,646 420,000 515,568 25,078 3
19665 14.33 878 590,461 280,000 301,983 8478 2
19676 10.06 1008 1,659,982 700,000 927,141 32,841 5
19687 10.78 982 1,305,184 560,000 720,770 24,414 4
19698 6.76 1173 1,268,408 560,000 677,446 30,962 4
19709 5.92 1232 1,683,744 700,000 937,154 46,590 5
197110 11.01 974 879,013 420,000 444,177 14,836 3
197211 7.83 1109 883,269 420,000 444,758 18,511 3
197312 27.35 675 207,819 140,000 66,612 1207 1
1974Total 18,266,579 7,840,000 10,003,838 422,741 56
19755 Note that in our calculation we have used the equivalent number in miles since the
1976rest of our units are in miles
1977S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1978178
1979influence area Ari
1980increases slightly to cover more region so as to ensure
1981that a fewer number of distribution centers are needed to minimize the
1982holding and reordering costs for the centers.
19837. Conclusion
1984This paper provides a first comprehensive design framework on
1985how to economically deploy drones to serve a disaster affected region.
1986The study determines the optimal locations for the distribution centers
1987with their corresponding emergency supply inventories and service
1988regions under stochastic demand and drone routing constraints so that
1989the overall system cost (including distribution center opening, inventory,
1990and transportation cost) can be minimized. We approximate
1991drone transportation cost by considering a number of routing specific
1992factors such as climbing, hovering, descending, turning, acceleration
1993and deceleration, rotational, and constant speed cost. The model will
1994Fig. 7. Graphical illustrations of TNC vs. Ari
1995and Qri
1996.
1997S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
1998179
1999provide a tool for the decision makers to use for fast deployment of
2000emergency supplies following a natural catastrophe. This problem is
2001??-hard and deals with densely distributed customer demand points
2002and non-linear drone transportation and emergency supply inventory
2003cost functions. Therefore, we formulate this problem in a continuous
2004model and propose a two-phase continuous approximation approach to
2005solve this model. We decompose the entire studied region into
2006subregions with slow varying functions and propose a continuous
2007approximation approach to determine the optimal influence area,
2008ordering quantity, and system cost to serve the disaster affected region.
2009Table 5
2010Impact of Ari
2011, Qri
2012, NDI, and TNC under different values of k.
2013k (ft) k (f)
2014200 400 600 800 200 400 600 800
2015Cluster
20161 7.99 7.92 7.84 7.38 9 9 9 10
20172 6.79 6.66 6.52 6.46 2 2 2 2
20183 7.39 7.24 7.09 7.04 13 13 14 14
20194 6.52 6.43 6.24 6.13 2 2 3 3
20205 14.79 14.62 14.51 14.33 NDI 2 2 2 2
2021Ari
20226 10.48 10.26 10.12 10.06 4 4 4 5
20237 11.46 11.21 10.94 10.78 3 3 4 4
20248 7.28 7.1 6.94 6.76 3 3 4 4
20259 6.43 6.21 6.04 5.92 5 5 5 5
202610 11.58 11.39 11.24 11.01 2 2 3 3
202711 8.23 8.08 7.96 7.83 2 3 3 3
202812 27.84 27.58 27.42 27.35 1 1 1 1
20291 1186 1174 1161 1135 3.123 3.195 3.21 3.277
20302 1259 11,236 1208 1194 0.531 0.568 0.575 0.585
20313 1243 1221 1202 1155 4.103 4.319 4.592 4.686
20324 1274 1261 1241 1216 0.802 0.823 0.865 0.879
20335 940 921 899 878 0.501 0.513 0.514 0.523
2034Qri
20356 1053 1039 1021 1008 TNC 1.522 1.535 1.586 1.595
20367 1033 1011 997 982 1.215 1.222 1.23 1.242
20378 1204 1191 1186 1173 0.981 0.985 1.001 1.157
20389 1280 1259 1248 1232 1.486 1.502 1.512 1.598
203910 1018 1004 999 974 0.701 0.733 0.748 0.766
204011 1129 1121 1118 1109 0.683 0.705 0.712 0.761
204112 696 689 684 675 0.101 0.102 0.104 0.117
2042Table 6
2043Impact of Ari
2044, Qri
2045, NDI, and TNC under different values of vc.
2046vc vc
2047(miles/h) (miles/h)
204835 45 55 65 75 35 45 55 65 75
2049Cluster
20501 6.3 6.9 7.3 7.7 8.1 11 10 10 9 9
20512 5.5 6.1 6.4 6.7 7 2 2 2 2 2
20523 6 6.6 7.1 7.3 7.7 11 10 10 9 9
20534 5.2 5.7 6.1 6.4 6.6 3 2 2 2 2
20545 12 13.3 14.3 15.1 15.8 NDI 2 2 2 2 1
2055Ari
20566 8.5 9.3 10 10.6 11.1 4 4 4 4 3
20577 9.1 10 10.7 11.3 11.8 3 3 3 3 3
20588 5.8 6.3 6.7 7 7.3 3 3 3 3 3
20599 5 5.5 5.9 6.1 6.4 4 4 4 3 3
206010 9.3 10.2 11 11.6 12.1 2 2 2 2 2
206111 6.6 7.3 7.8 8.2 8.5 2 2 2 2 2
206212 22.6 25.2 27.3 29.1 30.6 1 1 1 1 1
20631 1050 1099 1135 1163 1185 3.605 3.375 3.227 3.123 3.046
20642 1106 1157 1194 1222 1245 0.651 0.611 0.585 0.567 0.554
20653 1069 1119 1155 1183 1206 5.227 4.898 4.686 4.537 4.427
20664 1128 1179 1216 1245 1268 0.977 0.918 0.879 0.853 0.833
20675 804 846 878 903 923 0.597 0.553 0.523 0.503 0.487
2068Qri
20696 928 974 1008 1035 1056 TNC 1.8 1.676 1.595 1.539 1.497
20707 903 948 982 1008 1029 (×106) 1.405 1.306 1.242 1.197 1.163
20718 1086 1136 1173 1201 1223 1.29 1.209 1.157 1.121 1.094
20729 1143 1195 1232 1261 1284 1.773 1.667 1.598 1.55 1.515
207310 895 940 974 1000 1021 0.866 0.805 0.766 0.738 0.717
207411 1025 1074 1109 1137 1159 0.852 0.797 0.761 0.736 0.718
207512 614 648 675 697 714 0.136 0.125 0.117 0.112 0.107
2076S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
2077180
2078Table 7
2079Impact of Ari
2080, Qri
2081, NDI, and TNC under different values of p.
2082p p
20830 0.25 0.5 0.75 1.0 0 0.25 0.5 0.75 1.0
2084Cluster
20851 5.99 6.62 7.17 8.01 9.32 12 11 10 9 8
20862 5.11 5.84 6.29 7.01 7.96 2 2 2 2 1
20873 6.04 6.32 6.85 7.6 8.82 11 11 10 9 8
20884 5.29 5.58 5.98 6.61 7.51 2 2 2 2 2
20895 12.01 12.32 13.74 15.91 20.42 NDI 2 2 2 1 1
2090Ari
20916 8.36 8.86 9.72 10.99 13.41 4 4 4 3 3
20927 9.01 9.42 10.39 11.96 14.55 3 3 3 3 2
20938 5.89 6.11 6.59 7.32 8.42 3 3 3 3 2
20949 5.11 5.42 5.78 6.32 7.21 4 4 4 3 3
209510 9.22 9.62 10.61 11.99 14.92 2 2 2 2 2
209611 6.36 7 7.61 8.56 9.99 2 2 2 2 2
209712 22.05 22.89 25.9 30.99 44.56 1 1 1 1 1
20981 1054 1071 1113 1173 1269 3.305 3.274 3.244 3.194 3.112
20992 1111 1129 1173 1232 1320 0.591 0.588 0.587 0.583 0.573
21003 1076 1089 1134 1192 1286 4.754 4.749 4.706 4.645 4.543
21014 1139 1152 1196 1253 1339 0.882 0.881 0.881 0.877 0.866
21025 801 811 855 919 1042 0.581 0.558 0.535 0.505 0.463
2103Qri
21046 921 942 986 1059 1158 TNC 1.699 1.665 1.616 1.581 1.481
21057 896 911 959 1031 1135 (×106) 1.391 1.298 1.26 1.211 1.143
21068 1086 1111 1151 1212 1301 1.171 1.168 1.162 1.149 1.128
21079 1152 1171 1212 1262 1353 1.591 1.596 1.599 1.594 1.579
210810 876 905 951 1016 1128 0.823 0.801 0.777 0.746 0.702
210911 1028 1042 1088 1189 1247 0.785 0.785 0.767 0.752 0.729
211012 598 611 654 714 857 0.146 0.136 0.123 0.112 0.096
2111Fig. 8. Cost components under different k values.
2112Fig. 9. Cost components under different vc values.
2113S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
2114181
2115We use three disaster prone coastal counties of Mississippi i.e.,
2116Hancock, Harrison, and Jackson counties as a test bed to visualize
2117and validate our modeling results. Managerial insights are drawn to
2118determine how different key drone parameters (e.g., drone flight
2119height, speed), road access availability following a natural catastrophe,
2120and inventory holding and reordering cost impact the overall system
2121performance. Some key insights gained from this study are summarized
2122below:
2123• Influence area Ari
2124of distribution centers to serve a disaster affected
2125region is highly correlated with the demand density of that region.
2126• Increasing the drone flight height k drastically reduces their service
2127area Ari
2128and thus increases the overall system cost TNC for the
2129disaster relief operation (shown in Table 5 and Fig. 8).
2130• Increasing the drone speed vc drastically increases their service area
2131Ari
2132and thus can reduce the overall system cost TNC for the disaster
2133relief operation (shown in Table 6 and Fig. 9).
2134• Unit transportation cost for drones is more expensive than trucks
2135and thus adds an additional cost in the system if the road access p
2136drops significantly after a natural catastrophe (shown in Table 7 and
2137Fig. 10).
2138• Ordering quantity Qri
2139drops with increasing holding cost hr but
2140increases with increasing reorder cost Rr to offset the overall system
2141cost (shown in Tables 8 and 9).
2142This research also opens up a number of future research opportunities.
2143Note that the solution obtained from Algorithm 1 does not
2144guarantee optimality. Future research will devote to this research
2145direction by finding high quality solution approaches for our proposed
2146optimization model. This study ignores wind factors in approximating
2147the routing cost for drones. It will be interesting to see how the
2148inclusion of wind factors in the optimization framework impact the
2149overall delivery performance of drones. Moreover, in real settings when
2150a set of drones will be deployed to serve a disaster affected region then
2151it will be interesting to see how their scheduling and routing mechanisms
2152will be coordinated. In this work, although the model has been
2153Fig. 10. Cost components under different p values.
2154Table 8
2155Impact of Ari
2156, Qri
2157, NDI, and TNC under different values of hr.
2158hr hr
2159($/unit) ($/unit)
21601.00 5.00 10.00 15.00 20.00 1.00 5.00 10.00 15.00 20.00
2161Cluster
21621 7.38 7.39 7.41 7.42 7.42 8 8 8 8 8
21632 6.46 6.46 6.47 6.48 6.49 7 7 7 7 7
21643 7.04 7.05 7.07 7.07 7.08 8 8 8 8 8
21654 6.13 6.14 6.15 6.16 6.17 7 7 7 7 7
21665 14.33 14.36 14.37 14.39 14.4 NDI 15 15 15 15 15
2167Ari
21686 10.06 10.08 10.09 10.11 10.11 11 11 11 11 11
21697 10.78 10.8 10.81 10.82 10.83 11 11 11 11 11
21708 6.76 6.78 6.79 6.8 6.8 7 7 7 7 7
21719 5.92 5.93 5.94 5.95 5.96 6 6 6 6 6
217210 11.01 11.03 11.04 11.05 11.06 12 12 12 12 12
217311 7.83 7.85 7.86 7.87 7.88 8 8 8 8 8
217412 27.35 27.38 27.41 27.42 27.44 28 28 28 28 28
21751 1135 508 360 294 254 3.227 3.533 3.912 4.289 4.665
21762 1194 534 378 309 268 0.585 0.646 0.721 0.796 0.871
21773 1155 517 366 299 259 4.686 5.145 5.713 6.278 6.842
21784 1216 545 385 315 273 0.879 0.974 1.091 1.208 1.324
21795 878 393 278 227 197 0.523 0.555 0.593 0.632 0.671
2180Qri
21816 1008 451 319 261 226 TNC 1.595 1.718 1.869 2.02 2.17
21827 982 439 311 254 220 (×106) 1.242 1.333 1.445 1.557 1.668
21838 1173 525 371 303 263 1.157 1.274 1.418 1.562 1.705
21849 1232 552 390 319 276 1.598 1.774 1.992 2.208 2.425
218510 974 436 308 252 218 0.766 0.821 0.889 0.957 1.024
218611 1109 497 351 287 249 0.761 0.831 0.917 1.002 1.087
218712 675 302 214 175 151 0.117 0.121 0.127 0.132 0.137
2188S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
2189182
2190developed with utilitarian policy, some aspects of egalitarian policy
2191(e.g., fairness and time) can still be satisfied as long as enough trucks
2192and drones are available and the demand points fall within drone
2193range. In future work, we will try to understand how limiting the
2194availability of trucks and drones and incorporating both densely and
2195sparsely population centers will affect the results. We also do not
2196consider battery costs for the replacement right after trips. We assume
2197enough fully charged batteries will be available in the distribution
2198centers for drones. In future studies, we also want to incorporate this
2199factor into a mathematical model.
2200References
2201Abounacer, R., Rekik, M., Renaud, J., 2014. An exact solution approach for multiobjective
2202location? Transportation problem for disaster response. Comput. Oper.
2203Res. 41, 83–93.
2204Afshar, A., Haghani, A., 2012. Modeling integrated supply chain logistics in real-time
2205large-scale disaster relief operations. Socio-Econ. Plan. Sci. 46 (4), 327–338.
2206Ahmadi, M., Seifi, A., Tootooni, B., 2015. A humanitarian logistics model for disaster
2207relief operation considering network failure and standard relief time: a case study on
2208San Francisco district. Transp. Res. Part E: Logist. Transp. Rev. 75, 145–163.
2209Al-Tahir, R., Arthur, M., Davis, D., 2011. Low Cost Aerial Mapping Alternatives for
2210Natural Disasters in the Caribbean. Available from: 〈https://www.fig.net/resources/
2211proceedings/fig_proceedings/fig2011/papers/ts06b/ts06b_altahir_arthur_et_al_
22125153.pdf〉.
2213Araman, V.F., Caldentey, R., 2011. Revenue management with incomplete demand
2214information. Encycl. Oper. Res..
2215Balcik, B., Beamon, B.M., 2008. Facility location in humanitarian relief. Int. J. Logist.
2216Res. Appl. 11 (2), 101–121.
2217Barnes, G., Langworthy, P., 2003. The Per-mile Costs of Operating Automobiles and
2218Trucks. Available from: 〈http://www.conservancy.umn.edu/bitstream/handle/
221911299/909/200319.pdf?Sequence=1isAllowed=y〉.
2220Ben-Tal, A., Chung, B.D., Mandala, S.R., Yao, T., 2011. Robust optimization for
2221emergency logistics planning: risk mitigation in humanitarian relief supply chains.
2222Transp. Res. Part B 45 (8), 1177–1189.
2223Blumenfeld, D.E., Beckmann, M.J., 1985. Use of continuous space modeling to estimate
2224freight distribution costs. Transp. Res. Part A 19 (2), 173–187.
2225Burns, L.D., Hall, R.W., Blumenfeld, D.E., Dganzo, C.F., 1985. Distribution strategies
2226that minimize transportation and inventory costs. Oper. Res. 33 (3), 469–490.
2227Center for Disaster Philanthropy. The disaster life-cycle. Available from: 〈http://
2228disasterphilanthropy.org/the-disaster-life-cycle/〉, 2016.
2229Chen, A., Yu, Y., Ting-Yi, 2016. Network based temporary facility location for the
2230emergency medical services considering the disaster induced demand and the
2231transportation infrastructure in disaster response. Transp. Res. Part B: Methodol.,
223291, 408–423.
2233Cohen, R., 2014. Humanitarian Aid Delivered by Drones: A New Frontier for NGOs?.
2234Available from: 〈https://nonprofitquarterly.org/2014/07/16/humanitarian-aiddelivered-
2235by-drones-a-new-frontier-for-ngos/〉.
2236Coxworthe, B., 2016. Ehang 184 Drone Could Carry You Away One Day. Available from:
2237〈http://www.gizmag.com/ehang-184-aav-passenger-drone/41213/〉.
2238Daganzo, C.F., 1996. Logistics Systems Analysis. Springer, Berlin.
2239Dasci, A., Verter, V., 2001. A continuous model for production-distribution system
2240design. Eur. J. Oper. Res. 129 (2), 287–298.
2241Dronethusiast, 2015. EHang 184 is a Manned UAV You Will Never Get to Fly. Available
2242from: 〈http://www.dronethusiast.com/ehang-184-is-a-manned-uav-you-will-neverget-
2243to-fly/〉.
2244Eksioglu, B., Vural, A.V., Reisman, A., 2009. The vehicle routing problem: a taxonomic
2245review. Comput. Ind. Eng. 57 (4), 1472–1483.
2246Erlebacher, S.J., Meller, R.D., 2000. The interaction of location and inventory in
2247designing distribution systems. IIE Trans. 32 (2), 155–166.
2248Erlenkotter, D., 1989. The general optimal market area model. Ann. Oper. Res. 18 (1),
224943–70.
2250Federal Aviation Administration. Unmanned Aircraft Systems. Available from: 〈https://
2251www.faa.gov/uas/〉, 2016.
2252Franco, C.D., Buttazzo, G., 2015. Energy-aware coverage path planning of UAVs. IEEE
2253International Conference on Autonomous Robot Systems and Competitions.
2254Ganeshan, R., 1999. Managing supply chain inventories: a multiple retailer, one
2255warehouse, multiple supplier model. Int. J. Prod. Econ. 59 (1–3), 341–354.
2256Geoffrion, A.M., 1976. The purpose of mathematical programming is insight not
2257numbers. Interfaces 7 (1), 81–92.
2258Gonzales, D., Searcy, E.M., Eksioglu, S.D., 2013. Cost analysis for high-volume and longhaul
2259transportation of densified biomass feedstock. Transp. Res. Part A 49, 48–61.
2260Graves, S.C., Willems, S.P., 2003. Chapter 3: supply chain design: safety stock placement
2261and supply chain configuration. In: de Kok, A.G., Graves, S.C. (Eds.), Handbooks in
2262OR and MS, 11, 95–132.
2263Greenwood, F., 2015. Above and Beyond: Humanitarian Uses of Drones. Available from:
2264〈http://www.worldpoliticsreview.com/articles/16750/above-and-beyondhumanitarian-
2265uses-of-drones〉.
2266Hall, T., 2016. Starting Your Own Gas Station or Convenience Store in Minnesota.
2267Available from: 〈http://thompsonhall.com/starting-your-own-gas-stationminnesota-
2268business-attorney/〉.
2269Hamedi, M., Haghani, A., Yang, S., 2012. Reliable transportation of humanitarian
2270supplies in disaster response: model and heuristic. Procedia - Social. Behav. Sci. 54,
22711205–1219.
2272Ji, G., Zhu, C., 2012. A study on emergency supply chain and risk based on urgent relief
2273service in disasters. Syst. Eng. Procedia 5, 313–325.
2274Table 9
2275Impact of Ari
2276, Qri
2277, NDI, and TNC under different values of Rr.
2278Rr Rr
2279($/unit) ($/unit)
22800 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
2281Cluster
2282Ari 1 5.99 6.62 7.17 8.01 9.32 12 11 10 9 8
22832 5.11 5.84 6.29 7.01 7.96 2 2 2 2 1
22843 6.04 6.32 6.85 7.6 8.82 11 11 10 9 8
22854 5.29 5.58 5.98 6.61 7.51 2 2 2 2 2
22865 12.01 12.32 13.74 15.91 20.42 NDI 2 2 2 1 1
22876 8.36 8.86 9.72 10.99 13.41 4 4 4 3 3
22887 9.01 9.42 10.39 11.96 14.55 3 3 3 3 2
22898 5.89 6.11 6.59 7.32 8.42 3 3 3 3 2
22909 5.11 5.42 5.78 6.32 7.21 4 4 4 3 3
229110 9.22 9.62 10.61 11.99 14.92 2 2 2 2 2
229211 6.36 7 7.61 8.56 9.99 2 2 2 2 2
229312 22.05 22.89 25.9 30.99 44.56 1 1 1 1 1
2294Qri 1 1054 1071 1113 1173 1269 3.305 3.274 3.244 3.194 3.112
22952 1111 1129 1173 1232 1320 0.591 0.588 0.587 0.583 0.573
22963 1076 1089 1134 1192 1286 4.754 4.749 4.706 4.645 4.543
22974 1139 1152 1196 1253 1339 0.882 0.881 0.881 0.877 0.866
22985 801 811 855 919 1042 0.581 0.558 0.535 0.505 0.463
22996 921 942 986 1059 1158 TNC 1.699 1.665 1.616 1.581 1.481
23007 896 911 959 1031 1135 (×106) 1.391 1.298 1.26 1.211 1.143
23018 1086 1111 1151 1212 1301 1.171 1.168 1.162 1.149 1.128
23029 1152 1171 1212 1262 1353 1.591 1.596 1.599 1.594 1.579
230310 876 905 951 1016 1128 0.823 0.801 0.777 0.746 0.702
230411 1028 1042 1088 1189 1247 0.785 0.785 0.767 0.752 0.729
230512 598 611 654 714 857 0.146 0.136 0.123 0.112 0.096
2306S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
2307183
2308Jia, H., Ordóñez, F., Dessouky, M.M., 2007. Solution approaches for facility location of
2309medical supplies for large-scale emergencies. Comput. Ind. Eng. 52 (2), 257–276.
2310Khayal, D., Pradhananga, R., Pokharel, S., Mutlu, F., 2015. A model for planning
2311locations of temporary distribution facilities for emergency response. Socio-Econ.
2312Plan. Sci. 52, 22–30.
2313Kılcı, F., Kara, B.Y., Bozkaya, B., 2015. Locating temporary shelter areas after an
2314earthquake: a case for turkey. Eur. J. Oper. Res. 243, 323–332.
2315Langevin, A., Mbaraga, P., Campbell, J.F., 1996. Continuous approximation models in
2316freight distribution: an overview. Transp. Res. Part B 30 (3), 163–188.
2317Li, X., Ouyang, Y., 2010. A continuous approximation approach to reliable facility
2318location design under correlated probabilistic disruptions. Transp. Res. Part B 44,
2319535–548.
2320Lin, Y.H., Batta, R., Rogerson, P.A., Blatt, A., Flanigan, M., 2009. Application of a
2321Humanitarian Relief Logistics Model to an Earthquake Disaster. Available from:
2322〈http://www.acsu.buffalo.edu/~batta/TRB_Updated.pdf〉.
2323Luis, E., Dolinskaya, I.S., Smilowitz, K.R., 2012. Disaster relief routing: integrating
2324research and practice. Socio-Econ. Plan. Sci. 46 (1), 88–97.
2325Manopiniwes, W., Irohara, T., 2016. Stochastic optimization model for integrated
2326decisions on relief supply chains: preparedness for disaster response. Int. J. Prod.
2327Res., 1–18.
2328Marufuzzaman, M., Eksioglu, S.D., Hernandez, R., 2015. Truck versus pipeline
2329transportation cost analysis of wastewater sludge. Transp. Res. Part A 74, 14–30.
2330Masters, J., 2016. A Detailed View of the Storm Surge: Comparing Katrina to Camille.
2331Available from: 〈https://www.wunderground.com/hurricane/surgedetails.asp〉.
2332McCarthy, M.L., Zeger, S.L., Ding, R., Aronsky, D., Hoot, N.R., Kelen, G.D., 2008. The
2333challenge of predicting demand for emergency department services. Acad. Emerg.
2334Med. 15 (4), 337–346.
2335Measure-Red cross, 2015. Drones for Disaster Response and Relief Operations. Available
2336from: 〈http://www.issuelab.org/resources/21683/21683.pdf〉.
2337Melo, M.T., Nickel, S., Da, G., Saldanha, F., 2006. Dynamic multi-commodity capacitated
2338facility location: a mathematical modeling framework for strategic supply chain
2339planning. Comput. Oper. Res. 33 (1), 181–208.
2340Melo, M.T., Nickel, S., Da, G., Saldanha, F., 2011. Technical note? Exact analysis of a lost
2341sales model under stuttering poisson demand. Oper. Res. 59 (1), 249–253.
2342Miranda, P.A., Garrido, R.A., 2004. Incorporating inventory control decisions into a
2343strategic distribution network design model with stochastic demand. Transp. Res.
2344Part E 40 (3), 183–207.
2345Moshref-Javadi, M., Lee, S., 2016. The latency location-routing problem. Eur. J. Oper.
2346Res. 255 (2), 604–619.
2347Murat, A., Verter, V., Laporte, G., 2010. A continuous analysis framework for the solution
2348of location? Allocation problems with dense demand. Comput. Oper. Res. 37 (1),
2349123–136.
2350Murphy, R.R., 2014. Disaster Robotics, Intelligent Robotics and Autonomous Agents
2351Series. The MIT Press.
2352Newell, G.F., 1973. Scheduling, location, transportation, and continuum mechanics:
2353some simple approximations to optimization problems. SIAM J. Appl. Math. 25 (3),
2354346–360.
2355Noyan, N., 2012. Risk-averse two-stage stochastic programming with an application to
2356disaster management. Comput. Oper. Res. 39 (3), 541–559.
2357Nozick, L.K., Turnquist, M.A., 1998. Integrating inventory impacts into a fixed-charge
2358model for locating distribution centers. Transp. Res. Part E 34 (3), 173–186.
2359Nozick, L.K., Turnquist, M.A., 2001. Inventory, transportation, service quality and the
2360location of distribution centers. Eur. J. Oper. Res. 129 (2), 362–371.
2361Oran, A., Tan, K.C., Ooi, B.H., Sim, M., Jaillet, P., 2012. Location and routing models for
2362emergency response plans with priorities. Future Secur., 129–140.
2363Ozdamar, L., Ertem, M.A., 2015. Models, solutions and enabling technologies in
2364humanitarian logistics. Eur. J. Oper. Res. 244 (1), 55–65.
2365Pamintuan-Lamorena, M., 2014. Drones Used to Measure Radiation in Fukushima
2366Nuclear Plant. Available from: 〈http://japandailypress.com/drones-used-tomeasure-
2367radiation-in-fukushima-nuclear-plant-2743074/〉.
2368Pugliese, L.D.P., Guerriero, F., Zorbas, D., Razafindralambo, T., 2016. Modelling the
2369mobile target covering problem using flying drones. Optim. Lett. 10 (5), 1021–1052.
2370Pujari, N., Hale, T.S., Huq, F., 2008. A continuous approximation procedure for
2371determining inventory distribution schemes within supply chains. Eur. J. Oper. Res.
2372186 (1), 405–422.
2373Ransikarbum, K., Mason, S.J., 2016. Goal programming-based post-disaster decision
2374making for integrated relief distribution and early-stage network restoration. Int. J.
2375Prod. Econ. 182, 324–341.
2376Rath, S., Gutjahr, W.J., 2014. A math-heuristic for the warehouse location? Routing
2377problem in disaster relief. Comput. Oper. Res. 42, 25–39.
2378Rawls, C.G., Turnquist, M.A., 2010. Pre-positioning of emergency supplies for disaster
2379response. Transp. Res. Part B: Methodol. 44 (4), 521–534.
2380Rennemo, S., Rø, J., Kristina, F., Hvattum, L.M., Tirado, G., 2014. A three-stage
2381stochastic facility routing model for disaster response planning. Transp. Res. Part E:
2382Logist. Transp. Rev. 62, 116–135.
2383Rocky Mountain Institute, 2014. EV Charging Station Infrastructure Costs. Available
2384from: 〈http://cleantechnica.com/2014/05/03/ev-charging-station-infrastructurecosts/〉.
2385Romeijn, H.E., Shu, J., Teo, C., 2007. Designing two-echelon supply networks. Eur. J.
2386Oper. Res. 178 (2), 449–462.
2387Rutten, W.G.M.M., Van Laarhoven, P.J.M., Vos, B., 2001. An extension of the goma
2388model for determining the optimal number of depots. IIE Trans. 33 (11),
23891031–1036.
2390Sahebjamnia, N., Torabi, S.A., Mansouri, S.A., 2015. Integrated business continuity and
2391disaster recovery planning: towards organizational resilience. Eur. J. Oper. Res. 242
2392(1), 261–273.
2393Salman, F.S., Yücel, E., 2015. Emergency facility location under random network
2394damage: insights from the Istanbul case. Comput. Oper. Res. 62, 266–281.
2395Schmidt, M., Hartmann, W., Nyhuis, P., 2012. Simulation based comparison of safetystock
2396calculation methods. CIRP Ann. - Manuf. Technol. 61 (1), 403–406.
2397Shaw, J.M., 2015. The U.S. has More Natural Disasters than any Other Country in the
2398World. Available from: 〈http://www.marketwatch.com/story/the-us-has-morenatural-
2399disasters-than-any-other-country-in-the-world-2015-11-24〉.
2400Shen, Z.M., Coullard, C., Daskin, M.S., 2003. A joint location-inventory model. Transp.
2401Sci. 37 (1), 40–55.
2402Sheu, J., 2007. An emergency logistics distribution approach for quick response to urgent
2403relief demand in disasters. Transp. Res. Part E 43 (6), 687–709.
2404Sheu, J., Lan, L.W., Chen, Y., 2005. A novel model for quick response to disaster relief
2405distribution. Proceedings of the Eastern Asia Society for Transportation Studies.,
24065:2454–2462.
2407Specout, 2016. EHANG184. Available from: 〈http://drones.specout.com/l/328/EHang-
2408184〉.
2409Teo, C., Shu, J., 2004. Warehouse-retailer network design problem. Oper. Res. 52 (3),
2410396–408.
2411Thomas, A., Mizushima, M., 2005. Logistics training: necessity or luxury? Force. Migr.
2412Rev. 22, 60–61.
2413Tofighi, S., Torabi, S.A., Mansouri, S.A., 2016. Humanitarian logistics network design
2414under mixed uncertainty. Eur. J. Oper. Res. 250 (1), 239–250.
2415Tsao, Y., Mangotra, D., Lu, J., Dong, M., 2012. A continuous approximation approach for
2416the integrated facility-inventory allocation problem. Eur. J. Oper. Res. 222 (2),
2417216–228.
2418UNISDR: Center for Research on the Epidemiology of Disasters. The human cost of
2419weather related disasters. Available from: 〈http://www.unisdr.org/2015/docs/
2420climatechange/COP21-WeatherDisastersReport-2015-FINAL.pdf〉, 2015.
2421Wangand, N., Lu, J.C., 2006. Multi-level spatial modeling and decision-making with
2422applications in logistics systems. Technical report, Georgia Institute of Technology.
2423Yan, S., Shih, Y., 2009. Optimal scheduling of emergency roadway repair and subsequent
2424relief distribution. Comput. Oper. Res. 36 (6), 2049–2065.
2425Yi, W., Özdamar, L., 2007. A dynamic logistics coordination model for evacuation and
2426support in disaster response activities. Eur. J. Oper. Res. 179 (3), 1177–1193.
2427Yuan, Y., Wang, D., 2009. Path selection model and algorithm for emergency logistics
2428management. Comput. Ind. Eng. 56 (3), 1081–1094.
2429Yushimito, W.F., Jaller, M., Ukkusuri, S., 2012. A Voronoi-based heuristic algorithm for
2430locating distribution centers in disasters. Netw. Spat. Econ. 12 (1), 21–39.
2431Zhang, J., Li, J., Liu, Z., 2012. Multiple-resource and multiple-depot emergency response
2432problem considering secondary disasters. Expert Syst. Appl. 39 (12), 11066–11071.
2433S. Chowdhury et al. International Journal of Production Economics 188 (2017) 167–184
2434184