Drone-Delivery Network for Opioid Overdose - Nonlinear Integer Queueing-Optimization Models and Methods

The opioid overdose epidemic is ramping up as the number of incidents has quadrupled in the last 15 years. In the midst of this crisis, there are about 130 opioid overdose deaths each day [64]. The first wave of overdoses began in 1990s and was associated to prescription opioids. The 2010 second wave was due to heroin while the third one in 2013 was caused by synthetic opioids, in particular the illicitly manufactured fentanyl [19]. The ongoing COVID-19 pandemic has further exacerbated the opioid health crisis. According to the Centers for Disease Control and Prevention (CDC) [20], over 81,000 drug overdose deaths occurred in the USA between May 2019 and May 2020. Synthetic opioids are the primary driver behind the surge in overdose deaths. The death counts related to the synthetic opioid increased by 38.4% from the 12-month period ending in June 2019 to the 12-month period ending in May 2020. In addition to the substantial human toll, opioid overdose poses a great burden on the economy. The Society of Actuaries [74] estimates that the economic cost of opioid misuse amounts to at least $631 billions from 2015 to 2019 in the USA and includes the cost for healthcare services, premature mortality, criminal justice activity, and related educational and assistance program. The White House Council of Economic Advisers (CEA) estimates the 2015 total cost of opioid overdoses to $504 billion with $431.7 billion due to mortality costs and $72.3 billion due to health, productivity, and criminal justice costs [16].

Many overdose-related deaths are preventable if naloxone, an overdose-reversal medication, can be delivered within minutes following the breath cessation of the overdose victim [18]. Naloxone can be safely administered by a layperson witnessing an overdose incident [77] and its increased utilization has contributed to reducing the mortality rate of opioid overdoses as stressed in [50]. As opioids can slow down or even entirely stop breathing, permanent brain damage can occur within only four minutes of oxygen deprivation. This underscores the need of administering naloxone timely to prevent brain damages or death and explains why naloxone access is of the US Department of Health and Human Services’ three priority areas to combat the opioid crisis [44]. In order to mitigate the opioid crisis, [64] [64] advocate the development of a system that allows 911 dispatchers, trained as drone pilots, to swiftly deliver naloxone with drones to bystanders and to guide them to administer naloxone.

Drones, also known as unmanned aerial vehicles (UAVs), are an emerging technology that can accomplish time-sensitive medical delivery tasks by fastening the time-to-scene response, particularly in the areas lacking established transportation network. A test in Detroit showed that the DJI ’Inspire 2’ drones turned out to be systematically faster than ambulances – the traditional first responders – in delivering a naloxone nasal spray toolkit [75]. In order to make drone deliveries of naloxone a reality and to leverage their full capability, drones need to be strategically deployed and assigned in a timely manner to the randomly arriving overdose-triggered requests (OTRs) to account for the spatial and temporal uncertainty of opioid overdose incidents. The design of drone networks that sets up drone bases, pre-positions drones, and determines their assignment to overdose emergencies remains an open problem in the literature. In that regard, [18] [18] stress the need to develop drone network optimization models for opioid overdoses so that dispatchers “understand when and where to dispatch drones”.

The contributions of this study are fourfold spanning across modeling, reformulation, algorithmic development, and emergency medical service practice.

The remainder of this paper is organized as follows. The literature review of the various fields intersecting with this study is in Section 2. Section 3 derives the closed-form formulas for the queueing metrics used to assess the performance of the drone response network, presents the base formulation of the model, and analyzes its complexity. Section 4 is devoted to the reformulation method while Section 5 describes the algorithms. Section 6 describes the numerical study based on overdose real-life data from Virginia Beach, presents the computational efficiency of the method, and analyzes the healthcare benefits obtained with our model. Section 7 provides concluding remarks. We refer the reader to Appendix A for a list of the acronyms used in this paper and to Appendix B for a list of the mathematical notations.

EMSs struggle to provide a timely response to medical emergencies leading to patients passing away due to the delayed arrival of ground ambulances. Medical drones are viewed as an efficient solution to assist EMS professionals, including for overdoses. They can quickly deliver medicine to patients and take on-sight imagery to support EMS personnel and operations before the arrival of ambulances [66]. We review next the literature on the use of drones for medical emergencies. We focus on overdoses, the topic of this study, and cardiac arrests to which most drone-related studies have been devoted in the medical sphere [59]. We refer the reader to [70] for a recent review on the use of drones in health care.

[79] [79] propose a drone-delivery optimization model for overdoses solved with a genetic algorithm and tested on suspected opioid overdose data from Durham County, NC. The model reduces the response time to 4 minutes 38 seconds by using four drone bases and provides a 64.2% coverage of the county. [36] [36] propose a Markov decision process model that selects which drone should be dispatched to an overdose incident and then relocates the drone. Confronted with the curse of dimensionality, the authors develop a state aggregation heuristic to derive a lookup table policy. A simulation based on EMS data from Indiana reveals that the state aggregation approach outperforms the myopic policy used in practice, in particular when the overall request intensity gets higher, but also highlights the difficulty to estimate the value function. [64] [64] report that all participants in a pilot feasibility study based on 30 simulated opioid overdoses were able to administer the intranasal naloxone medication to the manikin within about two minutes after the 911 call and that 97% of them felt confident that they could do it successfully in a real event.

