The EMS vehicle patient transportation problem during a demand surge

Abstract

We consider a real-time emergency medical service (EMS) vehicle patient transportation problem in which vehicles are assigned to patients so they can be transported to hospitals during an emergency. The objective is to minimize the total travel time of all vehicles while satisfying two types of time window constraints. The first requires each EMS vehicle to arrive at a patient’s location within a specified time window. The second requires the vehicle to arrive at the designated hospital within another time window. We allow an EMS vehicle to serve up to two patients instead of just one. The problem is shown to be NP-complete. We, therefore, develop a simulated annealing (SA) heuristic for efficient solution in real-time. A column generation algorithm is developed for determining a tight lower bound. Numerical results show that the proposed SA heuristic provides high-quality solutions in much less CPU time, when compared to the general-purpose solver. Therefore, it is suitable for implementation in a real-time decision support system, which is available via a web portal (www.rtdss.org).

1 Introduction

This paper considers dispatching EMS vehicles in emergencies, e.g., a large scale pandemic viral attack such as COVID-19, when the demand for vehicles is high. It is developed in a way that makes it suitable for use during demand surge as well as during normal operations. As suggested in [1,2,3] and as is current practice, EMS vehicles may need to serve up to two patients in these situations.

The EMS vehicle patient transportation problem (EMSVPTP) we consider assigns vehicles to patients so they can be transported to hospitals they choose or those selected by EMS personnel. During normal operations, patients may choose hospitals based on their medical insurance considerations. During a medical surge, hospitals may be chosen by EMS personnel based on availability, wait time, drive time and patient condition. Each EMS vehicle can serve up to two patients. The objective is to minimize the total travel time among all vehicles subject to two time-windows—the time window within which a vehicle must arrive at a patient’s location and another time window within which the patient must be transported to the hospital. Two reasons have led to the use of total travel time as the objective function. First, an industry advisory board recommended that models for emergency response should not only be designed for use during a demand surge, but also during normal operations. This will ensure dispatchers trained to use a tool during normal operations can use it effectively when a demand surge occurs due to a pandemic. Second, minimizing total travel distance is a commonly used objective in software adopted by EMS departments in the U.S. and Canada, as they seek to minimize their overall cost while incorporating response time compliance as constraints. The proposed model has been implemented in a real-time decision support system, which is available on a web portal (www.rtdss.org).

The remainder of this paper is organized as follows. Section 2 reviews the relevant literature. Section 3 introduces the EMSVPTP and describes its computational complexity. Section 4 presents a mixed-integer programming (MIP) model of the EMSVPTP. A column generation method to obtain a lower bound for the EMSVPTP is presented in Sect. 5. Next, a simulated annealing heuristic method for the EMSVPTP is presented in Sect. 6. Numerical experiments and computational results with CPLEX [4] are provided in Sect. 7. Conclusions are drawn in Sect. 8.

2 Literature review

The problem of dispatching EMS vehicles has been considered in the literature from two main aspects: the dispatching decision to serve current patient requests; and coverage decision in anticipation of serving future patient requests. An overview of these papers can be found in [5,6,7].

Yang et al. [8] develop an integrated approach for dispatching EMS vehicles, police vehicles, and fire trucks. Their model can be used in real-time to decide a subsequent destination for each vehicle. The objective of their model is to minimize travel times. Penalties are incurred for not reaching a patient in the desired time window and not covering demand nodes. Other researchers have also considered the coverage and dispatching decisions jointly. For example, Andersson and Varbrand [9] study an EMS vehicle dispatching and dynamic relocation problem. Ibri et al. [10] study an integrated dispatch and coverage problem model and develop a method combining ant colony and tabu search techniques for its solution. Majzoubi et al. [11] propose a mathematical model to dispatch and relocate these vehicles when the demand for the EMS vehicles is high. In the current paper, we only deal with dispatching part of the problem and the focus is on utilizing the capacity of EMS vehicles during a demand surge.

Mathematically, the EMS dispatching problem can be considered as a real-time vehicle routing problem (VRP) which is addressed abundantly in the literature. Because each patient has pickup and drop-off points, it has characteristics that are similar to the VRP with pickup and delivery (VRPPD). Gendreau et al. [12] propose a tabu search heuristic for the VRPPD and use parallel computing for real-time fleet dispatch whenever there is a new request. Bent and Hentenryck [13] develop a hybrid algorithm for the VRPPD. They use simulated annealing for minimizing the number of vehicles and the travel cost. Further, Ropke and Pisinger [14] propose an extended large neighborhood search algorithm for the VRPPD.

