key: cord-0648336-qm5rznxi
authors: Ong, Joshua; Kulpanowski, David; Xie, Yangxinyu; Nikolova, Evdokia; Tran, Ngoc Mai
title: OpenEMS: an open-source Package for Two-Stage Stochastic and Robust Optimization for Ambulance Location and Routing with Applications to Austin-Travis County EMS Data
date: 2022-01-26
journal: nan
DOI: nan
sha: 30d922076372ce7832017149fe23921818638810
doc_id: 648336
cord_uid: qm5rznxi

Emergency Medical Systems (EMS) provide crucial pre-hospital care and transportation. Faster EMS response time provides quicker pre-hospital care and thus increases survival rate. We reduce response time by providing optimal ambulance stationing and routing decisions by solving two stage stochastic and robust linear programs. Although operational research on ambulance systems is decades old, there is little open-source code and consistency in simulations. We begin to bridge this gap by publishing OpenEMS, in collaboration with the Austin-Travis County EMS (ATCEMS) in Texas, an end-to-end pipeline to optimize ambulance strategic decisions. It includes data handling, optimization, and a calibrated simulation. We hope this open source framework will foster future research with and for EMS. Finally, we provide a detailed case study on the city of Austin, Texas. We find that optimal stationing would increase response time by 88.02 seconds. Further, we design optimal strategies in the case where Austin EMS must permanently add or remove one ambulance from their fleet.

The Emergency Medical Service (EMS) provides pre-hospital treatment and transportation. Receiving fast pre-hospital treatment is vital in scenarios such as cardiac arrest, stroke, and other severe trauma cases [29] . It is generally agreed upon that these calls need pre-hospital treatment within 8-10 minutes [36, 12, 7, 25] . If a call is not answered with in this threshold, it is considered a shortfall. To help EMS services achieve the best response time the operational research community has widely studied how ambulance stationing and routing can optimize response times. [7, 6, 27, 34, 11, 35, 39, 32, 22, 24] .

Strategic decisions can also help EMS departments save money. White et al. [41] found health care spending has dramatically increased in the US compared to other countries. Further, Heightman [20] approximates that it costs half a million dollars a year to staff one ambulance on a 24/7 basis. One way to combat increased spending is to use current EMS resources more efficiently. Optimization could help a city with a decreasing budget maintain performance while removing an ambulance from its operation.

The operational research community commonly constructs the ambulance problem into a two-stage optimization problem [7, 6, 27] . The first stage is the ambulance location problem (ALP). This answers the question: where should we station our ambulances around the city to achieve an optimal coverage? Ambulances can be stationed around the city in locations such as hospitals, fire stations, ect. In the ALP we want to maximize the coverage of the ambulances to respond to calls, using historic call data as empirical probability distributions.

There have been many case studies on ambulance stationing and routing [34, 11, 35, 39, 32, 22, 24] , but the methodologies regarding discrete-event simulations vary widely. Some are theoretical and include no simulations, some use private softwares, while only a few provide open-source tools [7] . In this work, we develop OpenEMS, a modular software with opensource code consisting of four well-separated parts: data handling, optimization, numerical computations of discrete time events, and visualization. This allows future research that innovates on one or several components to be directly compared with previously published work. We hope that this increased transparency would facilitate the translation from research results to field implementations by EMS departments. Our main contributions are as follows:

(1) A software pipeline at https://github.com/joshua-ong/ AmbulanceDeployment that can be applied to any city. The software includes: discretizing a set area into a grid, calibration, validation, optimization using state-of-the-art (SOTA) methods from Bertsimas et al. [7] , simulation, and visualizing improvements compared to historical performance.

(2) A case study of Austin EMS, including 246,809 thousand calls within 1039 square miles, demonstrating the impact of optimal decision making.