We now review the drone literature for out-of-hospital cardiac arrests (OHCA). Scott and Scott [72] propose two set covering formulations based on a tandem air-ground response to OHCAs. The models set up the location of the depots (where ground vehicles are stationed) and drone bases without taking into account congestion. Medical supplies are delivered with a ground-based vehicle from a depot to a drone location and the last-mile delivery to the OHCA location is carried out with a drone. Pulver et al. [68] determine where to set up drone bases in order to maximize the number of OHCA requests responded to within one minute. This coverage-type facility location model does not concern with the possible congestion of the network and assumes that the demand and the drone response times are deterministic. Pulver and Wei [67] propose a multi-objective set covering model with backup coverage. Decisions pertaining to the positioning of drones and their assignment to OHCAs are taken to maximize the backup and primary coverage of the network. Boutilier et al. [13] present a two-step approach to design a drone network delivering automated external defibrillators. In the first step, a set covering model is solved approximately to determine the number of drone bases needed to reduce the median response time. In the second step, an $M/M/K$ queueing model determines heuristically the number of drones at each base. The opened stations are fixed (determined in the first stage), and the arrival rate and expected service time are defined as fixed parameters determined ex-ante. Boutilier and Chan [15] propose two coverage-type models that minimize the number of drones needed to satisfy a response time target (i.e., decrease in expected value or conditional value-at-risk of response time). The models do not reward, nor incentivize faster response times. For example, if the target is to decrease the response time by, say, 20 seconds, the models consider a 20-second and an 80-second reduction in the response time as equally satisfactory. The queueing models underlying the formulations assume that the service and response times are independent from the utilization of drones. Mackle et al. [51] present a genetic algorithm to calculate the number of drone bases needed to satisfy a response time threshold. Bogle et al. [10] use a coverage-type facility location model to determine the number and location of drone stations so that a specified number of at-risk OHCA victims are reachable within a certain time. Cheskes et al. [23] use simulation to examine the feasibility and benefits of using drones for OHCAs in rural and remote areas.

The above models assume that the average service time is fixed and independent from the drone assignment decisions. Additionally, they are all covering models in which the objective function is to minimize the total costs or the number of drones (bases) but does not reward a quicker response. Our model contributes to fill in the gaps in the literature. Our strategic network model (i) is a survival network design model as it minimizes the response time which is the primary driver to increase the survival chance of overdosed patients; and (ii) represents the service time as a distribution-free decision-dependent random variable with service rate defined endogenously as a function of the drone dispatching policy.

We explain next how our model differs significantly from the OHCA models proposed in [15] which have, as our model, a queueing underlying structure. First, the objective in [15] to minimize the number of DBs is typical in coverage-type facility location models which, as noted in the literature, are unable to differentiate the impact of varying response times [32], or to reward shorter ones and whose suitability for time-critical medical events has been extensively questioned (see, e.g., [32, 42, 60]). In contrast, our model explicitly rewards quick response time and is hence conducive to maximizing the number of opioid overdose incident survivors. As the overdosed patients’ chance of survival rapidly deteriorates with the time-to-treatment, it is critical to use an objective function which provides a reward diminishing with the time needed to respond. The maximization of the QALY is an alternative objective function that could fulfill the true needs of a medical emergency response network. Second, the underlying queueing models in [15] assume that the service times are exponentially distributed with fixed rates (fixed average service time regardless of drones’ placement) and that the response time between any pair of drone bases and OHCAs incidents is fixed ex-ante and independent from drone utilization and network congestion. Such assumptions do not allow for an accurate representation of the queueing system as explained in a 2023 study [78]. The authors [78] demonstrate that it is critical to track the actual utilization of facilities and servers so that the known impact of utilization can be accurately reflected on the congestion level, waiting time, and response time. Along the same line, it was previously shown [60] that assuming that facilities’ workloads are known a priori and that each facility can provide the same service level irrelevant of their actual utilization significantly hurts the network’s efficiency and reliability. We do not resort to such simplifying assumptions and instead model the opioid overdose arrival rate, the service rate, and the average response times as unknowns dependent on the utilization of the resources and the assignment policy, and define them as decision variables whose values are determined endogenously by the model.

Queueing models can be categorized into two types: descriptive queueing models which conduct an after-the-fact analysis to analyze how a predefined system configuration has performed while optimization-based (prescriptive) queueing models deal with congestion and are used for decision-making purposes (e.g., number of servers, location) [53]. We propose a new optimization-based queueing model that belongs to the class of stochastic network design queueing models with congestion [7]. This model class can be further decomposed into two sub-categories depending on whether the servers are fixed (immobile) requesting the customers to travel to the facility, or mobile and travel to the demand location. Since we consider the delivery of naloxone with drones, our literature review is focused on mobile servers, which is much less extensive than the one for immobile servers (see, e.g., [37]).