Another specific type of VRP that is related to the EMSVPTP is the dial-a-ride (DAR) problem which has also been studied in the past. A survey of models and algorithms for the DAR problem is provided in Cordeau and Laporte [15]. A more recent survey of the DAR literature is provided in Ho et al. [16]. Using the DAR framework, Beaudry et al. [17] consider a problem for transporting patients on a hospital campus. Their problem considers the dynamic arrivals of new transportation requests. The objective in [17] is to minimize a weighted sum of three criteria: total travel time, lateness, and earliness. The authors use insertion movements and tabu search to develop a two-phase heuristic procedure in solving their particular DAR problem. Kergosien et al. [18] consider the problem of transporting patients between care units. They use a tabu search method similar to the method used in [12]. Parragh [19] develops two mathematical models for the DAR problem similar to those in [15] and proposes a branch-and-cut algorithm with a variable neighborhood search to solve the problem. Hu and Chang [20] propose a revised branch-and-price for the DAR problem, while Liu et al. [21] present a branch-and-cut algorithm for the problem. Markovic et al. [22] develop a Mobile Resource Management System for the DAR problem, designing a customized heuristic for dynamic dispatching. They consider a case study in Maryland and show their proposed solution can reduce the expenses for that case by 18%.

3 Computational complexity

Krumke et al. [23] study a general vehicle dispatching problem with at most two requests (VDP2), which simply dispatches vehicles so that the total transportation cost is minimized. Note that unlike a k-customer VRP (e.g., Haimovich et al. [24]), all vehicles in VDP2 are geographically dispersed and thus are not necessarily based at one depot. Furthermore, Krumke et al. [23] show that the VDP2 is NP-complete. In analyzing the computational complexity for the EMSVPTP, we demonstrate that essentially VDP2 is reducible to the EMSVPTP in polynomial time.

First, we formally introduce the VDP2 in [23]. We then show that a slight variation of VDP2, namely, a vehicle dispatching problem with two requests and a common destination depot (VDP2CDD), is NP-complete. Finally, we establish the proof that the EMSVPTP is NP-complete through a reduction from VDP2CDD.

Vehicle dispatching problem with two requests (VDP2): Given a set of requests R, a set of vehicles U ($|U|\le 2|R|$), cost function $c:\ (U\cup R)\times (U\cup R)\mapsto R^+$, the vehicle dispatching problem with two requests (VPD2) finds a tour $t_u=(u, r_{u,1}, r_{u,l(u)})$ for each vehicle $u\in U$ that serves $l(u)\le 2$ requests, such that each request is served in exactly one tour and the total cost of all tours is minimized.

Vehicle dispatching problem with two requests and a common destination depot (VDP2CDD): Given a set of requests R, a set of vehicles U ($|U|\le 2|R|$), a common destination depot d, cost function $c:\ (U\cup R)\times (U\cup R\cup \{d\})\mapsto R^+$, the vehicle dispatching problem with two requests and a common destination depot (VDP2CDD) finds a tour $t_u=(u, r_{u,1},r_{u,l(u)},d)$ for each vehicle $u\in U$ that serves $l(u)\le 2$ requests, such that each request is served in exactly one tour, all tours end at d, and the total cost of all tours is minimized.

Lemma 1

The VDP2CDD is NP-hard.

Proof

Consider any VDP2 instance (R, U, c). Construct a VDP2CDD instance (R, U, c, d), where $c(r,d)=a\ne 0$ for all $r\in R$ and a is a constant. As described, the construction of VDP2CDD can be done is polynomial time. In addition, any optimal solution to VDP2 solves the above constructed VDP2CDD and vice versa. In Krumke et al. [23], it is shown that the decision version of the VDP2 is NP-complete, thus VDP2 is NP-hard. Therefore, VDP2CDD is also NP-hard.

To prove EMSVPTP’s complexity, we restate the EMSVPTP as follows. $\square$

EMS Vehicle Patient Transportation Problem (EMSVPTP): Consider a set of patients P, a set of their corresponding requested hospitals $H=\{h(p)|\forall p\in P\}$ ($|H|=|P|$), and a set of vehicles V ($|V|\le 2|P|$) each with a capacity of k patients. Define $N=P\cup H$ and the cost function $s:\ (V\cup N)\times (V\cup N)\mapsto R^+$. Then, the EMS vehicle patient transportation problem, denoted as EMSVPTP(V, P, H, s), finds a tour $t_v=(v, n_{v,i_1}, n_{v,i_2}, \ldots , n_{v,i_{2l(v)}})$ for each vehicle $v\in V$ that serves at most $l(v)\le k$ patient requests, i.e., picks up each patient and drops each patient off at the requested hospital, such that each patient is served in exactly one tour and the total cost of all tours is minimized.