The rest of the paper is organized as follows. Section 2 gives a brief overview of ambulance dispatch research. Next, we include the necessary theory of the optimization models in Section 3. We detail the simulation model in Section 4. Next, we have a case study of the city of Austin in section 5. Next, we include a socio-economic of Austin in section 6. Next, there is a user code sample in section 7. We conclude in Section 8.

There exist many case studies across various countries showing that EMS systems' heuristics are not optimal. Morohosi [34] , Eaton et al. [14] , and Burkey et al. [11] use popular coverage linear programs to study Tokyo, Japan; Austin, Texas; and 4 states in the southern USA respectively. Eaton et al. [14] , and Burkey et al. [11] also use coverage to examine social equality. Nilsang et al. [35] included a case study in Bangkok, Thailand where they utilize social media to design a covering model. They found they could reduce the number of ambulance locations by 81.6 % without reducing coverage. Swalehe et al. [39] included a case study in Eskisehir Province, Turkey using GIS (Geogrpahic Information System) methods and found they could reduce the response time from 10 minutes to 5 minutes for 77.6 % of ambulance demand areas. McCormack et al. [32] included a case study in London, England using genetic algorithms and testing the algorithm on a calibrated and validated simulation model. Doener et al. [13] includes a broader case study of Austria including coverage lin-ear programs for ambulance systems but also other related areas like inventory routing and disaster relief. Nonetheless, none of these papers include open-source code.

These case studies can also answer important hypothetical questions. McCormack et al. [32] used simulations to find how performance would decrease if there were resource gaps that day. Ingolfsson et al. [22] tested the policy of a 'single start station', where all ambulances start at the same station using the ProModel software [19] . Although the policy of having all ambulances start at the same station seems inefficient, it has the benefits of having pooled resources in one location. Inakawa et al. [24] used simulation on the city of Seto, Japan to find that adding a new location was twice as effective as adding two ambulance stations. Aboueljinane et al. [3] used simulation in Val-de-Marne department of France using the ARENA software [18] . In it they recommended repositioning teams to increase coverage and found that it improved response time by 3.5%.

Implementing novel algorithms into an EMS department has proven to be difficult. For example, in 1992, the London Ambulance Service (LAS) tried to implement a computerized system including a Computerized Dispatch System (CAD) (rather than human operators), and automatic vehicle stationing [15] . According to Adamu et al. [15] the software had a £ 7.5 million budget, but only lasted 9 days before it was closed due to failure and increased mortality. The software was abandoned after 9 days as implementations did not have time to learn and fail, because there are human lives at stake. It further shows the significance of testing novel algorithms on confident simulations.

A notable success began with Mason and Henderson's [21] simulation model BartSim, which later matured into a company [30] called Optima. Optima currently has case studies in New Zealand, Australia, Denmark, United Kingdoms, Canada, and the United States. However, Optima is a commercial software designed for EMS operators. In contrast, OpenEMS is open-source and designed to support researchers making the latest advances in the ambulance allocation and routing problem with data validations while incurring minimal technical setup and no charge.

Some of the earliest work began when ReVelle et al. [37] proposed the Maximum Available Location Problem (MALP), inspired by the more general Maximum Coverage Problem (MCP). In MALP they optimize the location of stations for ambulances to maximize the coverage, where a caller is covered if there is an ambulance that can reach the call within a fixed time. This is called the ambulance location problem (ALP). Later Batta et al. [5] introduced Maximum Expected Covering Location Problem (MEXCLP), which used an ambulance's "busy fraction" as a Bernoulli random variable to approximate the expected coverage. Other classical models include the Double Standard Model (DSM) by Gendreau et al. [16] . These remain the important benchmarks in many case studies [7, 34, 14, 11] .

Further research from the operational research community has since been applied to solving the ambulance problems. Beraldi et al. [6] solved the ALP with chance constraints to model routing ambulances to an unknown demand. For more information on two stage stochastic formulation, we refer the reader to [4] and [9] . A two-stage robust formulation was first implemented by Bertsimas et al. [7] who applied a data-driven model to generate scenarios. A two-stage distributional robust formulation was implemented by Liu et al. [27] .