The first queueing optimization model with mobile server [9] determines the optimal location of a single facility in order to either maximize the service coverage or to minimize the response cost. Modelling the network as an $M/G/1$ system, the authors consider that, if the server is busy when the demand arrives, the demand can be either rejected (demand is lost and covered by backup servers at a cost) or enter a queue managed under the first-come-first-serve discipline. Building on this study, Chiu and Larson [24] consider a single facility operating as an $M/G/K$ loss queueing system and show that, under a mild assumption, the optimal location reduces to a Hakimi median, which minimizes the average travel time to a client. Batta and Berman [5] study an $M/G/K$ system and develop an an approximation-based approach to minimize the average response time. While the above studies consider a single facility, the queueing probabilistic location set-covering model proposed in [53] considers several facilities, each operating as an $M/M/K$ queueing system. The authors determine the minimum number of facilities needed so that the probability of at least one server being available is greater than or equal to a certain threshold. Later, operating under the auspices of an $M/G/K$ system, Marianov and Revelle [52] study the probabilistic version of the maximal covering problem to site a limited number of emergency vehicles with the goal of maximizing availability when a call arrives. Considering an $M/M/K$ system, [1] [1] determine the number of servers needed to minimize the setup costs to open facilities and operate servers, the travel costs, and the queueing delay costs. In the EMS context, [13] [13] propose a decoupled approach to design a drone network for OHCAs. Each DB operates as an $M/M/K$ system in which the expected service time and the arrival rate are fixed parameters. The coverage-type queueing models in [15] minimize the number of drones needed to satisfy a response time target. They assume that the service rates and response times are fixed and independent from the utilization of the drones, their positioning at DBs, and the congestion in the network. The queueing models in [15] extend the one in [54] for fixed servers and in which the objective is to minimize the number of facilities so that the waiting time or queue size does not exceed a set threshold. An ongoing study [46] considers an $M/G/1$ queueing system for responding to cardiac arrests and proposes new optimality-based bound tightening techniques. In contrast to the above health care studies, this work considers that the service times, the delays, the response times, and the arrival process of (overdose) requests assigned to drones are random variables whose parameters (mean) are endogenized. In particular, the response times and the service rates are decision variables whose values are determined by the model and depend on assignments.