Note that in the above definition for EMSVPTP, $(i_1,i_2,\ldots ,i_{2l(v)})$ is a permutation for $\{1,2,\ldots , 2l(v)\}$. Further, if patient $p_0$ is the $i_t$th stop for vehicle v in its tour and its associated hospital $h_0=h(p_0)$ is the $i_q$th stop for the vehicle, i.e., $p_0=n_{v,i_t}$, $h_0=n_{{v,i_q}}$, then $i_q>i_t.$

Theorem 1 establishes the NP-hardness result for the EMSVPTP with $k=2$ from the reduction of the VDP2CDD problem in Lemma 1 in polynomial time.

Theorem 1

The EMS patient transportation problem with $k=2$ is NP-hard.

Proof

Consider any instance of VDP2CDD (R, U, c, d). To construct an equivalent EMSVPTP with $k=2$, let $V=U$, $P=R$, $s=c$, and $H=\{d,d,\ldots ,d\}$, i.e., all the destination hospitals are at the same location. Clearly, the EMSVPTP(V, P, H, s) problem can be constructed in polynomial time.

Let $t^*=\{t^*_u\}_{u\in U}$ be an optimal solution to VDP2CDD. Then, for each vehicle u, the optimal tour serves either one or two requests.

• (i) If $u\in U$ serves one request , i.e., $t^*_u=(u,r^*_{u,1},d)$, then the corresponding solution to the EMSVPTP with $k=2$ is $w^*_v=(v,p^*_{v,1}=r^*_{u,1},h=d).$

• (ii) If $u\in U$ serves two requests , i.e. $t^*_u=(u,r^*_{u,1},r^*_{u,2},d)$, the corresponding solution to the EMSVPTP with $k=2$ is $w^*_v=(v,p^*_{v,1},p^*_{v,2},h=d)$.

Note that (ii) is true because both requests have a common destination. We show that $\{w^*_v\}_{v\in V}$ is optimal to the above-constructed EMSVPTP with $k=2$ by contradiction. If there exists a solution $\bar{w}=\{ \bar{w}_{v\in V}\}$ with total cost $\sum _v \bar{s}_v \le \sum _v {s^*_v}$, then solution $\bar{w}$ will induce an alternative solution ${\bar{t}}$ using (i)–(ii) with lower cost for VDP2CDD. This contradicts that $t^*$ is an optimal solution to VDP2CDD. Thus, any optimal solution to VDP2CDD solves the above constructed EMSVPTP. Similarly, it can be shown that any optimal solution $w^*$ to the EMSVPTP with $k=2$ solves the VDP2CDD. Thus, the VDP2CDD is equivalent to the constructed EMSVPTP. From Lemma 1, EMSVPTP with $k=2$ is NP-hard.

Due to the above complexity results for EMSVPTP, we focus on developing efficient heuristic methods for solving EMSVPTP in a real-time decision support system. Because the current paper deals with the EMSVPTP with $k=2$ as a prototype problem, in subsequent sections we refer to the EMSVPTP with $k=2$ simply as EMSVPTP without explicitly specifying $k=2$. $\square$

4 The mathematical formulation for the EMSVPTP

To formulate a mixed integer linear program for the EMSVPTP, we define the following sets and mapping functions:

Sets
K	Set of vehicles
$K^H\subseteq K$	Set of vehicles en route to hospitals
P	Set of patients
H	Set of hospitals
O	Set of virtual depots for vehicles
Mapping functions
$o(k)\in O$	A function mapping vehicle k to its virtual depot
$p(k)\in P$	A function mapping vehicle k to the patient
	currently on board, $\forall k\in K^H$

In the above notation, O refers to a set of virtual depots for the vehicles. It is assumed that the current location of an EMS vehicle is that of its (virtual) depot. Based on the defined sets and mapping functions, and following the notation in the three-index formulation in [15], we define a network $G=\{N,A\}$, where $N=O\cup P\cup H\cup \{d\}$ is the set of all nodes including virtual depots, patients, hospitals, and a common virtual destination node d for all vehicles. After the last patient is dropped off at a hospital, it is assumed that the EMS vehicle travels to the virtual destination node d, for which the travel time is assumed to be zero.