Other studies integrate a tool called the approximate dynamic programming (ADP). Some of the first papers on applying ADP to the ambulance problem were from Maxwell et al. [31] and Schmid et al. [38] . Later Jagtenberg et al. [23] used ADP to show that the closest ambulance routing policy is not necessarily optimal. For a more comprehensive overview, we refer readers to the review by Aboueljinane et al. [2] and Tassone et al. [40] .

In this section, we lay out the theoretical formulation of two-stage stochastic and robust formulations appeared in Bertsimas et al. [7] to solve ALP and ARP. We adopt their notation in this paper. In a two-stage optimization problem, there is the first stage here-and-now decision which must be made under uncertainty and the wait-and-see variable that happens afterwards.

In our context, at the first stage we must select how many ambulances are stationed at each location. Given the set of ambulance stations I; let x i be the number of ambulances at location i ∈ I. This is the first stage decision, as we must station the ambulances before we know the realization of the demand. In the second stage the calls are revealed and we must determine the optimal routing of stationed ambulances to these calls. Let J be the set of regions, then each call can be assigned to a region j ∈ J. Let y ij be the routing of ambulance stationed at i to region j. Let τ be the longest time it takes for an ambulance to serve a call and return to their station. Then we can group calls within periods of τ where we assume an ambulance can serve one call in each period. A complete table of notation is included in table 1.

Set of ambulance locations J Set of demand regions E Set of feasible edges from I to J M Set of observed ambulance calls n Total number of ambulances τ

Total time (minutes) for an ambulance to serve a call and be available again y ij Number of ambulances routed from location i to region j x i Number of ambulances at location i d j

Number of ambulance calls from region j B I Incidence matrix representing incoming edges B J Incidence matrix representing outgoing edges Table 1 .: Summary of notation for ambulance location and dispatch.

The stochastic formulation objective is to minimize the expectation of shortfalls. To solve this problem we use Monte-Carlo sampling to construct a discrete distributionP(d) where we assume each sample of calls is equally likely. The stochastic formulation is written as:

subject to e x ≤ n (1b)

The objective 1 contains a slack variable that is bounded from below later. Constraint 1b limits the number of ambulances stationed to be less than or equal to the total number of ambulances. Constraint 1c iterates through each Monte-Carlo sample (m ∈ M ), for each sample we limit a stationed ambulance to only be sent once per time period, τ . Constraint 1d makes the slack variable greater than or equal to the shortfalls, which is equal to calls minus calls answered. (Note this can be thought of as a maxflow problem of a bipartite graph, where x and d are the vertices of the bipartite graph and y are the edge values.)

In the robust optimization formulation, we want to minimize the worst-case scenarios. To do this we first model demand as a Poisson distribution, γ. (Single is the demand of a single grid cell. Local is the demand of adjacent grid cells. Regional is the demand within 10 minutes of a grid cell. Global is the demand of all of Austin.) We then construct the uncertainty set D(α) as the α Value-at-Risk of these Poisson distributions. Finally, we take the worst scenarios from and minimize over them in the robust formulation like:

where:

The formulation is the same as the stochastic formulation except for the objective. For more information we refer the reader to the original formulation by Bertsimas et al. [7] on how to dualize 2 and make it tractable with the column-constraint method. Zeng and Zhao [42] first designed the column-constrained method for two-stage robust optimization. The column-constraint method has the same idea as the Benders-style but with faster heuristic implementation.

In this section, we describe the different components of the simulation model in the general setting. This includes the steps to replicate this on a different case study. Notably, we measure the confidence of the simulation response time so that an EMS department can understand the possible variance in real-world applications. First, we explain the state space of a general Emergency Medical System (EMS) and the correspondence of these events to the simulation. Next,we explain how the simulation represents the real world by discretizing it into a grid. We then verify the model by simulating historical actions and comparing them with what happened to validate that the model works in normal conditions. Finally, we explain the process of tuning the linear program parameters for best results.