As it will be shown in Section 3.3, the optimization problem studied is a nonlinear fractional integer problem in which the objective function is a ratio of two nonlinear functions each depending on integer variables. The class of problems closest to ours is that of fractional linear 0/1 (binary) problems for which significant improvements have been made over the last years [11, 58]. Such problems arise in multiple applications, such as chemical engineering, clinical trials, facility location [48], product assortment [17], and revenue management (see [12] for a review). We also note a few fractional mixed-integer queueing-optimization problems (see, e.g., [37, 31, 39]. Such problems pose serious computational challenges due to the pseudo-convexity and the combinatorial nature of the fractional objective. They minimize a single or a sum of ratios in which the denominator and numerator are linear functions of binary variables. When more than one ratio term is in the objective function, the problem, even unconstrained, is NP-hard [65].

One approach to tackle fractional linear 0-1 programs is to move the fractional terms to the constraint set and to then derive MILP reformulations, which requires the introduction of continuous auxiliary variables and big-M constraints (see, e.g., [47]). However, MILP solvers have difficulties – in particular when the number of ratio terms increases – due to the significant lifting of the decision and constraint spaces and the looseness of the continuous relaxation induced by the big-M constraints. Within this family, the MILP reformulation proposed in [11] and based on the binary expansion of the integer-valued expressions in the ratio terms contains many less bilinear terms and hence less linearization variables and constraints. This formulation scales much better but can generate weak continuous relaxations, which hurts the convergence of the branch-and-bound algorithm. Mixed-integer second-order cone reformulations have also been proposed (e.g., [73]). They are based on the submodularity concept and the derivation of extended polymatroid cuts [3] to obtain tighter continuous relaxations. However, it is reported [58] that state-of-the-art optimization solvers still struggle to solve moderate-size mixed-integer conic problems and that their performance degrades quickly as size increases. Building upon the links between MILP and mixed-integer conic reformulations, [58] [58] derive tight convex continuous relaxations while limiting the lifting of the decision and constraint spaces.

The above methods, while providing very valuable insights, can not be directly used to the problem tackled here as the latter differs from the fractional 0-1 problem along several dimensions. First, the denominator and the numerator of each ratio term are nonconvex functions and each involve binary as well as general integer decision variables. Second, the objective function is the sum of a very large number (i.e., several thousands) of ratio terms. Third, the integer variables are multiplied by fractional parameters which prohibits the use of the MILP method (see [11]). Fourth, the nonlinear terms are not simply bilinear. The numerator of each ratio term includes exponential terms and higher-degree polynomial terms while each denominator includes polynomial and factorial terms.

The problem studied in this paper is called the Drone Network Design Problem for Opioid Overdose Response (DNDP). Consider an EMS provider that seeks to design a drone network to deliver naloxone, i.e., an opioid antagonist, to an opioid overdose incident with the objective to minimize the response time and thereby to maximize the chance of survival of the overdosed patient. Given a set of candidate locations (e.g., fire, police, EMS stations), some of them are selected to be set up as DBs where the available medical drones can be deployed. The DNDP model is a bilocation-allocation problem that simultaneously determines the locations of the DBs, the capacity of those (i.e., number of drones at each DB), and the response policy determined by the assignment of drones to OTRs in order to minimize the response time. The response time is defined as the sum of the queueing delay time experienced if drones are busy when requested for service and the drone flight time from a DB to an OTR location. The network is capacitated to account for the severely limited availability of medical resources. Considering the opioid epidemic, Buckland et al. [18] stress that the amount of medical resources (drones, ambulances, paramedics) is “substantially smaller than the number of patients who require immediate care”.

The following notational set is used in the formulation of the model (see Appendix B). Let $I$ be the set of OTR locations and $J$ be the set of potential DB locations. Let $A_{j}$ be a vector of parameters $(a_{j},b_{j},c_{j})$ representing the coordinates of DB $j$ in the earth-centered, earth-fixed coordinate system while $A_{i}$ is a vector of parameters $(a_{i},b_{i},c_{i})$ representing the location $i$ of an OTR. The parameter $d_{ij}=\|A_{i}-A_{j}\|$ is the euclidean distance between DB $j$ and OTR location $i$. A number $p$ of drones with speed $v$ can be deployed at up to $q$ $(q\leq p)$ open DBs across the network. Due to battery and autonomy limitations, each drone has a limited coverage defined by its catchment area with radius $r$. A drone can only service OTRs that are within its catchment radius [13], since the drone must have battery coverage to return to its base. Accordingly, we define the sets $J_{i},i\in I$ (resp., $I_{j},j\in J$) which include all the DBs (resp., OTR locations) that are within $r$ of OTR location $i$ (resp, DB $j$). The binary decision variable $x_{j}$ is 1 if a DB is set up at location $j$, and is 0 otherwise. The binary variable $y_{ij}$ takes value 1 if an OTR at $i$ is assigned to DB $j$ and value 0 otherwise. The general integer decision variable $K_{j}$ is the number of drones deployed at DB $j$. By convention, we denote the upper and lower bound of any decision variable $x$ by $\overline{x}$ and $\underline{x}$, respectively.

The stochastic nature of the occurrence of overdoses and the resulting uncertainty about the arrival rates of OTRs at DBs along with the uncertain service times can cause delays, requests being queued, and waiting times until a drone can be dispatched. In order to account for those, we model the network as a collection of $M/G/K_{j}$ queues in which each DB $j$ operates as an $M/G/K_{j}$ queueing system where the capacity $K_{j}$ of DB $j$ is a bounded general integer decision variable representing the number of drones to be deployed at DB $j$. The occurrence of opioid overdoses at any location $i$ follows a Poisson distribution with arrival rate $\lambda_{i}$. The service time of drones is a random variable with general distribution and with known first and second moments. The arrival rate $\eta_{j}$ of OTRs at DB $j$ can be inferred to be a Poisson process since it is a linear combination of independent Poisson variables (i.e., weighted sum of assigned OTRs; see (5)). The arrival rate of OTRs at DBs is unknown ex-ante and is defined endogenously via the solution of the optimization problem. The same applies to the drone service time whose expected value is also endogenized. It follows that the arrival process, i.e., demand at DBs and the service times of drones are endogenous sources of uncertainty with decision-dependent parameter (expected value) uncertainty as coined in [40]. If a drone cannot be dispatched on the spot after reception of an OTR, the request is placed in a queue depleted in a first-come-first-served manner. After providing service, drones travel back to the DB to be cleaned, recharged, and prepared for the next trip [64].

We now derive steady-state closed-form expressions for the queueing metrics - response time, queueing delay, service time – needed to formulate the DNDP problem. We use the notations $S_{j}$, $R_{i}$, and $Q_{j}$ to represent the random variables respectively denoting the total service time at DB $j$, the response time for an overdose at location $i$, and the queueing delay at DB $j$.

Different from immobile servers, the travel time to and from the scene are included in the service time for mobile servers [8]. The service time for any OTR at $i$ serviced by DB $j$ [9] is $S_{ij}=\beta\frac{d_{ij}}{v}+\alpha_{i}+\epsilon_{i}$, where $\beta$ is a constant that allows for different travel speeds, and $\alpha_{i}$ and $\epsilon_{i}$ are independent and identically distributed (i.i.d) random variables representing the on-scene service time and the drone’s reset time (to be charged and prepared for the next demand), respectively. Let $\xi_{i}=\alpha_{i}+\epsilon_{i}$, with expected value $\mathbb{E}[\xi_{i}]$. The expected value $\mathbb{E}[S_{ij}]$ of the service time conditional to DB $j$ responding to an OTR at $i$ is:

$\displaystyle\mathbb{E}[S_{ij}]$

$\displaystyle=\beta\frac{d_{ij}}{v}+\mathbb{E}[\xi_{i}]\ .$

(1)

The expected value of the total response time for $j$ is the sum of the expected service times for all the OTRs serviced (i.e., $y_{ij}=1$) by the drones stationed at $j$ and depends on the dispatching variables $y_{ij}$

$\displaystyle\mathbb{E}[S_{j}]$

$\displaystyle=\frac{\sum_{i\in I_{j}}\lambda_{i}y_{ij}\mathbb{E}[S_{ij}]}{\sum_{i\in I_{j}}\lambda_{i}y_{ij}}$

(2)

while the second moment of the total service time at DB $j$ is:

$\displaystyle\mathbb{E}[S_{j}^{2}]$

$\displaystyle=\frac{\sum_{i\in I_{j}}\lambda_{i}y_{ij}\mathbb{E}[S_{ij}^{2}]}{\sum_{i\in I_{j}}\lambda_{i}y_{ij}}\ .$

(3)

The service time $S_{j}$ follows a general unspecified distribution with known first and second moments and the opioid overdoses occur at each location $i$ according to a Poisson process with rate $\lambda_{i}$. Accordingly, each DB $j$ is modelled as an $M/G/K_{j}$ queueing system in which a variable and upper bounded number $K_{j}$ of drones can be stationed and whose expected queueing delay can be approximated [63, 69] as

$\displaystyle\mathbb{E}[Q_{j}]\approx\frac{\eta_{j}^{K_{j}}\mathbb{E}[S_{j}^{2}]\mathbb{E}[S_{j}]^{K_{j}-1}}{2(K_{j}-1)!(K_{j}-\eta_{j}\mathbb{E}[S_{j}])^{2}\big{[}\sum_{n=0}^{K_{j}-1}\frac{(\eta_{j}\mathbb{E}[S_{j}])^{n}}{n!}+\frac{(\eta_{j}\mathbb{E}[S_{j}])^{K_{j}}}{(K_{j}-\eta_{j}\mathbb{E}[S_{j}]}\big{]}}$

(4)

which represents the expected waiting time until a drone is ready to be dispatched after reception of an OTR. The arrival rate of OTRs at DB $j$ is given by

$$\eta_{j}=\sum_{i\in I_{j}}\lambda_{i}y_{ij}\ ,\ j\in J\vspace{-0.105in}$$

(5)

which shows that the arrival rate is unknown ex-ante and depends on the assignment decisions $y_{ij}$, thereby highlighting the decision-dependent parameter uncertainty of the OTR arrivals at DBs.

The expected response time in steady-state for an OTR at location $i$ is:

$\displaystyle\mathbb{E}[R_{i}]=\sum_{j\in J_{i}}\left(\mathbb{E}[Q_{j}]+\frac{d_{ij}}{v}\right)y_{ij}\ .\vspace{-0.107in}$

(6)

We present now the base formulation of the DNDP problem, which belongs to the family of nonlinear fractional integer problems, and analyze its complexity. Before presenting its formulation, we recall the most distinctive features of the proposed DNDP model. First, the minimization of the average response time is a survival objective function as it contributes to increasing the survival chance of the victims of opioid overdoses. The chance of survival is a monotone decreasing function of the response time (see, e.g., [4, 27]). By minimizing response time and rewarding quick care, our model effectively maximizes the chance of survival of overdosed victims, which is viewed in the health care literature as a survival optimization model. Second, the queueing-optimization model accounts for congestion and decision-dependent uncertainty as the service and arrival rates characterizing the Poisson arrival process and service times of OTRs at DBs are determined endogenously. Third, the capacity – number of drones positioned of each $M/G/K_{j}$ DB $j$ is not fixed ex-ante, but is a bounded decision variable.

The fractional nonlinear integer base formulation $\mathbf{B-IFP}$ of problem DNDP is:

	$\displaystyle\mathbf{B-IFP:}\min$	$\displaystyle\;\frac{1}{\sum_{l\in I}\lambda_{l}}\Bigg{[}\sum_{i\in I}\sum_{j\in J_{i}}\frac{\lambda_{i}d_{ij}y_{ij}}{v}+$
		$\displaystyle\sum_{i\in I}\sum_{j\in J_{i}}\frac{\lambda_{i}y_{ij}\sum_{l\in I_{j}}(\lambda_{l}y_{lj}\mathbb{E}[S_{lj}^{2}])(\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])^{K_{j}-1}}{2(K_{j}-1)!(K_{j}-\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])^{2}\big{[}\sum_{n=0}^{K_{j}-1}\frac{\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])^{n}}{n!}+\frac{(\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])^{K_{j}}}{(K_{j}-1)!(K_{j}-\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])}\big{]}}\Bigg{]}$		(8a)
	$\displaystyle s.to\ \$	$\displaystyle\sum_{i\in I_{j}}\lambda_{i}y_{ij}<K_{j}\frac{\sum_{i\in I_{j}}\lambda_{i}y_{ij}}{\sum_{i\in I_{j}}\lambda_{i}y_{ij}\mathbb{E}[S_{ij}]},\ \ j\in J$		(8b)
		$\displaystyle\sum_{j\in J_{i}}y_{ij}=1,\ \ i\in I$		(8c)
		$\displaystyle y_{ij}\leq x_{j}\ \ i\in I,j\in J_{i}$		(8d)
		$\displaystyle\sum_{j\in J}x_{j}\leq q$		(8e)
		$\displaystyle x_{j}\leq K_{j}\leq Mx_{j},\ j\in J$		(8f)
		$\displaystyle\sum_{j\in J}K_{j}=p$		(8g)
		$\displaystyle x_{j},y_{ij}\in\{0,1\},\ \ i\in I_{j},j\in J$		(8h)
		$\displaystyle K_{j}\in\mathbb{Z}_{+},\ \ j\in J$		(8i)