Assuming there are m EMS vehicles and n patients, we index all nodes in N as follows:

$$\begin{aligned} \begin{array}{l} P=\{1,2,\ldots ,n\},\ H=\{n+1,n+2,\ldots ,2n\},\\ O=\{2n+1,2n+2,\ldots ,2n+m\},\ d=\{2n+m+1\}. \end{array} \end{aligned}$$

Note that there is a distinct hospital node corresponding to each patient node. In reality, two or more of the hospital nodes may correspond to the same hospital, but we formulate our model so that the transfer of each patient to a hospital has a unique arc associated with it, following the convention used in general VRPs with unique pick-up and drop-off points.

In the EMSVPTP, there are six possible vehicle movements at any given time as listed below.

• 1. From the virtual starting depot to a patient;

• 2. From first patient to a second patient;

• 3. From a patient to a hospital;

• 4. From a hospital to pick up a second patient;

• 5. From first patient’s hospital to the second patient’s hospital and vice versa;

• 6. From hospital to the common virtual destination.

Therefore, the set of valid arcs A is defined as

$$\begin{aligned} \begin{array}{rl} A=\{(i,j)| &{} i\in O, j\in P\cup \{d\}; \\ &{} i\in P\cup H, \ j\in P\cup H, \text{ for } i\ne j \text{ and } i\ne n+j;\\ &{} i\in H, j=2n+m+1\}, \end{array} \end{aligned}$$

where the first set of (i, j)’s for $i\in O,\ j\in P$ represents the first type of EMS vehicle movement, the third set of (i, j)’s represents the sixth type of vehicle movement, and the second set represents the remaining four types of movements.

Let binary variable $x^k_{ij}$ (where $(i,j)\in A$ and $k\in K$) be 1 if vehicle k goes from node i to j; 0 otherwise. Also, let variable $y_{ik}$ represent the time that vehicle k arrives at node i. Let $t_{ij}$ be the travel time on arc $(i,j)\in A$ and $\psi _i$ be the desired time by which node i must be visited. Then, similar to [15], the EMSVPTP can be formulated as follows:

$$\begin{aligned} \ \ \ \min&\ \ \ \sum _{(i,j)\in A}{\sum _{k \in K}{t_{ij}x^k_{ij}}}&\end{aligned}$$

(1)

$$\begin{aligned} \text {s.t.}&\ \ \ \sum _{j \in N}\sum _{k\in K} x^k_{ij} = 1,&\forall i\in P \end{aligned}$$

(2)

$$\begin{aligned}&\ \ \ \sum _{i \in N}x^k_{2n+k,i}=\sum _{i \in N}x^k_{i,2n+m+1}=1,&\forall k\in K \end{aligned}$$

(3)

$$\begin{aligned}&\ \ \ \sum _{j\in N} x^k_{ij}- \sum _{j\in N} x^k_{j,n+i} =0,&\forall i\in P,\ \forall k\in K \end{aligned}$$

(4)

$$\begin{aligned}&\ \ \ \sum _{j\in N}x^k_{ji} - \sum _{j \in N}x^k_{ij}=0,&\forall i\in P\cup H,\ \forall k\in K \end{aligned}$$

(5)

$$\begin{aligned}&\ \ \ x^k_{2n+k,p(k)}=1,&\forall k \in K^H \end{aligned}$$

(6)

$$\begin{aligned}&\ \ \ \sum _{i\in N}\sum _{j\in P} x^k_{ij}\le 2,&\forall k\in K \end{aligned}$$

(7)

$$\begin{aligned}&\ \ \ y_{ik}+t_{i,i+n}\le y_{i+n,k},&\forall i\in P,\ \forall k\in K \end{aligned}$$

(8)

$$\begin{aligned}&\ \ \ y_{jk}\ge y_{ik}+ t_{ij}-M(1-x^k_{ij}),&\forall (i,j)\in A, \forall k\in K \end{aligned}$$

(9)

$$\begin{aligned}&\ \ \ y_{ik}\le \psi _i,&\forall i\in N, \forall k\in K \end{aligned}$$

(10)