A general framework of the EMS procedure for answering an emergency call is detailed in figure 1 . The first state includes from the time the call was made to when it was assigned to an ambulance, where an operator must collect enough information like the seriousness of the injury and location. The time the ambulance is en route to when it arrives at the call is labeled the response time. This is arguably the most important metric for survival rate of the patient. From the time the ambulance arrived to when it departs the scene captures the time it took to administer pre-hospital treatment and load the patient. The time departed from the scene to the hospital is drive time. The ambulance is now available to serve another call, if there are no new calls, it will return to its station. Our simulation model includes the same steps as the state space just described. Each step is an event that is put into a priority queue with respect to time. An overview of the code is included in pseudo-code 1. In the beginning, all calls are input into the priority queue at their respective times. The calls then propagate the other events following the state space. We do not include the setup time, because it is only a fraction of the whole ambulance travel time and is not impacted by our strategic decision. The travel times are determined using Open Street Map, which was introduced by Haklay and Weber [17] . The intervention time is modeled as a lognormal distribution. They have been shown to effectively fit ambulance processes like in Lam et al. [26] and McLay and Mayorga [33] . 

A grid construction is necessary for both solving the linear programs and modeling the simulation. The main idea is to discretize space by representing the continuous space with a grid, like in figure A3 . If an event happens somewhere within a grid cell, we assume it happens at the center of the grid cell. Discretizing the ambulance space is common in both linear programming theory [7] as well as other simulation models [32] .

To construct the grid, first we inscribe the case study region into a rectangle. Then we split the grid into uniformly sized rectangle members. Note there is a tradeoff between how much resolution (number of grid cells) to include and how much computation is required. If the grid is finer it has more accurate travel times but longer computations and vice versa. Once the grid is established, any call within a rectangle member will be considered to be at the center. Then the travel time between any two points becomes the distance of the members of the grid cells they belong to.

We compute travel times with Open Street Map (OSM). In general, we find computed travel times are slower than the reported travel time. Brown et al. [10] found that ambulances using lights and sirens traveled on average 105 seconds faster than ambulances without. This is to say that on average the special privileges of an ambulance allow it to go faster by bypassing traffic. Thus we recommend adjusting grid times to be more accurate. Regression will depend on the data available and other apriori knowledge, but can be relatively simple including temporal variables like day-of-week and hour or spatial variables like the number of grid cells away the path is.

The verification step builds confidence that the model reflects the real world. It evaluates how well the calibration step compares to real world measurements. The verification model inputs historical actions and compares the simulated times with the reported times. Because there is a lot of randomness in the drive time to do with paths and traffic, we compare the overall trends for n = 1000 calls or about two weeks. We then look at the mean response time error on average. To measure confidence, we then invoke Multiple Replication Procedure (MRP) on multiple batches of n-calls and measure the variance among batches.

As suggested by Bertsimas et al. [7] , the α hyperparameter is determined experimentally for the robust linear program. This is done using cross-validation. In step one we choose a list of alpha values, for example [0.1, .01, .001]. In step two we split the dataset into a k-fold of training and testing sets with an 80/20 split. We optimize the robust optimization on the training set and simulate it with our model on the testing set. We then remove any alpha values that are saturated (every station includes exactly one ambulance) and choose the alpha value which gives the best performance.

In this section, we provide a case study of Austin, Texas, as an application of our open-source OpenEMS. Austin-Travis County EMS (ATCEMS) 1 has 37 full-time ambulances that on average serve 338 calls per day. Data in this study includes 246,809 calls from the 2019-2020 period. ATCEMS serves 2.2 million residents within a 1,039 square mile space. Since Austin is one of the biggest cities in the United States, it is an interesting and representative subject to study.