The objective function (8a) minimizes the average response time of the network, a key feature given the time-critical nature of opioid overdoses. Given that the queueing formulas used to represent the response time assume that the system is in steady-state, we introduce the nonlinear constraints (8b) to ensure the stability of the queueing system [57] as they prevent the arrival rate at any DB from being equal or larger than the service rate. Constraint (8c) makes sure that each OTR is serviced by exactly one drone while (8d) requires that an OTR can only be assigned to an open DB within the drone catchment area. Constraint (8e) limits from above the number of open DBs with the parameter $q$. Constraint (8f) ensures that drones can only be placed at an opened DB and that an opened DB must have at least one drone. The parameter $M$ defines the maximum number of drones that can be stationed at any DB $j$. As suggested in [61] and as implemented in [6, 67, 36] for medical drone networks, we set $M$ equal to two as drones are scarce resources and are typically spread to expand the coverage of the network. Additionally, our preliminary tests based on real-life data show that giving the possibility to place a third drone at a DB does not reduce the response times. Constraint (8g) ensures that all available drones are deployed. The binary and general integer nature of the decision variables are enforced by (8h) and (8i) with $\mathbb{Z}_{+}$ denoting the set of nonnegative integer numbers. Remark 2 summarizes several key features of the model.

The validity of the above statements is demonstrated in Appendix C.2. We also propose in Appendix D an illustration for a small network and present the formulation of the objective function in $\mathbf{B-IFP}$.