Recall that $x^k_{ij}$ is defined only on permissible arc $(i,j)\in A$ for vehicle k. Then, in the above formulation (EMSVPTP), the objective function (1) minimizes the total travel time for all the vehicles. Constraint (2) states that each patient must be served by a vehicle. Constraint (3) ensures each vehicle begins at the virtual depot and returns to the virtual destination node. Constraint (4) ensures the same vehicle will pick up a patient and transport them to the associated hospital. In other words, if vehicle k serves patient i, vehicle k has to transport that patient to hospital node $n+i$ as well. Constraint (5) is the flow balance equation at either a patient or a hospital node, ensuring the vehicle that comes into node i must leave that node. Constraint (6) ensures that a vehicle en-route to a hospital with a patient is assigned to that patient. Of course, this vehicle may then serve one more patient, as specified by constraint (7). Finally, constraints (8)–(10) ensure the desired time window constraints are satisfied at both patient and hospital nodes.

The above model can be used in a real-time environment for an EMS department with a CAD (computer-aided dispatcher) and AVL (automatic vehicle locations). When latitude and longitude of a vehicle are provided, the origin-destination matrix $t_{ij}$ can be readily calculated using real-time routing applications. To further reduce processing time, the time between street segment centroids can be pre-calculated and the vehicles can be snapped back to these segments when we receive the AVL data.

The mixed-integer program (1)–(10) is implemented using C# via CPLEX 12.10 [4]. Because the NP-hardness results from Sect. 3, it is not surprising that general-purpose solvers such as CPLEX take excessive amount of CPU time for even small-sized instances, as reported in detail in Sect. 7. Thus, the next section proposes an efficient heuristic method for obtaining a lower bound and a high-quality solution faster.

5 Column generation

In this section, we present a column generation approach for obtaining a tight lower bound for the EMSVPTP. We then exploit this method to solve the and find a good feasible solution. First, we present the basics of applying the column generation technique to the general VRP similar to many papers in the literature (e.g., Desrochers et al. [25], Xu et al. [26], Fillet [27], and Ropke and Cordeau [28]). Then we discuss a special pricing algorithm for solving the subproblem in the column generation method for the EMSVPTP.

In the literature, many researchers have used the column generation technique to solve the VRP through a reformulation as a set covering problem (e.g., Xu et al. [26]). When considered as a set covering problem, the VRP is to determine the set of routes with the minimum total travel cost given the entire set of all possible routes.

More formally, let $a_{iku} = 1$, if route u includes a visit of node i by vehicle k; 0 otherwise. Let $c_{ku}$ be the cost of route u for vehicle k. In addition, let $\varOmega$ be the set of all feasible routes and $\varOmega _1\subseteq \varOmega$ be an arbitrary subset of $\varOmega$. A route is feasible if it satisfies the required time window constraints. Suppose decision variable $\theta _{ku}=1$ if route u is served by vehicle k; 0 otherwise. Then, the VRP can be written as the set covering problem as follows:

$$\begin{aligned} \text {(VRP-SC)}&\min&\ \ \ \sum _{k \in V} \sum _{u \in \varOmega }{ c_{ku}\theta _{ku}}\ \end{aligned}$$

(11)

$$\begin{aligned}&\text {s.t.}&\sum _{k \in V} \sum _{u{\in }\varOmega }a_{iku}\theta _{ku} \ge 1, \ \ \ \forall i \in N \end{aligned}$$

(12)

$$\begin{aligned}&\sum _{u{\in }\varOmega } \theta _{ku} \le 1, \ \ \ \forall k \in V\end{aligned}$$

(13)

$$\begin{aligned}&\theta _{ku}\in \{0,1\} \end{aligned}$$

(14)

The column generation approach begins with a subset of routes $\varOmega _1\subset \varOmega$ to solve the following (restricted) master problem of (VRP-SC).

$$\begin{aligned} \text {(MP)}&\min&\ \ \ \sum _{k \in V} \sum _{u \in \varOmega _1}{ c_{ku}\theta _{ku}}\&\end{aligned}$$

(15)

$$\begin{aligned}&\text {s.t.}&\sum _{k \in V} \sum _{u{\in }\varOmega _1}a_{iku}\theta _{ku} \ge 1, \ \ \ \forall i \in N&\end{aligned}$$

(16)

$$\begin{aligned}&\sum _{u{\in }\varOmega _1} \theta _{ku} \le 1, \ \ \ \forall k \in V&\end{aligned}$$

(17)

$$\begin{aligned}&\theta _{ku}\in \{0,1\}, \ \ \ \forall k \in V&\end{aligned}$$

(18)

Consider the linear relaxation of the master problem (MP). Let $\lambda _i$ and $\omega _k$ be the dual variables associated with constraints (16) and (17), respectively. We solve the subproblem of minimizing the reduced cost for all routes $u\in \varOmega$. If the minimum reduced cost is zero, then the optimal route is found; otherwise, we update $\varOmega _1$ by including the route with the minimum reduced cost and resolve the master problem (MP). Therefore, the subproblem is commonly referred to as the “pricing” problem. Rather compactly, the subproblem can be written as:

$$\begin{aligned} \text {(PP)}: \min _{u\in \varOmega } \{-\sum _{i\in N} { a_{iku}\lambda _i} + \omega _k +c_{ku}\} \end{aligned}$$

Let $b_{ijku}$ = 1 if route u uses vehicle k to visit arc (i, j) and 0 otherwise. Then, $a_{iku} = \sum _{j\in N} b_{ijku}$ and $c_{ku} = \sum _{(i,j)\in (N\cup O)\times N}b_{ijku} t_{ij}=\sum _{i,j\in N}b_{ijku}t_{ij} + \sum _{j \in N} {b_{o(k)jku}t_{o(k)j}}$. Thus, the objective function in (PP) can be rewritten as

$$\begin{aligned} c_{ku} - \sum _{i\in N} {a_{iku}\lambda _i}+\omega _k = \sum _{i,j\in N} b_{ijku} (t_{ij}-\lambda _i) + \sum _{j \in N} {b_{o(k)jku}t_{o(k)j}} + \omega _k \end{aligned}$$

(19)

In other words, the travel cost between nodes is $(t_{i,j}-\lambda _i)$ for $i,j \in N$ and $(t_{o(k),j} + \omega _k)$ for $(o(k),j) \in O\times N$. Therefore, the subproblem (PP) can be considered as an elementary shortest path problem with resource constraints (ESPPRC) (see e.g., Desrochers et al. [25] and Ropke and Cordeau [28]), where the resource constraints correspond to the time window constraints in the underlying route selection.

Finally, we make two notes regarding the implementation of the pricing problem (PP). First, because each vehicle cannot serve more than two patients in our problem, the maximum length of a route (excluding depots) is four. Therefore, we simply enumerate all routes and add the best route to $\varOmega _1$. Second, as detailed in Sect. 7, our implementation adds M routes/columns (with the M most negative reduced costs) at a time, where M is a pre-defined parameter representing the maximum number of columns to add in the column generation algorithm. This process is continued until we are not able to find any negative cost routes.

6 A simulated annealing algorithm for the EMSVPTP

In this section, we present a meta-heuristic method for the EMSVPTP using the simulated annealing approach. We define a pairwise exchange type of neighborhood function/search in Subroutine 1, followed by Subroutine 2 which describes how we optimally sequence the nodes visited by a vehicle after the selections of nodes are determined for all vehicles. Finally, Subroutine 3 conducts a specialized local search attempting to introduce more diverse solutions at later stages (with lower temperatures) of the SA procedure. We refer to the combination of this local search with the pairwise exchange neighborhood search as the ‘hybrid’ method in this and subsequent sections.

In particular, Subroutine 1 performs a pairwise type of exchange to obtain a neighborhood solution from any given solution. It first randomly chooses two vehicles in the current solution and randomly selects one patient from each vehicle. The patients are then swapped between the two vehicles. Note that a randomly selected vehicle may not have any patient, thus there are three possible outcomes for the exchange/swap: 1) swapping two patients between the two chosen vehicles, each has at least one patient; 2) moving one patient from a vehicle with at least one patient to the other that has no patient; 3) doing nothing when neither vehicle has a patient. Subroutine 1 below excludes all swaps with outcome 3).

Subroutine 2 OptimizeSequence determines the best sequence of all nodes to be visited by vehicle k when it serves two patients. In this case, at most six possible sequences exist because there are, at most, two patients and two hospitals to be visited by a vehicle, and each patient is visited before the associated hospital. In particular, let $1\le \rho \le 6$ be the sequence index, $p_1$ and $p_2$ be the two patients and $h_1$ and $h_2$ be their respective hospitals. Then sequences 1 through 6 are shown in the table below. For example, in sequence 3, vehicle k visits patient $p_1$ first, followed by that patient’s hospital $h_1$, then patient $p_2$, and finally the hospital $h_2$ of patient 2.

Sequence index $\rho$	The actual sequence
1	$(k,p_1,p_2,h_1,h_2)$
2	$(k,p_1,p_2,h_2,h_1)$
3	$(k,p_1,h_1,p_2,h_2)$
4	$(k,p_2,p_1,h_1,h_2)$
5	$(k,p_2,p_1,h_2,h_1)$
6	$(k,p_2,h_2,p_1,h_1)$