Here we give an outline of what publicly available data exists. The main data entries included were the date of event, longitude, and latitude of call, and reported ambulance response time. (There also exists other entries such as reported ambulance assign time, reported ambulance on-site time, reported ambulance hospital time, ect.) We then graphed all calls and constructed the temporal distribution in figure 2 . From this figure, we noticed that the most consistently high demand calls are from 8:00-20:00 and from Monday-Friday. We call these times peak hours because they have more consistently high call demand than the offpeak hours. Further, we remove Saturday and Sunday because the times of high demand are less consistent. For the following results, we limit our experiments to these peak hours as they are the times which Austin EMS system is most stressed. Offpeak hours tend to be slower and easier to manage. We also plotted the data geospatially in 3. Like most cities, a majority of calls happen in the downtown region. Austin also includes a sprawling suburb area that has relatively fewer calls but takes a long time to drive between. This makes Austin hard to solve using heuristic planning to handle both the high-demand area of downtown while still being able to handle the occasional suburban calls. 

The calibration step is meant to calibrate the grid time to the reported travel time. To do this we tried several regression models to fit our grid time with the reported travel time. In general, we find that there is a low correlation because while grid time is deterministic, the reported travel time is dependent on many unknown factors such as traffic. We tried including several other variables like time of day and the spatial distance from two points but we found there were marginal benefits.

For this case study, we found that trimming the top and bottom 1 percent of data was necessary to avoid outlier reported travel times like 0 minutes (which is impossible) or calls that took more than an hour. Then we fitted the log of grid time to the log of reported time. The r 2 value was .1472. One possible reason is that the variance is mostly not dependent on the average grid time, but on unobserved factors such as traffic. We computed the adjusted grid times and compared them to the reported grid times in figure 4 . From this we observe that the adjusted times are squeezed with shorter tailends than before. From the validation we argue this produces a fair mean response time, but loses some resolution of extreme events as well as skewing times around 10 minutes which makes it harder to calculate shortfalls. 

The results of the verification are shown in figure 5 . Each data point represents the mean response time of 1000 calls or approximately 2 weeks' worth of data. Before the grid time had an error of −186.23 ± 70.13 seconds, while with the adjusted grid time the error becomes 2.41 ± 35.31. From this, we see the calibration reduces the error. Thus we could expect the simulation mean response time to be close to the actual events within a quantified confidence. Figure 5 .: This is a plot of 20 batches where each batch is one data point.

Next, we determine the best α value experimentally with a 3-fold cross-validation for the robust optimization. Results are in the appendix A1. The set of optimal stationing was then tested on the same hold-out set. Interestingly, a majority of the optimizations have the same optimal value. This is the case except for when α becomes too big (.1) where it may not predict as many scenarios and when α becomes too small (.0001) where it may predict for too extreme scenarios. From this cross-validation we choose α = .01 for the following experiments.

The code was primarily written in Julia [8] with auxiliary scripts in Python. "Julia For Mathematical Programming" (JuMP) [28] was used to interface Gurobi 9.0 with Julia. The computations are done on a standard computer. The station has 16 GB of RAM and a 2.5 GHz processor.

To compare with real data we use 40 ambulances, similar to Austin EMS during the day. Simulations are run for 1000 calls. Then 12 batches of 1000 call simulations are run to establish confidence using Multiple Replication Procedure (MRP). (Note if a simulation is run on too few samples it is biased to perform better because all of the ambulances are available in the beginning. Then it will reach an equilibrium. If shortfalls become too numerous in a certain region, they can plug the system and give a bias towards worse performance.) The time it takes to serve is based on a lognormal distribution with a mean 3.65 and variance 0.3. Calls are filtered to peak hours to be more homogenous (see subsection 5.2). In this simulation, the closest ambulance is always sent to the call. All parameters are summarized in table 2.

Value number of ambulances 40 number of calls 1000 number of batches 12 serve time distribution log-normal(µ = 3.65, σ = 0.3) α .01 Table 2 .: Summary of experimental parameters.