The derivation of an MILP reformulation is relatively complex and we split it into two main steps (see Proposition 2 and Theorem 3) to ease the exposition.

Proposition 2 proposes an equivalent reformulation taking the form of a fractional nonlinear binary problem R-BFP with nonconvex continuous relaxation. The proof is given in Appendix C.3.

The following comments are worth noting. A first difference with $\mathbf{B-IFP}$ is that $\mathbf{R-BFP}$ has only binary decision variables and does not include any general integer variables $K_{j}$. The second difference and advantage of $\mathbf{R-BFP}$ over $\mathbf{B-IFP}$ is that the polynomial terms in the reformulation $\mathbf{R-BFP}$ are of lower degree than those in $\mathbf{B-IFP}$. Third, in $\mathbf{R-BFP}$, there is no decision variable $K_{j}$ appearing in the upper limits of the summation operations, such as $\sum_{n=0}^{K_{j}-1}(\sum_{l\in I_{j}}\lambda_{l}y_{lj}\mathbb{E}[S_{lj}])^{n}$ in $\mathbf{B-IFP}$. Fourth, there is no exponential term of form $y_{ij}^{K_{j}}$ in $\mathbf{R-BFP}$. Fifth, the factorial terms $(K_{j}-1)!$ in $\mathbf{B-IFP}$ which, besides being nonlinear, can also be undefined if $K_{j}=0$, are no longer present in $\mathbf{R-BFP}$ in which they are replaced by the fixed parameters $(m-1)!$. The undefined issue is resolved since $m$ is at least equal to 1. The reformulation $\mathbf{R-BFP}$ is valid for any value assigned to $M$.

The constraints (8c)-(8e), (8h), and (9b)-(9f) representing the feasible area of problem R-BFP define collectively a linear binary feasible set which, to ease the notation, is thereafter referred to as:

$$\mathcal{B}=\left\{(x,y,\gamma)\in\{0,1\}^{|J|+\sum_{i\in I}|J_{i}|+M|J|}:\eqref{assignment}-\eqref{eq_q};\eqref{binary};\eqref{steady-state-2}-\eqref{binary2}\right\}\ .\vspace{-0.025in}$$

(10)

We give in Appendix D the objective function of problem $\mathbf{R-BFP}$ for a small network. Table 9 in Appendix E specifies the number of constraints and integer decision variables of each type in the base formulation $\mathbf{B-IFP}$ and its reformulation $\mathbf{R-BFP}$.

Although simpler than $\mathbf{B-IFP}$, $\mathbf{R-BFP}$ is still a challenging optimization problem as it involves polynomial and fractional terms and its continuous relaxation is nonconvex. We shall now demonstrate in Theorem 3 that the MILP problem $\mathbf{R-MILP}$ is equivalent to $\mathbf{R-BFP}$ and $\mathbf{B-IFP}$. The proof is split into two main steps that starts with the moving of the fractional terms from the objective function into the constraint set and continues with the linearization of the polynomial terms. To shorten the notations in the proof, we replace $\mathbf{E}[S_{lj}]$ and $\mathbf{E}[S_{lj}^{2}]$ by $\tilde{S}_{lj}$ and $\tilde{S}^{2}_{lj}$, respectively. We first consider the medical drone network design problem DNDP in which $M=2$ (i.e., up to two servers per DB). Theorem 4 will later show that an MILP reformulation can be derived for any value of the parameter $M$, or, in other words, for any $M/G/K$ queueing system in which the capacity of each system is unknown and variable.

We recall two well-known linearization methods for bilinear terms used in the proof of Theorem 3:

• •

  Product of binary variables [35] :
  The polynomial set $\{(y_{1},y_{2},z)\in\{0,1\}^{2}\times[0,1]:z=y_{1}y_{2}\}$ is equivalent to the MILP set:

  $$\mathcal{F}_{y}^{z}=\Big{\{}(y_{1},y_{2},z):z\geq 0,\;z\geq y_{1}+y_{2}-1,\;z\leq y_{i},\;i=1,2\Big{\}}\ .\vspace{-0.1205in}$$

  (11)

• •

  Product of a bounded continuous variable by a binary one [56]:
  Let $0\leq\underline{x}<\bar{x}$. The polynomial set $\{(y,x,z)\in\{0,1\}\times[\underline{x},\bar{x}]\times[0,\bar{x}]:z=yx\}$ can be equivalently represented with the MILP set:

  $$\mathcal{M}_{yx}^{z}=\Big{\{}(y,x,z):z\geq\underline{x}y,\;z\geq\bar{x}(y-1)+x,\;z\leq\bar{x}y,\>z\leq\underline{x}(y-1)+x\Big{\}}\ .$$

  (12)

The proof of Theorem 3 is given in Appendix C.4 and is decomposed into two main parts. The linearization part first demonstrates that the nonlinear problems B-IFP and R-BFP are MILP-representable. The compaction part shows that the number of linearization constraints and variables can be significantly reduced. This result is attained by first showing that $\omega^{(2)}_{ij}=\mu^{(2)}_{ij}=\mu^{(2)}_{ij}\cdot\gamma^{(2)}_{j},i\in I,j\in J_{i}$, which is in turn used to prove that the two sets of linear inequalities $(\gamma^{(1)}_{j},\mu^{(1)}_{ij},\omega^{(1)}_{ij})\in\mathcal{M}_{\ \gamma\mu}^{{}^{\prime}\omega},\ i\in I,j\in J_{i}$ (13h) and