Using this notation, Subroutine 2 examines all applicable sequences l through u ($l\le u$) for vehicle k in serving patient j, with X being the associated solution. If improvement is found, the subroutine updates X with the improved sequence. Note that for vehicle $k\in K\setminus K^H$, all six sequences are feasible, hence $l=1$ and $u=6$. However, for vehicle $k\in K^H$, because only sequences 1 through 3 are feasible, $l=1$ and $u=3$.

In addition to the pairwise exchange neighborhood search in Subroutine 1, we propose a nested vehicle search (NVS) algorithm to be used in later stages of the SA process when the temperature is relatively low. In essence, the NVS in Subroutine 3 attempts to diversify the solution space when the SA process begins to concentrate its search. In particular, beginning with the current incumbent solution in the SA process, the NVS first randomly select a patient (j) and cancels the assignment of patient $p_1$ to its current hospital $h_1$. Second, a local greedy search of vehicle is performed for this patient j. This search checks all vehicles with three possible scenarios. In the first scenario where vehicle i is assigned to only one patient, step 1.2.1 checks if patient j can be reached by vehicle i within the specified time window. In addition, Subroutine 2 will be applied to optimize the sequence for this vehicle to visit two patients considering six alternative sequences. In the second scenario where vehicle i is transporting a patient to a hospital, step 1.2.2 checks all three possible sequences when performing sequence optimization. This is because the current patient has to be the first patient, thus making $\rho =4$, 5 and 6 infeasible in Subroutine 2. In the third scenario where a vehicle i is assigned to serve two patients, let $j_1$ be the immediate patient and $j_2$ the second patient. Then step 1.2.3 checks if swapping current patient j with the second patient $j_2$ can lead to a better solution. In addition, step 1.2.3 examines if $j_2$ can be served by either a vehicle which is assigned to a patient (step 1.2.1) or a vehicle en route to a hospital (step 1.2.2). If so, the swap will be accepted. After considering all three scenarios for patient j, the best assignment is obtained. This process repeats for all patients other than j before a best local solution is identified.

7 Numerical experiments

To evaluate the proposed simulated annealing heuristic, we used C# to implement the SA heuristic, as well as the mixed-integer model for EMSVPTP via the ILOG/Concert technology by CPLEX 12.10 [4]. The CPU times reported here are from an Acer Intel(R) Core(TM) i7- 6500 U CPU dual-core processor with 2.5 GHz and 2.6 GHz and 8 GB RAM on a 64-bit Windows operating system.

Preliminary pilot studies are conducted to fine-tune the parameters in the proposed SA heuristic. From the pilot runs (see Majzoubi [11]), the best parameter setting recommendations are as follows. The initial temperature $T_{init}=20$ which equates to approximately 80% probability of accepting an uphill transition initially. The choice of the exponential cooling function is $T_{k}=T_0\alpha ^k$ with $\alpha =0.9$. The total epoch length and total length are set to $L=10{,}000$ and $L_{total}=100{,}000$, respectively. For the column generation method, we set the maximum number of columns to add in one iteration to be $M=800$.

In our numerical experiments, all test instances are randomly generated in a $50 \times 50$ mile grid. The location of patients and vehicles are chosen randomly using a uniform distribution within the grid.

To evaluate the efficacy of allowing up to two patients in a vehicle, we run 20 instances with the numbers of patients and vehicles varying from 15 to 25. We compare two dispatching strategies: (a) a traditional dispatching where the closest vehicle is dispatched to serve a single patient, and (b) the proposed method allowing an EMS vehicle to serve up to two patients. For the first strategy, the average time for a patient to arrive at their hospital is approximately 28.5 min whereas the same service time reduces to 19.6 min under the second strategy. This demonstrates the effectiveness of the proposed model. A more thorough analysis regarding the benefit of increasing EMS vehicle’s capacity to two is done in [11].

To measure the quality of our proposed algorithm, we first compare the performance of the SA heuristic for EMSVPTP against that of the MIP model in Sect. 4 by CPLEX [4] on 20 test instances. These test instances are relatively small-sized, with the numbers of patients and vehicles ranging between 15 and 25. Multiple instances are created for each combination of patient and vehicle sizes. The results are compared and summarized in Table 1 which displays the numbers of patients and EMS vehicles for each group of instances in column 2; the objective function values and the computer processing unit (CPU) times for SA in columns 3 and 4; and the objective function value as well as CPU times for CPLEX in columns 5 and 6. Column 7 represents the percentage gap between SA and CPLEX. Also reported in Table 1 are the average objective function values and CPU times for SA and CPLEX.