Results from the simulation can be summarized in table 3. We found the stochastic optimization performed better than robust optimization. Although, the robust optimization sometimes outperformed the stochastic optimization when there were more ambulances. This may be because robust optimization is more conservative, waiting for an anticipated scenario that does not happen in these batches. Both results perform better than the reported call time though. Namely, the stochastic optimization is expected to have a 88.02 second faster mean response time than heuristic decision making. The ambulance stationing was visualized in figure 6 . The purple circles are proportional to how many ambulances should be stationed at each location. We note that the robust opti-mization will concentrate more ambulances in one station. We believe this is to handle some anticipated extreme scenarios that the robust optimization generated. Implementation of these policies is easy, as they only need to change the ambulance stationings to the optimal ones we calculated. Table 3 .: The results of MRP batches on the simulation. 

Simulations can also help a city understand what to do when they have to permanently add or remove an ambulance from their fleet. An ambulance may be added if a city's population is growing and the department has new funds to employ a new ambulance engendering the question, where should they put it? The opposite scenario is when a department budget is cut and they must remove an ambulance, which one can they possibly remove to have the least impact? Thus optimization can help a city with a small budget maintain performance while removing an ambulance from operation.

We conduct the optimal stationing as the number of ambulances changes in figure 7 . As can be seen, both the adding and the removal of ambulances only cause minor changes. This is good in the sense that if an ambulance must be removed performance does not decrease that much. But disappointingly adding ambulances also does not increase performance. These differences we found may be magnified though, as the grid time adjustment we saw earlier tends to "squeeze" the times together. 

To test whether socio-economic factors play a role in influencing the EMS response time across different neighborhoods across Austin, we applied the following protocol:

Dataset For each census tract, we calculated the mean reported travel time of all incidents located within the census tract as the dependent variable. As for the independent variables, we included the minimum grid time and the average grid time within each census tract. This was calculated in two steps. First we calculated the minimum and average grid times for each grid cell. Then we took the weighted average (based on the call frequency) of grid cells within each census tract to represent them. Furthermore, we included 19 socio-economic independent variables extracted from the 2018 CDC Social Vulnerability Index dataset [1] . The variables included in this study are listed in Table B1 in the appendix.

We applied 5-fold cross-validation. For each fold, we included 80% as the training set and the remaining 20% as the test set. Before training predictive models, we normalized each variable of the training set by subtracting the mean and then dividing by the standard deviation. We then applied the same to the test set using the sample mean and standard deviation from the training set.

Models We applied 5 models in total. The baseline model was simply the average of all the census tracts' mean reported travel time. Because we believe the reported travel time is correlated with the grid time, we then fitted 3 linear regressions models: the first model used minimum grid time as the sole independent variable; the second model used average grid time as the sole independent variable; the third model included minimum grid time and average grid time as its independent variables. The last model is a LASSO (Least Absolute Shrinkage and Selection Operator) regression model incorporating all 19 socio-economic factors besides minimum grid time and average grid time.

Average MSE on the test set To compare these 5 models, we calculated the mean squared error (MSE) on the test set for each fold. Then we took the average among these 5 folds. If the average MSE of the LASSO model incorporating socio-economic factors is lower than all the other models, then we have evidence to believe that socio-economic factors contribute to differences in reported travel time among various census tracts. Otherwise, we believe that the reported travel time solely depends on geo-location factors reflected by the grid time and is not explicitly affected by socio-economic factors.

As reported in table 4, we found that the lowest average MSE was achieved by the linear regression model fitted on the average grid time. Thus, we conclude that socio-economic factors may not explicitly affect the travel time for ambulances in Austin.