$$(\gamma^{(m)}_{j},\mu^{(m)}_{ij},\omega^{(m)}_{ij})\in\mathcal{M}_{\gamma\mu}^{\omega}\;,\;i\in I,j\in J_{i},m=1,\ldots,M\vspace{-0.075in}$$

(14)

are equivalent (i.e., the superscript in (14) is $(m),m=1,\ldots,m$ while it is $(1)$ in (13h)). In Section 6.4, we assess the computational benefits of the compaction approach by comparing the results with the compact MILP reformulation R-MILP and the larger-size MILP reformulation L-MILP

$${\bf L-MILP}:\max\eqref{obj_lin}:\eqref{COM1}-\eqref{COM4};\eqref{OLD-R-MILP}\vspace{-0.1in}$$

(15)

in which (14) replaces (13h) in R-MILP. Observe that the MILP reformulation L-MILP contains $\sum_{i\in I}|J_{i}|$ more decision variables and $4\ \sum_{i\in I}|J_{i}|$ more linearization constraints than R-MILP.

We shall now demonstrate in Theorem 4 that the linearization approach proposed above for an $M/G/K$ queueing system with variable number of servers (capacity) $K$ limited from above to 2 can be extended to any $M/G/K$ system with finite variable capacity $K$.

The proof of Theorem 4 is given in Appendix C.5 and shows that the objective function and each constraint (48) can be linearized regardless of the value of $M$. Note that T1, T2, and T3 in (16b) are polynomial expressions. More precisely, T1 and T2 involve polynomials of degree $(m+2)$ with monomials involving the product of one single continuous bounded variable by up to $(m+1)$ binary variables while T3 includes polynomials of degree $(m+1)$ with monomials defined as products of $(m+1)$ binary variables. Additionally, the polynomials in T1, T2, and T3 have a nested structure (see, e.g., [25, 34]) for which customized valid inequalities could be derived to strengthen the formulation.

A lazy constraint is an integral part of the actual constraint set and is unlikely to be binding at the optimum. Instead of incorporating all lazy constraints in the formulation, they are grouped in a pool and are at first removed from the constraint set before being (possibly) iteratively and selectively reinstated on an as-needed basis. The targeted benefit is to obtain a reduced-size relaxation or outer approximation problem, which is quicker to solve. One should err on the side of caution when deciding which constraints are set up as lazy. Indeed, the inspection of whether a lazy constraint is violated is carried out each time a new incumbent solution for the outer approximation problem is found and the computational overheads consecutive to the reintroduction of violated lazy constraints can be significant.

Within this approach, the reduced-size relaxation is solved at each node of the tree. Each time a new incumbent solution is found, a verification is made to check whether any lazy constraint is violated. If it is the case, the incumbent integer solution is discarded and the violated lazy constraints are (re)introduced in the constraint set of all unprocessed nodes of the tree, thereby cutting off the current solution.

We define here the linearization constraints in the sets $\mathcal{F}_{y}^{z}$ and $\mathcal{M}^{\tau}_{zU}$ as lazy constraints. The motivation for this choice is threefold and guided by: the results of preliminary numerical tests, the very large number $\sum_{j\in J}|I_{j}|\cdot(|I_{j}|-1)/2$ of such constraints, and the fact that the inequalities in the set $\mathcal{M}^{\tau}_{zU}$ are not always needed, i.e., they only play a role if two drones are placed at the same DB.

We introduce next the notations used when using the reformulation $\mathbf{R-MILP}$ within the algorithm. The exact same procedure is applicable for $\mathbf{L-MILP}$ and only requires changing the name of the feasible set. Let $\mathcal{O}$ denote the set of open nodes in the tree. We call iteration a node of the branch-and-bound tree at which the optimal solution of the continuous relaxation is integer-feasible. Let $\mathcal{C}$ be the entire constraint set of problem $\mathbf{R-MILP}$ (see Theorem 3), $\mathcal{L}_{k}$ be the set of lazy constraints at node $k$, $\mathcal{V}^{L}_{k}$ be the set of violated lazy constraints at $k$, and $\mathcal{A}_{k}:=C\setminus\mathcal{L}_{k}$ be the set of active constraints at $k$, i.e., the set of constraints of the reduced-size outer approximation problem $\mathbf{OA-MILP}_{k}$. The composition of the sets vary across the algorithmic process.

The outer approximation framework is designed as follows. At the root node ($k=0$), we have:

	$\displaystyle\mathcal{L}_{0}:=\{\mathcal{F}_{y}^{z}\cup\mathcal{M}^{\tau}_{zU}\}.$		(17)
	$\displaystyle\mathcal{A}_{0}:=\{\mathcal{B};\eqref{U}-\eqref{V};\mathcal{M}^{\mu}_{yU};\mathcal{M}^{\omega}_{\gamma\mu}\}.$		(18)
	$\displaystyle\mathcal{V}_{0}^{L}:=\emptyset.$		(19)

The linearization approach for R-MILP involves the introduction of big-M constraints which typically lead to loose continuous relaxations. To tighten the continuous relaxation, we derive new valid inequalities and optimality cuts.