From Table 1, it can be seen that for the small-sized EMSVPTP instances, the average CPU time for SA is 0.9 seconds compared to 118.38 seconds for CPLEX while the quality of the SA solution is rather high with an average gap of 0.09%. More specifically, in four of the 20 instances, the CPLEX solution is better than those provided by SA. On the other hand, for 16 examples, CPLEX and SA found the same optimal solution. In all of these instances, less CPU time is required for the SA than CPLEX. In particular, SA was 130 times faster than CPLEX, on average. Thus, we conclude that the proposed SA algorithm is efficient and robust, because for most instances it was able to provide the same solution as CPLEX in much less CPU time.

Next, we generated medium-sized EMSVPTP instances where the numbers of vehicles and patients range between 70 and 80, and between 55 and 95, respectively. These 50 instances are too large to be solved by CPLEX, therefore we compare the heuristic solution obtained using SA to the lower bound obtained using the column generation method described in Sect. 5.

In addition to the proposed algorithms, we compare our results with a greedy heuristic algorithm for bench-marking purposes. Subroutine 4 presents an outline of this algorithm.

Table 1 Comparison between the SA and CPLEX solutions

Table 2 Comparison between the SA and the lower bound from column generation

The results of these algorithms are provided in Table 2. All associated master and pricing problems for the set covering problem and all routes for the column generation algorithm were solved by CPLEX via Concert technology and C#. In Table 2, the size of the instances are (70, 55) for instances 1–5, (70, 65) for instances 6–10, (70, 75) for instances 11–15, (70, 85) for instances 16–20, (70, 95) for instances 21–25, (80, 55) for instances 26–30, (80,65) for instances 31–35, (80, 75) for instances 36–40, (80, 85) for instances 41–45 and (80, 95) for instances 46–50 wherein the first argument is the number of vehicles and the second argument is the number of patients i.e., (# of Patients, # of Vehicles). Columns 2 and 3 show the objective function values and CPU times for the column generation algorithm. Columns 4 and 5 show the objective values and CPU times for the proposed hybrid simulated annealing algorithm. The gap between the SA and the lower bound solutions are provided in Column 6. Finally, Columns 7, 8, and 9 show the same information for the greedy algorithm introduced earlier in this section.

Table 3 Comparison between the pure SA and hybrid SA

As shown in Table 2, the average gap of the solutions to the lower bound is 2.45% indicating that the SA algorithm was able to provide a solution very close to the lower bound. In contrast, the greedy algorithm is on average 6.28% worse than the lower bound. Thus, Table 2 shows that the proposed SA algorithm provides good quality solutions for large scale EMSVPTP instances.

In another experiment, we tested our algorithm for larger sized problems, specifically, more than 200 patients and 150 vehicles. Because none of the algorithms besides our proposed SA algorithm was able to find useful solutions, we do not show those results.

Finally, we compare the results obtained using the pure and hybrid SA methods. The pure SA method refers to the SA method using purely pairwise exchange neighborhood function in Subroutine 1, and the hybrid SA method refers to the one using Subroutine 1 and the nested vehicle search neighborhood function in Subroutine 3. As mentioned in Sect. 6, the tests on the hybrid method in Majzoubi et al. [11] suggest that the benefit of including the nested vehicle search is only significant after the temperature in SA is relatively low. Thus the hybrid SA method activates the nested vehicle search when the algorithm is in the last 5% of its iterations.

Table 3 shows the results of this hybrid method on the same 50 instances as in Table 2. Overall, this table suggests that the hybrid method can provide better solutions for 33 out of 50 instances with an average improvement on the objective value of 0.5%. For instance 23, the improvement is greater than 10%, showing the value of the nested local search algorithm.

8 Conclusions

In this paper, we consider the EMSVPTP problem during a demand surge. It minimizes the total time traveled by all vehicles while ensuring all patient requests are served within pre-specified time window constraints. To better utilize the EMS vehicles’ capacity, each EMS vehicle can serve up to two patients. We present a mixed-integer programming formulation for the EMSVPTP problem, similar to the vehicle routing problem with pickup and delivery, based on the three-index formulation in [15]. We show the problem is NP-hard, and thus propose a simulated, we propose a simulated annealing heuristic method for its solution for real-world applications. Our numerical results show that the SA method is efficient in solving large-scale instances for EMSVPTP, when compared to exact solvers such as CPLEX. Future research could involve developing a customized exact algorithm for EMSVPTP using branch-and-price techniques.