Variable ( Table 4 .: Comparing the 5 regression models.

OpenEMS is open-source and published at https://github.com/joshua-ong/ AmbulanceDeployment. A user guide of the software package can be seen in figure 8 ; an overview of the software variables are in table 5. 

Description | regions |

The number or regions is how many uniform discrete regions the continuous space is made into. This can be though of as the resolution of the grid. grid

An object containing the information of the discrete grid. This includes a time matrix between any two grid cells. Integer vector in Z |Stations| where index i is the optimal number of ambulances in station i. ncalls number of calls to simulate nbatches number of batches to simulate Table 5 .: List of variables for user guide.

Open Street Map module creates the discrete grid described in section 4.2. Inputs include the number of regions, as well as the minimum/maximum longitudes/latitudes that bound the city of interest. The module returns the grid object which can be saved as a .json file. This grid object includes both the data about the grid cells as well as useful functions to handle that data. For socioeconomic interest, it can also contain a mapping of census tracts to grid points.

> grid = create grid(num regions,min lattitude,max lattitude ,min longitude,max longitude) Data Pre-processing module collects the city EMS call data with headers for date of call, latitude and longitude of the emergency call for training the linear program and simulations. Further headers for reported travel time, reported service time, and latitude and longitude of the ambulance at the time of the call are needed for calibration and verification. The data is then processed into an adjacent neighborhood matrix, coverage matrix, LP train call matrix, calibration train call matrix, and test call matrix and saved. The user may consider customizing training and testing sets. For example, if the city's data changed drastically because of the Covid-19 pandemic, they could train before Covid-19 and test during Covid-19.

> (adjacent nbhd,coverage,LP train call,callibration train call,test call) =data pre processing(city data = "path/to/city data.csv",grid)

Calibration module includes fitting a regression between grid times and reported travel times to calibrate the model. There is an example in section 5.3. Datasets should be trained with the callibration train set. We included linear and log-log linear regressions, but we suggest the reader look at other regressions to find what works best for their dataset. The result is the adjusted travel time.

Verification module measures the accuracy of the adjusted travel time. There is an example in 5.4. The test calls vector includes the location and time of the call, as well as the reported response time and longitude and latitude of the ambulance at the time of the call. This is so we can verify the ambulance reported travel path time with our grid. The calibration object takes grid time and outputs adjusted grid time. The outputs are the differences in adjusted grid time and reported time. This can be measured as a score or as a plot.

> verification(test call,callibration = log callibration)

Linear Program Solvers module solves the two stage stochastic and robust programs detailed in section 3. Besides the pre-processed data, the other input is the α value that determines how conservative the robust model is. The function then returns the solved optimal stationing x stoch * and x robust * . > x stoch * , x robust * = generate models(α = .1,LPtrain calls,adjacent nbhd, coverage,callibration=log callibration,grid)

Simulation module performs numerical computations using the test calls, optimal stationing, and adjusted grid time. There is an example in 5.8. The inputs were all computed in previous sections. The output is the computed response time. From this one can compute their metrics.

> simulated response time = simulation(x,test calls,grid,n calls = 1000,lognormal params = (µ,σ),callibration = log callibration)

The results will require unique analysis for each case study. All of the visualizations from this paper are included and can be reproduced.

This paper presents an open-source project OpenEMS, an end-to-end pipeline to optimize ambulance strategic decisions. It includes two-stage stochastic and robust programming to optimize ambulance stationing and routing. In this open-source package, we include many of the details that prevent state of the art research from being applied on a case study level. This includes constructing a realistic grid using open street maps, calibrating these travel times to account for ambulances' special privileges, and verifying they work how we think they should work. Further, we provide practical guidance to hypothetical questions like the case where an EMS department has to add or remove one or more ambulances from their fleet. This is a common question.

We then applied this software to a case study on the city of Austin, Texas. In it, we found we could increase the mean response time by 88.02 seconds if the stationing was changed. We also showed the stochastic optimization performed better than the robust optimization. We also found the ideal locations to add or remove an ambulance to increase or maintain performance respectively.

There is a lot of avenues for future work. Now that this code is open-source, one could simply try some new mathematical program and simulate for new results, for example, the distributionally robust optimization (DRO) that is growing in popularity in the operation research community. Along the lines of the simulation, we want to include more features such as sensitivity analysis with the choice of the grid resolution. There are also more questions about robust optimization; one main problem is the non-convex nature of the VaR constraint in the formulation. This leads to some unexpected results. For example, alpha = .1 and alpha = .01 could work fine, but this does not give any guarantees on interpolations between those values.

Funding details Figure A1 .: Ambulance location schemes for the stochastic optimization (left column) and robust (α = .01) optimization (right column). Rows start with 30 ambulances, increase by 5 ambulances per row, and end with 50 ambulances. The purple circle indicates there is an ambulance station there. The radius of the circle is proportional to the number of ambulances in the area where the biggest circle is 5 and the smallest circle is 0. Figure A2 .: A scatter plot of all call data. Note they cannot represent density. Figure A3 .: A visualization of the discrete grid over Austin county.

Centers for disease control and prevention/ agency for toxic substances and disease registry/ geospatial research, analysis, and services program

A review on simulation models applied to emergency medical service operations

Reducing ambulance response time using simulation: The case of val-de-marne department emergency medical service

Two-Stage Stochastic Integer Programming: A Brief Introduction

The maximal expected covering location problem: Revisited

Designing robust emergency medical service via stochastic programming

Robust and stochastic formulations for ambulance deployment and dispatch

Julia: A fresh approach to numerical computing

Introduction to Stochastic Programming

Do warning lights and sirens reduce ambulance response times? Prehospital Emergency Care

A location-based comparison of health care services in four u.s. states with efficiency and equity

The affordable care act and ambulance response times

Health care logistics, emergency preparedness, and disaster relief: New challenges for routing problems with a focus on the austrian situation

Determining emergency medical service vehicle deployment in austin

The turnaround of the London Ambulance Service Computer-Aided Despatch system (LASCAD)

Solving an ambulance location model by tabu search

Openstreetmap: User-generated street maps

Introduction to arena [simulation software

Simulation modeling and optimization using promodel

Resource overkill: we can do more for less

Ambulance Service Planning: Simulation and Data Visualisation

Simulation of single start station for edmonton ems

Dynamic ambulance dispatching: is the closest-idle policy always optimal?

Effect of ambulance station locations and number of ambulances to the quality of the emergency service

Mohanavalli Rajagopal Lakshmanan, Yih Yng Ng, and Ong. Marcus Eng Hock. Ambulance deployment under demand uncertainty

Ambulance deployment under demand uncertainty

Distributionally robust optimization of an emergency medical service station location and sizing problem with joint chance constraints

Computing in Operations Research Using Julia

Emergency response time and pre-hospital trauma survival rate of the national ambulance service

Simulation and real-time optimised relocation for improving ambulance operations

Approximate dynamic programming for ambulance redeployment

A simulation model to enable the optimization of ambulance fleet allocation and base station location for increased patient survival

A model for optimally dispatching ambulances to emergency calls with classification errors in patient priorities

A case study of optimal ambulance location problems

Locating an ambulance base by using social media: a case study in bangkok

Eight minutes or less: does the ambulance response time guideline impact trauma patient outcome?1 1Selected Topics: Prehospital Care is coordinated by Peter T. Pons, MD, of Denver Health Medical Center

The maximum availability location problem

Solving the dynamic ambulance relocation and dispatching problem using approximate dynamic programming

Dynamic ambulance deployment to reduce ambulance response times using geographic information systems: A case study of odunpazari district of eskisehir province, turkey

A comprehensive survey on the ambulance routing and location problems

Health care spending growth: How different is the united states from the rest of the OECD?

Solving two-stage robust optimization problems using a column-and-constraint generation method

Three-fold crossvalidation on 5 alpha values to find when α saturates

We also want to thank the Senior Project team from University of Texas at Austin who helped code sections of this open-source project: Ethan Santoni-Colvin, Michael Hilborn, Sudeep Narala, Xander Tedjo, Will Worthington

The authors report there are no competing interests to declare.

Appendix B. List of variables extracted from CDC SVI dataset