A valid inequality does not rule out any feasible integer solutions but cuts off fractional solutions feasible for the continuous relaxation problem. In essence, it pares away at the space between the linear and integer hulls, thereby providing a tighter formulation. In contrast to lazy constraints, a valid inequality is inserted in the formulation if the optimal solution of the continuous relaxation at any node is fractional and violates this valid inequality. We shall now derive several types of valid inequalities.

The valid inequalities (20) reflect the fact that, if either one of the variables $U_{j}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{(1)}}}$ or $U_{j}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{(2)}}}$ related to the expected delay is positive, then at least one drone is positioned at DB $j$, and that, if DB $j$ is not active, then the variables $U_{j}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{(m)}}},m=1,2$ are equal to 0.

The next valid inequalities (21) reflect that the queueing delay at a DB $j$ is a decreasing function of the number of drones positioned at this DB.

The proofs of Proposition 5 and Proposition 6 are given in Appendix C.6 and Appendix C.7.

Besides valid inequalities, we also derive optimality cuts (see Proposition 7) which cut off integer feasible solutions that are not optimal or in a way that not all optimal solutions are removed. The proposed optimality cuts (22) state that the opening of a DB at location $j$ is required if and only if either one or two drones are positioned at $j$. The proof of Proposition 7 can be found in Appendix C.8.

The incorporation of the valid inequalities and optimality cuts in the outer approximation branch-and-cut algorithm OA-B&C has a significant computational impact as shown in Section 6.4.1. The outer approximation branch-and-cut algorithm OA-B&C is structured as follows.

First, the set $\mathcal{B}$ of optimality cuts (22) are added to the pool of lazy constraints (17):

$$\mathcal{L}^{\prime}_{0}:=\mathcal{L}_{0}\cup\mathcal{B}.\vspace{-0.05in}$$

(23)

Second, the valid inequalities are derived up-front, incorporated into a pool of user cuts, and applied and checked dynamically each time the nodal optimal solution $X^{*}_{k}$ is fractional through a user cut callback implemented in Gurobi. Let $\mathcal{U}_{0}$ be the set of user cuts (valid inequalities) at the root node: $\mathcal{U}_{0}:=\{\eqref{VI3}-\eqref{VI4}\}$ while $\mathcal{V}^{U}_{k}$ is the set of valid inequalities violated by the fractional optimal solution of the continuous relaxation at node $k$.

The valid inequalities modifies the outer approximation framework as follows. If $X^{*}_{k}$ is fractional, the user callback is applied, and the violated – if any – valid inequalities in the current user cut pool are added to the active constraint set of each open node $o\in\mathcal{O}$ (i.e., to cut off $X^{*}_{k}$) and removed from $\mathcal{U}_{o}$:

$$\mathcal{U}_{o}\leftarrow\mathcal{U}_{o}\setminus\mathcal{V}^{U}_{k}\quad\text{and}\quad\mathcal{A}_{o}\leftarrow\mathcal{A}_{o}\cup\mathcal{V}^{U}_{k}.\vspace{-0.085in}$$

If no inequality in $\mathcal{U}_{o}$ is violated by $X^{*}_{k}$, then two branching constraints are entered to cut off $X^{*}_{k}$ and the next open node is processed. Note that, if $X^{*}_{k}$ is integer feasible, the user cut callback is not applied. The algorithm stops when the set of unprocessed nodes becomes empty. The pseudo-code of the algorithmic method OA-B&C is in Appendix F.

The dataset used in the tests describes all the opioid overdose requests (actual and suspected opioid incidents) in the city of Virginia Beach and is publicly available [26]. The data were collected through multiple sources, including the OpenVB data portal, Freedom of Information Act requests to the government of Virginia Beach, and public reports, and were validated by the Virginia Beach officials.

The dataset contains all dispatch records for OTRs from the second quarter of 2018 to the third quarter of 2019, which amounts to a total of 733 data points (overdoses). Each record has four fields including the time at which the request was received by the EMS, the response time, i.e, time between the reception of the request and the arrival of the EMS personnel on the scene, and the location (i.e., latitude and longitude) of the request. The dataset also provides the location of the twenty-six established EMS facilities (i.e., fire, police, and EMS stations) that can be selected as drone bases. As in [13], we consider that drones can travel at a speed of $27.8$ meters per second (m/s) and take 10 seconds to take off and to land. We use 25 minutes as the expected non-travel service time which includes the time to administer naloxone at the overdose location and to recharge and prepare the drone for the next assignment.

We have conducted a grid search in order to identify the minimal size of the drone network to be able to respond to all OTRs. This reveals that the drone network should include $p=11$ drones and does not require more than $q=10$ DBs. Thereafter, we refer to base scenario the case in which the drone response network allows for the opening of ten DBs and the deployment of eleven drones.

In this section, we study the relationships between response time, survival chance, delivery mode, and compare the drone-network constructed with our approach with the ambulance-based EMS network currently in use in Virginia Beach. To assess the two EMS systems under similar experimental settings, we use the historical real-life data for the ambulance-based EMS system in Virginia Beach and the data obtained from the simulator (presented next in Section 6.2.1) for the drone system. Using the results calculated numerically with our optimization model may have provided a less balanced comparison as the mathematical model may ignore some practical issues that may exist in realty. Section 6.2.2 analyzes the reduction in the response time attributable to the drone network, the applicability and robustness of the proposed approach, and the spatio-temporal adjustments in the OTRs and drone network. Section 6.2.3 investigates the increase in the chance of survival of an overdose victim and the expected number of lives saved. Section 6.2.4 quantifies the impact on the QALY and performs a cost analysis.