AMBULANCE LOCATION UNDER TEMPORAL VARIATION IN DEMAND USING A MIXED CODED MEMETIC ALGORITHM

. Emergency medical services (EMS) are among the most important services in any society due to their role in saving people’s lives and reducing morbidities. The location of ambulance stations and the allocation of ambulances to the stations is an important planning problem for any EMS system to ensure adequate coverage while minimising the response time. This study considers a mixed-integer programming model that determines the ambulance locations by considering the time of day variations in demand. The presented model also considers heterogeneous performance measures based on survival function and coverage for different patient types with varying levels of urgency. A memetic algorithm based-approach that applies a mixed chromosome representation for solutions is proposed to solve the problem. Our computational results indicate that neglecting time-dependent variation of demand can underestimate the number of ambulances required by up to 15% during peak demand. We also demonstrate the effectiveness of the proposed solution approach in providing good quality solutions within a reasonable time.


Introduction
Emergency medical services (EMS) are an integral part of modern healthcare systems and are responsible for providing pre-hospital care and transportation services to both emergency and non-emergency patients.Providing a timely response to emergency calls requiring urgent medical care can save lives and reduce morbidity in patients.The response time, i.e. the time taken to respond to emergency calls, plays a vital role in patient outcomes, necessitating that the ambulances reach the location within the least possible time.The need to respond to calls in a timely manner leads to the ambulance location problem that requires locating ambulance stations and allocating ambulances to these stations such that all demand zones can be reached within a specified time limit.Determining the optimal location of ambulance stations is a strategic-level planning problem, while the optimal allocation of ambulances to these stations is a tactical planning problem.Despite the difference in planning levels, these problems have often been studied simultaneously in the literature [56].The two-level nature of the ambulance location is a critical problem since it influences other tactical and operational level problems such as crew scheduling, relocation, and routing of ambulances.
EMS systems often aim to meet a predefined response time threshold while responding to emergency calls [73].Therefore, response time-based measures are the most commonly used performance measure in ambulance location problems as they are also easy to evaluate [45].Coverage, which measures the expected number of calls to be served within a pre-specified response time, is used in the majority of ambulance location studies [56].However, some studies have also considered other performance measures such as survival function [30] and equity [12].The survival function evaluates EMS system performance based on the probability of patient survival, while equity considers whether demand from all demand zones is equally served.Since an EMS system could serve calls of different urgency levels, location models considering heterogeneous performance measures have been developed [35].
Ambulance location depends on various input factors such as demand arising out of each region, travel time required to reach each location in the region, and response time requirement that needs to be met to achieve the target service performance levels.An important factor that requires consideration while locating ambulances is the probabilistic nature of the problem because ambulances might be busy when a call arrives [22].Many ambulance location models estimate an area-specific busy probability assuming all ambulances in a given area have equal busy probability [38].Some researchers have relaxed this assumption to consider server-level busy probability [7].Another important factor that needs to be considered in determining the location of ambulances is the uncertainty in input parameters [4].Demand from each region exhibits temporal variation, i.e. it usually varies over the day and on various days of the week [11].Considering an average demand over a period of time may result in underestimating the number of ambulances required to achieve the required service level, especially during peak demand periods.Taking into account these variations in demand can improve the overall allocation of ambulances.
Facility location models, including coverage-based location models, are NP-hard [15,16], due to which many heuristic and metaheuristic approaches have been proposed to solve these models [17,34,51].Furthermore, many real-life urban EMS systems serve a large region with a large population.For example, the case considered in our test instances is the city of Mumbai, with a population above 12 million and an area of 619 square km (www.mcgm.gov.in).For such large urban EMS systems, the demand from the region would need to be divided into a large number of zones for improved accuracy and realism, which will result in large-scale optimisation problems [51].Solving such large-scale optimisation problems will require significantly high computational power and time.Although the location of ambulance stations is a strategic problem, ambulance allocation is a tactical problem [56].Therefore, an ambulance location problem is not just required to be solved while initially setting up an EMS system but frequently on a medium-term basis to redistribute ambulances based on the demand forecast.Thus, it becomes necessary to develop effective solution approaches that provide quality solutions within a reasonable time for such realistic problems.
Determining optimal location-allocation of ambulances is thus an important planning problem for the operation of an EMS system due to their direct impact on the survival probability of patients and the overall performance of the EMS system.Incorporating the temporal variation in demand is significant due to the impact on the number of ambulances and ambulance stations required, affecting the busy probability and utilisation of ambulances.As ambulances serve different types of patients, it is also necessary to consider different performance measures to evaluate the performance of the EMS system.Further, ambulance location decisions also affect the operational-level planning of EMS, such as relocation of ambulances and crew scheduling.Based on these observations, we formulate the following research questions to address in this study.
(i) What is the impact of temporal variation in demand on the number of ambulances and ambulance stations located compared to considering average demand?(ii) How does temporal variation in demand affect the number of ambulances and ambulance stations located compared to considering maximum demand throughout the day?(iii) How do the busy probability and server utilisation of ambulances vary as demand varies over the day?(iv) Can a heuristic approach be developed to solve the ambulance location problem with time-dependent variation and server-level busy probability for obtaining better solutions within less time than commercial solvers and the conventional genetic algorithm (GA) for realistic large-scale problem instances?
To address the above research questions, we consider an ambulance location problem that accounts for both temporal variations in demand and server-level busy probability.We also consider heterogeneous performance measures for different patient types, including coverage and survival probability-based measures.A memetic algorithm (MA) based solution approach is developed to solve large-size real-life instances of the problem.The overall contributions of this work can be summarised as follows.
-Present a model for the location of ambulance stations and allocation of ambulances to the stations, which accounts for temporal variation in demand and server-level busy probability.-Develop and implement a novel mixed coded MA-based solution approach to solve large-scale problems, utilising a combination of integer and real coding for solution representation.-Analyse the impact of considering variations in demand over the time of day on ambulance allocation decisions.
The organisation of the remainder of this paper is as follows.Section 2 reviews the existing literature in the area of ambulance location.Subsequently, Section 3 discusses the problem background and the model formulation, followed by Section 4 detailing the proposed solution methodology.In Section 5, a detailed summary of the computational results is presented.Finally, Section 6 presents the conclusions and scope for further research.

Review of literature
A considerable number of papers that study ambulance planning models can be understood from the numerous review papers published related to the topic [4,9].Brotcorne et al. [9] review ambulance location problems considering static, probabilistic and dynamic models.Aringhieri et al. [4] review more recent EMS planning problems across the entire emergency care pathway.Reuter-Oppermann et al. [56] classify EMS planning problems by focusing on the interdependencies between different planning levels and problems.Bélanger et al. [6] review EMS planning problems at the tactical and operational levels, focusing on the interaction of various decisions related to location, relocation, routing and dispatching of ambulances.These review papers provide a comprehensive overview of the overall research in the area of ambulance location.We present an overview of the papers more closely related to the present study.
The location set covering model (LSCM) and maximum covering location problem (MCLP) constitute the earliest major ambulance location models in the literature.Toregas et al. [64] propose LSCM that minimises ambulances required while covering all demand zones.The MCLP, proposed by Church and ReVelle [14], on the other hand, tries to maximise the possible coverage given a fixed number of ambulances.Daskin [22] develops a maximum expected covering location problem (MEXCLP) that addresses the unavailability of ambulances by explicitly considering the busy fraction of ambulances while maximising the expected coverage value.Larson [37] introduces a hypercube queuing model that accounts for the busy probability of ambulances.Davis [23] presents a simulation-based approach that tackles temporal variations in demand and travel time, multiple casualty accidents, and variations in service time at each incident location.The double standard model (DSM) addressing the issue of ambulance unavailability by maximising the demand covered by two ambulances for each location is proposed by Gendreau et al. [31].Various extensions and improvements to these models have been developed and presented in the literature are discussed subsequently.
Hogan and Revelle [33] introduce backup coverage models that maximise the backup coverage for demand zones while enforcing single coverage as a constraint.They analyse the model that considers a trade-off between primary and backup coverages, which they extend to include multiple coverages.Mandell [40] proposes a covering model for a two-tiered ambulance system where advanced life support (ALS) and basic life support (BLS) can serve all calls.A region is considered covered if an ALS is available for service within a pre-specified response time.Revelle and Hogan [57] introduce a probabilistic model called the maximum availability location problem (MALP) to extend the LSCM to account for the busy probability of ambulances.The MALP-I assumes that the busy probability is the same for all ambulances, which is relaxed to develop the MALP-II by considering areaspecific busy probability.Batta et al. [5] develop an adjusted MEXCLP model that assumes that ambulances are independent.Marianov and Revelle [42] present an extension to the MALP that relaxes the assumption that the busy probability of ambulances is independent of each other.They propose a queueing MALP (Q-MALP) model that calculates server-level busy probability using queueing theory to make the model more realistic.
McLay [44] develops the MEXCLP2 model that extends the MEXCLP to account for two server types.
The demand for emergency care follows a definite pattern and varies continuously over the day, with low demand during the night and peak demand during day time [10,11,41].Disregarding this variation in demand while determining the location of ambulances can result in inaccuracies in estimating the number of ambulances required to meet the demand.Repede and Bernardo [55], while developing a decision support system for the Louisville EMS system, present an extension to the MEXCLP to consider time variation in demand.They solve the model for 35 scenarios of the EMS system with the number of ambulances ranging from 5 to 11. Rajagopalan et al. [52] develop a dynamic available coverage location model, a multi-period model applicable for dynamic demand, which accounts for significant changes in demand across different time clusters.Schmid and Doerner [58] develop a mixed-integer linear programming (MILP) model for multi-period planning of ambulance locations while ensuring a pre-specified coverage level throughout the planning horizon.They apply the model for a dataset based on Vienna with potential locations varying from 16 to 163 and 16 available ambulances.Degel et al. [24] take a data-driven optimisation and develop a model that determines the flexible location and fleet size based on empirically determined required coverage, which accounts for the time-varying nature of demand.They demonstrate their approach using the case study of Bochum (Germany), with 163 demand zones, 21 potential ambulance stations and 14 ambulances.Van Den Berg and Aardal [67] also develop a probabilistic model to account for the time-dependent nature of demand and travel time that, while maximising expected coverage, minimises the number of stations and relocations required across different periods.They apply their model to a randomly generated instance with 500 demand zones and 50 potential ambulance stations.
Coverage is the most widely used objective function for evaluating EMS location models.However, it is ineffective in discerning the impact of differences in response times [45].Erkut et al. [30] observe that coveragebased models may locate stations at sites just within the response time limit while increasing coverage to additional regions resulting in decreased survival probability of patients.To overcome this limitation, they introduce a maximum survival location problem to maximise the survival function for patients, which measures the impact of response time on the probability of survival.Erkut et al. [30] solve problems with 180 demand points, 16 potential stations and up to a maximum of 16 ambulances.Knight et al. [35] propose a maximal expected survival location model for heterogeneous patients to account for multiple outcome measures for various patient types and present an approximation methodology that iteratively finds the best solution.They solve the problem based on EMS in Wales with 18 demand nodes, 11 ambulance stations and 36 ambulances.An important objective considered in the ambulance literature is equity in service levels between different zones of a region.Chanta et al. [12] introduce a minimum -envy location problem (MELP) to maximise service equity among all zones.MELP models equity using envy, a function of the distance of a demand zone from its nearest and backup stations.Chanta et al. [13] improve the -envy model introduced in Chanta et al. [12] by considering envy as a function of the difference between survival probabilities to capture more accurately the difference in patient outcomes among different zones.
Leknes et al. [38] formulate a MILP model to consider heterogeneity in demand while accounting for various outcome measures for multiple patient types.They tested the proposed model on instances based on Trondheim and Malvik with 67 demand zones and 44 potential station locations and Sør-Trøndelag with 139 demand zones and 76 potential locations.Yoon and Albert [72] present an EMS system with priority queues for different patient types by employing multiple types of patients to develop a MILP model for ambulance deployment.The developed model is studied using a dataset based on Hanover County, with 270 demand nodes, 16 potential stations and up to 15 ambulances.El Itani et al. [29] extend the MEXCLP to consider the problem of utilising additional private ambulances to improve coverage.Andersson et al. [3] extend the model presented by Leknes et al. [38] to consider various strategic scenarios, including the closure of emergency rooms, usage of designated vehicles for non-urgent patients, and time-dependent variation in demand.Boutilier and Chan [8] consider the problem of locating and routing ambulances in a lower-middle-income country setting of Dhaka, Bangladesh.They use field data and prediction models to handle uncertainty in the data and apply a simulation-based approach to answer policy-related questions.Nelas and Dias [50] present a location model for an EMS system with multiple types of ambulances that serve multiple types of patients, explicitly considering the substitutability between different types of ambulances for different care services and the assignment of ambulances to specific calls.They applied their model to a dataset based on the Coimbra district (Portugal) for locating 35 ambulances among 34 potential ambulance stations.Naji et al. [49] propose a two-server dynamic covering location model to consider two types of patients and two types of servers.Yoon et al. [74] present a stochastic model for joint location and dispatching of ambulances under demand uncertainty.
Due to the NP-hard nature of the ambulance location problem, solution approaches based on various heuristic and metaheuristic approaches have been proposed in the literature.Gendreau et al. [31] develop and apply a solution approach based on tabu-search to a static ambulance location model with double coverage.Doerner et al. [27] propose an ant colony optimisation-based solution approach to extend the double coverage model and compare it with the tabu search metaheuristic.Rajagopalan et al. [51] compare four metaheuristic approaches, including an evolutionary algorithm, tabu search, simulated annealing, and a hybridised hill-climbing algorithm for solving the MEXCLP.Toro-Díaz et al. [65] present a GA-based optimisation framework for the location and dispatching of ambulances simultaneously.Toro-Díaz et al. [66] propose a tabu search-based heuristic for a model that considers fairness in large-scale EMS systems.Zhen et al. [76] develop a simulation-optimisation framework with GA to deploy and relocate ambulances.Akdogan et al. [2] construct a GA to solve the ambulance location model based on an approximate queueing model that minimises the system's response time.Kaveh and Mesgari [34] present an improved biogeography-based optimisation algorithm for solving the MCLP, which they apply to solve the ambulance location problem for a real dataset in Tehran.
Table 1 compares selected articles in the literature based on some key features of the models.The table shows that only Knight et al. [35], Leknes et al. [38], and Andersson et al. [3] consider both server-level busy probability and heterogeneous performance measures (coverage and survival function).However, Knight et al. [35] and Leknes et al. [38] do not consider temporal variation in demand, while Andersson et al. [3] consider stations to be fixed in their time-dependent model.Among the papers that consider relocation of ambulances based on temporal variation, none consider server-level busy probability or survival-based objective.Table 1 also shows the different solution approaches applied to solve the ambulance location problems.Based on our literature review, we identified a need to analyse the impact of temporal variation in demand on the number of ambulances located while incorporating server-level busy probability and heterogeneous performance measures.Therefore, we propose a time-dependent maximum expected performance location problem for heterogeneous patients that accounts for the time-dependent variation in demand and estimates station-specific busy probabilities.The proposed problem also considers an objective function based on the combination of coverage and survival function for different patient types.A mixed-integer non-linear programming (MINLP) model is presented and converted to a MILP model by linearising the non-linear constraints.A memetic algorithm-based approach is proposed to solve the model as the presented model is difficult to solve using commercial optimisation solvers.

Problem description
The problem we address in this work is characterised by a set of demand zones (nodes), where each zone represents a neighbourhood from where calls for ambulances are received.The calls received from demand zones can be of three types (A, B, and C) with varying levels of urgency [35].Type A calls are time-critical emergency calls, where reaching the patient location within the shortest possible time is important for patient survival.Type B calls are emergency calls that need immediate transportation but are not life-threatening, and type C calls require non-emergency transportation.The arrival rate for each type of call from each zone is assumed to Tabu search heuristic McLay [44] Commercial solver Schmid and Doerner [58] Variable neighbourhood search Knight et al. [35] Iterative approach Degel et al. [24] Commercial solver Van Den Berg and Aardal [67] Commercial solver Leknes et al. [38] Commercial solver Yoon and Albert [72] Branch and Benders cut Kaveh and Mesgari [34] Improved biogeography-based optimisation Boutilier and Chan [8] Simulation and heuristic Nelas and Dias [50] Commercial solver Yoon et al. [74]

Simulation optimisation This work
Memetic algorithm be known.The planning horizon is divided into a set of periods, and each period is associated with different call arrival rates for each zone.A set of potential sites are available where ambulance stations can be located.A fleet of identical ambulances is available, and these ambulances need to be assigned to the selected ambulance stations, where multiple ambulances can be assigned to each station.Demand for ambulances from a zone can be assigned to two (primary and secondary) nearest ambulance stations.Primary and secondary stations are ranked as one and two, respectively, to indicate the priority of stations in sending an ambulance to the demand zone.If all ambulances at the primary station are busy, the demand will be satisfied from the secondary station.Thus, the proposed problem consists of three simultaneous decisions related to ambulance planning: (i) to select the optimal sites for ambulance stations from the available potential sites, (ii) to assign demand for ambulances from each zone to primary and secondary stations, and (iii) to allocate ambulances to each station.The objective function of the model is formulated to maximise the sum of the survival function for type A calls and the coverage for type B and type C calls for all ambulance station locations.

Assumptions
The proposed problem makes the following assumptions to formulate the model.
(i) All ambulances are identical and can serve all call types.(ii) The complete demand for a zone is assumed to be served by ambulances from either a primary or secondary station since other stations covering the demand zone will be located considerably far away in most cases [38].(iii) The service time for an ambulance includes the time taken to return to the base station, travel time from the station to the patient location, then to the hospital, and time spent at the patient location (On-scene time).(iv) To simplify the model, the maximum number of ambulances assigned to a station is assumed to be the same for all stations.However, it can be easily relaxed by defining the input parameter  max separately for each location  as   max .

Mathematical formulation
In this section, we formulate the proposed problem as an MINLP model to determine the optimal location of stations, assignment of demand zones to stations and ambulance allocation to stations by maximising the weighted sum of survival probability and coverage for all stations.The non-linearity of the model arises due to station-specific service rates and busy probability.The non-linear equations are then linearised in Section 3.2.4 using Special Ordered Set 2 (SOS2) type variables.

𝐼
The objective function in (3.1) maximises the total performance measure value of the station locations in all periods.The performance measures are defined as a function of the mean travel time   from station  to zone .For type A calls, a survival function-based performance measure given by the formula is used [25].However, any similar monotonically decreasing function can be utilised, as the optimal location of stations is insensitive to the parameters of the survival function [30].For type B and type C calls, a coveragebased binary performance measure is used, i.e. () = 1 or 0, depending on whether   is less than or greater than the response time threshold defined.The value of  ℎ in the objective function is calculated by multiplying () by a weightage value of 4 for type A calls and 1 for call types B and C to prioritise type A calls [35].The performance weight values ( ℎ ) for all three types of calls are summarised in Table 2.

Constraints
Constraints (3.2) and (3.3) ensure that the number of stations and ambulances allocated in any period is less than the maximum number allowed, respectively.Constraint (3.4) limits the number of ambulances within the maximum limit at each station.Constraints (3.5)-(3.9)are coverage constraints, where constraint (3.5) ensures that one primary station is assigned to every demand zone.Constraint (3.6) ensures that demand can be allocated to a station with rank 1 only if it is the primary station for the corresponding zone during that period.Constraint (3.7) ensures that the primary and secondary stations for any demand zone are different in every period, while constraint (3.8) ensures that the primary station has a higher demand proportion allocated compared to the secondary station.Constraint (3.9) ensures that the total demand from each zone is covered in each period, i.e. the total proportion of demand allocated to both primary and secondary stations is equal to 1.
Equation (3.10) represents the aggregate arrival rate associated with station  during period .In equation (3.10),   is used to calculate the arrival rate for a primary station since all calls arrive at the primary station first.In contrast,  2 is used for a secondary station since only calls diverted from the primary station arrive at the secondary station.Equation (3.11) represents the aggregate service rate associated with station  during period .Equation (3.12) represents the busy probability of ambulances at a station as a function of arrival rate, service rate and number of ambulances at that station during each period.The value of   is calculated by assuming each station as an  / /-loss queueing system using the Erlang loss formula.Constraint (3.13) ensures that the probability of availability of ambulances at a station should be greater than the proportion of demand allocated to the station.Constraints (3.14) and (3.15) represent the balance between the number of ambulances being relocated to different stations from each station from one period to the next.A dummy station is considered to account for the unused ambulances from one period to another.

Linearisation of the formulation
The formulation presented above has two non-linear equations, equations (3.11) and (3.12), which are linearised using SOS2 type variables [38,68].The notation of parameters and variables used for linearisation is presented in Table 3. To The objective function (3.1) subject to constraints (3.2)-(3.10),and constraints (3.12)-(3.30),along with constraints (3.32) and (3.33), present the complete mathematical formulation of the proposed ambulance location problem.The presented model incorporates time-dependent variation in demand and allows the relocation of ambulances from one period to another to tackle the variation in demand.Another key feature of the proposed formulation is incorporating the server-level busy probability of ambulances to accurately capture the impact on actual coverage and survival probability of patients.However, this makes the model non-linear, which is linearised with the help of SOS2 type variables and breakpoint values.The linearised model consists of the linear objective function and constraints, making the model solvable by commercial MILP solvers.However, the problem is complex because of a large number of binary variables and integer variables used to represent the selection of ambulance stations, the primary station of demand zones, the number of ambulances allocated and the number of relocations.Additionally, the linearisation of constraints adds many SOS2 type variables, a collection of binary variables with a special structure.Thus, the proposed problem is more realistic and complex.
As both the objective function and constraints are linear, the problem remains feasible for all positive values of demand and travel time.However, the linearisation of busy probability and service rate introduces additional artificial variables and parameters that need to be determined to maintain the feasibility of the model.The breakpoint parameter   represents the numerator of the service rate.Therefore, the range of the parameter needs to be chosen such that the problem is feasible, as given in equation (3.34).Similarly, the range of the parameter   is given in equation (3.35).The parameter   represents the total arrival rate for each station in each period and therefore has the same range of values as   , as shown in equation (3.36).The parameter   represents the service rate of the ambulance stations for each period, given by the range in equation (3.37).
Within the above ranges of the parameters and limits of the assumptions highlighted in Section 3.1, the proposed model remains feasible and can be applied to any EMS system with a given set of demand zones and potential sites for location of ambulance stations.

Solution methodology
Our initial experiments show that the problem is not easily solvable for medium and large-sized instances.In our trial runs, using randomly generated data for instances with more than 30 zones, we could not get feasible solutions even after 12 h of using the commercial solver.This can be explained based on the observation that considering multiple periods in a day increases the solution space exponentially.The basic location-allocation problem of locating  servers in  available sites has a solution space of   [51].Considering  different periods increases the possible solution space to mr  .This rapid increase in solution space combined with the introduction of individual station-level busy probability and relocation-related constraints increases the complexity of the problem.Especially as the number of demand zones, potential locations and periods increases, the computational time required to obtain even a good quality feasible solution increases rapidly.Therefore, we develop a metaheuristic-based solution approach to solve the problem.Several authors [43,76] have applied various metaheuristic approaches to solve ambulance location problems.In this work, we adopt a MA-based approach, which improves the GA by embedding a local search routine.
MA is a population-based metaheuristic that combines natural adaptation with individual learning acquired by the members of population [36].MA expedites the process of finding the global optima and avoids premature convergence by balancing the exploration of the search space using GA while enabling the exploitation of the current neighbourhood using local search.MA has been successfully applied in solving various combinatorial problems such as scheduling [70], vehicle routing [47], assignment problems [46], supply chain network design [69] and location routing [1,48].MA has also been used in solving healthcare planning problems such as ambulance routing [28,75] and home healthcare planning [32].Memetic algorithms presented in the literature generally use binary or integer coding for solution representation.Real-coded memetic algorithms have also been presented to solve optimisation problems with continuous variables [36,39].Like a location routing problem that combines strategic level facility location problem [19,20] and tactical level routing problem [61], an ambulance locationallocation problem combines strategic and tactical level decisions.The problem presented in this work involves a combination of decisions, including station location, demand allocation and ambulance allocation.Therefore, we adopt a mixed coded solution representation in the proposed algorithm to combine binary coding for station location and real coding for demand allocation decisions.However, mixed coding using a combination of binary and integer representation has been used in MA for location routing problems [26,71].
The overall solution framework based on the memetic algorithm is presented in Figure 1.Initial solutions are generated using a heuristic approach that randomly allocates demand from different zones to stations, ensuring the diversity of solutions.Individual solutions are evaluated based on the fitness function corresponding to the value of the objective function of the problem.Initial solutions are ranked based on their fitness values, and a selection procedure is applied to choose parent chromosomes that are subject to crossover and mutation to produce offspring.A local search is applied to the offspring solutions to improve the quality of these solutions.A pre-specified set of such offspring is generated during an iteration to generate a pool of offspring by repeatedly selecting a pair of parents from the initial population.These offspring replace the worst-ranked solution among the current generation.The complete procedure is repeated until a fixed number of iterations or convergence is reached.Subsequent sections detail the major phases involved in the application of MA.

Chromosome representation of a solution
One of the most critical points in the design of MA is the representation of each solution in the form of a chromosome.A mixed coded representation is proposed where each chromosome depicting a solution is represented as a combination of two vectors, namely station vector and demand vector, to accommodate both continuous and discrete variables in the model.The station vector S(x) represents the station assigned to each demand zone (which is a discrete variable), and demand vector D(x) represents the fraction of demand from a zone that is being served by the corresponding station (which is a continuous variable).The number of ambulances allocated to a station depends on the total demand allocated to that station.The proposed chromosome representation is illustrated in Figure 2. The first element in the station vector represents that station 2 is the primary station and station 3 is the secondary station for zone 1 during period 1.Similarly, station 2 is the primary station and station 5 is the secondary station for zone 1 during period 2. The demand vector represents that 85.5% of calls from zone 1 during the first period are expected to be covered by primary station 2, while 14.5% of calls are covered by secondary station 3.

Selection, crossover and mutation
Selection is an important phase of the MA process.It is necessary to select good chromosomes from the existing population to generate good quality offspring while maintaining diversity in the solution pool to explore a larger solution space.Various selection procedures exist, such as tournament selection, roulette wheel selection, rank selection, and Boltzmann selection [53].The procedure employed in this work is adapted from the ranking selection procedure suggested by Reeves [54].We first order all the individuals from the current population in ascending order of their objective function value, considering their feasibility.The solutions with no constraint violations are ranked higher than those with violations (i.e., infeasible solutions) irrespective of their objective function value.The probability of selection of each individual is then given by the formula  () = 2  ( +1) , where  is the position of the chromosome in the ordered list of objective values, and  is the population size.The advantage of this procedure is that while being simple to implement, it gives a higher preference to the solution with better objective function value and lesser constraint violations.The solution with the best fitness value has  times higher probability than the lowest-ranked solution and two times the probability of the median value.
The two individual chromosomes selected during the selection phase are recombined to produce new offspring during the crossover phase.We apply a single point crossover separately on the station vector and the BLX- crossover for the demand vector.A single point is randomly chosen for the station vector, and the parent chromosomes are swapped along this point to produce station vectors for the child chromosome.The BLX- crossover for the demand vector part of the chromosome is performed as follows.
Step 1.Consider  1 and  2 to be the two demand proportion values in the th location of parent 1 and parent 2, respectively.Step 2. Determine  min = min{ 1 ,  2 } and  max = max{ 1 ,  2 }.
Step 5.If the random number lies outside the range of the demand proportion value, i.e. [0, 1], then select another random number.
The mutation operator in MA is used to prevent the solution from converging to local optima by reducing the convergence rate.In our problem, we apply two separate operators for mutation on the station and demand vectors.For the station vector, we randomly remove one station location from all existing station allocations and replace it with a station chosen from the set of potential locations that are not allocated.This allows the introduction of a new station allocation into the solution, thus diversifying the solution.For the demand vector, we simply replace the existing value of a given chromosome element with a random value in the feasible range for that element.

Local search
Local search is an optimisation technique that tries to generate a local optimal solution by exploring the neighbourhood of a given solution.A neighbourhood of a solution is explored by performing a few modifications to the initial solution.We apply multiple neighbourhood exploration procedures to explore the neighbourhood more efficiently and avoid rapid convergence to some local optima.Different neighbourhood operators applied in the proposed algorithm are detailed below.

For station vector
Swap stations: we select two random elements representing two zones on the station vector and swap the stations allocated at these locations on the chromosome to generate a new solution.
Replace: in this operator, we select an element on the station vector, and it is replaced by a station that is allocated to at least one of the other zones.

For demand vector
Since the demand vector is a continuous variable, a fixed step random search is applied to determine the local optima in this case.Once the local search is applied to the station vector, we change an element of the current solution by a small step in a random direction to generate a new solution.If the new solution is better than the original one, it is updated, and the next improvement is applied.

Computational results
Computational experiments were carried out on the dataset generated based on the urban location of Mumbai in India.The results were used to analyse the performance of the proposed solution approach and the impact of temporal variation in demand on the ambulance location decisions.We first describe the details related to various input parameters in Section 5.1, followed by the summary of results and related discussion in Sections 5.2 and 5.3.

Generating input data
This section presents the details related to various input parameters such as potential sites for station locations, demand, travel time, and service time for the dataset generated.

Demand zones and base locations
To make the dataset realistic, the region was divided into smaller demand zones from where the calls for ambulances arise.Potential base locations for ambulances were selected by considering various public places such as railway stations, shopping malls, public hospitals, and schools/colleges.Figure 3 shows all the different selected potential locations for the region of Mumbai under consideration.These locations were considered since the potential locations need to provide some basic features such as parking space for ambulances, security for ambulances and equipment, and electrical supply for recharging or operating equipment.The coordinates of the selected locations were obtained using the QGIS software.A subset of these potential locations and demand zones were considered to generate smaller test instances of the problem.

Call demand
We considered five randomly simulated values for total demand to create five test instances.The total demand is then divided among 144 zones such that demand is approximately proportional to the population of each zone.The smaller problem instances are generated by considering only a subset of the zones.The variation in demand data over a day is considered by breaking the day into four periods: morning (6.00 am to 12.00 pm), afternoon (12.00 pm to 06.00 pm), evening (6.00 pm to 12.00 am) and night (12.00 am to 06.00 am).The total demand is divided into three types, type A, type B, and type C calls, based on their requirements for response time within which they need to be covered.The demand was also separated by the day of the week as weekdays and weekends, as weekends tend to have slightly lower demand and less pronounced peak demand.

Travel and service time
The mean service time needed to serve any demand location from an ambulance station depends on the travel time from the station to the patient location, on-scene time, travel time from the scene to the nearest hospital and then back to the base station.To obtain a better approximation of the travel time, a large number of random points were simulated using the QGIS software for the location of calls.The time taken to reach these locations and then to the nearest hospital is determined based on the distance between locations.The time spent at a patient location is assumed to be a constant value of 15 min [38].The total service time was calculated by adding the on-scene time to all travel time components from leaving a station until the ambulance returns to the station.
The computational tests were performed on a PC with Intel(R) Core (TM) i5-4570T CPU @2.90 Hz and 8 GB of memory.The MILP model used for comparison was implemented using IBM ILOG CPLEX in concert technology with Java.The proposed MA was implemented using MATLAB 2018b.Parameter settings for the proposed MA and the GA used for comparison were obtained using multiple trial runs.The optimal population size was 12, with a crossover probability of 0.90% and a mutation probability of 0.05%.The selection procedure was repeated to obtain an offspring pool of 8 solutions in each iteration.A maximum time limit of 7200 s is considered for all computational experiments as both MA and GA converged within this time for even large instances.The same time limit is also used for CPLEX, as trial runs longer than 7200 s did not result in any meaningful improvement in the objective function value for medium and large-size instances.

Performance of the proposed approach compared to CPLEX and GA
To validate the effectiveness of our proposed solution approach, we compared the performance of the proposed approach with the exact solutions obtained using CPLEX and a GA-based approach.Various small-size instances were developed from our initial dataset by varying the number of demand zones from 5 to 15, station locations from 10 to 30, and considering five different total demand values for each case.Table 4 summarises the results obtained for these test instances from CPLEX, GA and MA.From Table 4, we observe that although CPLEX could find an optimal solution for most instances, the time taken to reach the optimal solution increases drastically as the problem size increases.While it only takes less than 41 s to reach optimality for all instances with five demand zones, CPLEX could not find the optimal solution within the specified time limit of 7200 s for any instances with 15 zones.Both GA and MA could find optimal solutions for some of the small-size instances.As problem size increases, both GA and MA converge to a solution within 2% of the CPLEX solution but take significantly less CPU time than CPLEX.For most instances, the difference between GA and MA is less than 1%, with MA slightly outperforming GA in all instances.
Medium-size instances of the problem were generated by varying the number of demand zones from 20 to 60, potential station locations from 40 to 124, and considering five different total demand values for each case.Table 5 summarises the results for the medium-sized instances for comparing the proposed approach with CPLEX and GA.From Table 5, we observe that CPLEX was able to find feasible solutions within the time limit for only three instances out of the total 25 instances, while both MA and GA were able to find feasible solutions for all instances.As there is inherent uncertainty in MA and GA, five trials of each instance were solved.The average solution value and the best solution out of the five trials are reported in Table 5.Similarly, the table also reports the percentage improvement in the best objective value found using MA compared to the exact approach and GA.The proposed MA-based approach found a slightly better solution than GA for almost all instances, with the difference varying in the range 0-3%.
Large-size problem instances were considered with more than 80 demand zones and 150 potential locations.Table 6 summarises the results for the large-size problem instances reporting the best solution found and the average of five trials found using MA and GA.CPLEX is not included in the comparison as it was unable to find a feasible solution for any of the instances within the given time limit.Among both GA and MA, there was a clear difference in the best solution found, as the memetic-based approach consistently converged to a better solution than GA.The difference between both algorithms varied up to 15%, with MA outperforming GA in 21 out of 25 instances.This difference is significant compared to the difference observed between both the algorithms in medium and small-size instances.Thus, MA outperforms GA for both medium and large-size instances in the best solution found, and the improvement is also higher as the problem size increases.Figure 4 presents a typical convergence behaviour of MA compared to GA for one of the instances where both algorithms achieved the optimal solution provided by CPLEX. Figure 4 shows that MA converges to the best solution in significantly fewer iterations than GA.A similar difference in performance was observed in other instances for both algorithms.This difference could be attributed to the local search embedded in MA, which increases the CPU time for each iteration but improves the objective function value.The increase in time required can be observed in Figure 5, which shows the convergence of GA and MA with respect to CPU time for the same problem instance.It is observed that even though MA takes almost 67% fewer iterations to achieve the same objective value but the difference in CPU time for both algorithms is negligible.

Impact of temporal variation in demand
To illustrate the impact of temporal variation in demand, we compare results obtained from three different models that locate ambulances based on three scenarios.Model 1 takes into account the temporal variation in demand and determines the optimal locations for stations and the allocation of ambulances to each station in each period.In model 2, the call arrival rate is considered equal to the average call arrival rate in a day.In model 3, we assume the same call arrival rate as in the peak period during all periods.Both model 2 and model 3 determine a single set of station locations and ambulance allocations throughout the day.The comparison of   these models is carried out based on the number of stations and ambulances located and the variation in the busy probability of ambulances.
The summary of variation in calls considered across each period for different patient types is presented in Table 7. Maximum demand occurs during the day, and minimum demand is considered during the night on weekdays and weekends.Overall demand is lower during weekends compared to weekdays.Figure 6 plots the percentage difference between the actual call arrival rate and the average call arrival rate in each period.A larger percentage difference indicates that the expected call rate will be higher than the planned demand, which will lead to a decrease in ambulance availability.It can be seen that the maximum difference between the average and actual call arrival rates occurs during weekday afternoons (12.00 pm to 6 pm), as the overall demand is highest during this period.Almost 25% of the calls could be underserved during this peak demand period, while the expected overall under-coverage is about 12% of the total calls across all periods.Table 8 presents the number of stations and ambulances located using model 1 for the input demand given in Table 7.As expected, the estimated number of ambulances required varies with the variation in demand, and the highest number of ambulances is allocated during the peak period.Also, the number of stations located follows a pattern similar to the variation in the number of ambulances.The maximum number of stations allocated is 51, with the allocation of 80 ambulances to these stations, while at the lowest demand period, 47 stations and 66 total ambulances are allocated.Variations in the number of ambulance stations can be implemented in the EMS system using temporary station locations during peak demand periods [55].Under this approach, ambulances can initially be located at permanent stations and then relocated to a temporary station to manage   observed that model 2 again results in an underestimation of ambulances during the peak period, while model 3 locates almost the same number of ambulances as in model 1 during the peak period.Another important observation is that model 1 allocates significantly fewer ambulances during low-demand periods than the other two models.Figure 8 also shows the variation in the average arrival rate of calls for each period.The number of ambulances allocated follows a pattern similar to the variation in call arrival rate, indicating model 1 assigns ambulances proportional to the demand arriving during each period.Figure 9 depicts the comparison between the mean busy probability obtained using models 1 and 2, along with the percentage difference between both values.A positive difference in the figure indicates that model 1 has a higher busy fraction, and a negative difference indicates otherwise.The busy fraction was calculated assuming each station to be an  / /-loss queueing system.It is observed that model 1 results in a lower mean busy probability value during peak periods while slightly higher busy probability during the periods of low demand compared to model 2. Figure 10 shows the comparison between models 1 and 2 based on mean server utilisation across different periods.A pattern similar to that of busy probability is observed with slightly lower utilisation for model 1 during higher demand periods and higher utilisation during lower demand periods.This is expected as model 1 allocates fewer ambulances during the low demand period and more ambulances during other periods.This shows that using average demand throughout the day will result in significantly higher utilisation of ambulances and busy probability during peak periods.As maximum demand occurs during peak periods, this will lead to increased total lost calls.
Figure 11 shows the comparison of busy probabilities for models 1 and 3.As model 3 assumes a maximum call arrival rate for all periods, a higher number of ambulances are allocated, which results in a lower busy fraction during all periods.This low busy fraction comes at the expense of lower resource utilisation, as shown in Figure 12. Figure 12 compares the mean server utilisation for both cases throughout the day, and it is clearly seen that the utilisation is significantly higher in all periods for model 1.Thus, locating ambulances assuming maximum demand as in model 3 definitely provides a good estimate of the maximum number of stations and ambulances required with lower busy probabilities.However, taking into account temporal variation leads to better resource utilisation with a slightly higher busy probability during the lower demand periods, such as night and weekends.This is an important consideration since keeping all ambulances operational during low-demand periods will require additional ambulance crew members, resulting in significantly higher operational costs.
From Figure 8, it can be observed that model 1 requires 67 ambulances to be located during weekday nights, while model 3 results in 81 ambulances, which is a significant difference.Thus, model 1 with temporal variation provides a balance between ambulance availability during peak demand periods and resource utilisation during low demand periods.Table 9 presents the number of stations and ambulances located using model 1 during the peak demand period compared to model 2 (average demand) for large problem instances.It also presents the percentage difference between both cases.The table shows that the number of ambulances and ambulance stations located is significantly higher for model 1 during peak demand.Model 3 is not considered in comparison since the ambulance locations and the number of ambulances obtained using model 3 are similar to values obtained during the peak period using model 1.The number of ambulances required during peak demand is 3.5-15% higher than that determined assuming average demand (model 2), with an average difference of 9.6% between both cases.Table 9 also shows the percentage difference between the maximum demand and average demand values for each instance.It is interesting to note that the difference between the numbers of stations located for both cases was significantly lower and varies in the range 1-5% in most instances, with an average difference of 3.6%.Therefore, we can conclude that temporal variation in demand significantly affects the number of ambulances required, although the impact is lower on the number of stations required.
Based on our results, we can make various key conclusions that can provide insights to decision-makers (DM) of EMS systems.Our model can be applied to EMS systems where multiple ambulances are located at the same station.The model allows DM to determine the optimal sites for the ambulance stations and the number of ambulances allocated to each station.Additionally, DM can also determine the preference order of stations for each zone, i.e. the rank of stations.Overall, the proposed model allows DM to determine the best allocation of available ambulances to maximise the coverage and survival probability of patients of different types.Our results show that considering temporal variation in demand provides a better estimate of the number of ambulances required compared to considering average demand.At the same time, considering temporal variation results in increased utilisation of ambulances compared to considering only the peak demand.Therefore, DM should not only focus on average demand or maximum demand over the day but also consider the hourly variation while estimating the ambulances required.An important observation from our result is that the number of ambulance stations located and the optimal location of the stations do not significantly vary over the day.However, the number of ambulances allocated to stations varies up to 15% with the variation in demand.Thus, DM can fix the optimal location of stations while varying the number of ambulances to balance the total coverage achieved and the number of relocations.
Busy probability of ambulances is affected by the variation in demand and is observed to be high during the afternoon and evening compared to the night and morning.Therefore, DM can utilise the proposed model to solve for the different maximum number of ambulances to determine the optimal number of ambulances that minimise the overall busy probability.Further, the proposed model provides the number of ambulances required during each period, which allows DM to determine the number of EMS crews required to be assigned to each station.This also allows DM to determine the schedule of the crews to operate the required number of ambulances during each period.Ultimately, we can say that our study will enable DM in the decisionmaking process to determine the optimal allocation of ambulance stations, ambulances, and the preference order of stations that maximise the performance of an EMS system while improving busy probability and server utilisation.Additionally, our model provides a framework for the operational level decisions of relocation of ambulances and scheduling of ambulance crews.We also demonstrate the effectiveness of the MA compared to the exact solver and GA for solving different size instances of the proposed problem.The exact solver is unable to provide feasible solutions for the problem within a reasonable time, while the memetic algorithm converges to a good quality solution within significantly less computational time compared to both the exact solver and GA.

Conclusions
In this paper, we present a location model for ambulances that accounts for the time-dependent nature of EMS demand along with server-level busy probability.An MINLP model is developed that accounts for the spatio-temporal variation in demand while considering heterogeneous performance measures, namely coverage and survival function, for different call types.The presented non-linear model is then linearised using SOS2 variables and breakpoint values to transform into a MILP model.We divided the day into four blocks of six hours with different call arrival rates to account for temporal variation with the minimum call rate during the night (12.00 am to 6.00 am) and maximum call rate during the day (12.00 pm to 6.00 pm).Similarly, the variation in demand during different days of the week is also considered by dividing a week into weekdays and weekends with different call arrival rates for both cases.To test the model, we generated various test instances by dividing an urban region into different demand zones and identifying potential station locations within the region.A memetic algorithm-based solution approach is proposed since commercial MILP solver (CPLEX) could not even find feasible solutions for most medium and large-size problem instances within a reasonable time.
Our results indicate that the solution generated by the proposed memetic algorithm is within 3% of the optimal solution for all the instances for which CPLEX was able to solve to optimality and with significantly less run time than CPLEX.The proposed MA-based approach also outperformed GA in terms of the best objective value for most test instances.We compared the results obtained using the time-dependent model with two models considering average and peak call rates.Computational results show that neglecting time-dependent variation and considering an average call arrival rate can result in an underestimation of ambulances required by up to 15% during peak demand.This can lead to a loss of significant demand due to a higher busy probability during the peak period while there is excess capacity during the low demand periods.Contrarily, considering a peak arrival rate for all periods provides a better estimate of the maximum number of stations and ambulances required and results in lower busy probability but has significantly low resource utilisation for all periods.Thus, considering temporal variation balances ambulance availability during peak periods and resource utilisation for the periods with low demand.Our results illustrate the need for considering the time-varying nature of demand and also demonstrate the effectiveness of the proposed solution approach.
Although this work addresses various research gaps in ambulance location models to develop a model that is applicable for more realistic and complex EMS systems, there are some limitations of our work that can be considered for future research.In this work, we have not considered the impact of variations in demand due to seasonal factors, which could play an important role in locating ambulances.Uncertainty in travel time, especially in urban regions, is another aspect that can be explored.Incorporating uncertainty in data through fuzzy logic or grey numbers and leveraging big data to model stochastic parameters can also be applied to develop decision support systems [18,63].Equity is another important issue in EMS planning, as heterogeneity in demand can result in an unequal allocation of resources.Incorporating equity-based objective function along with coverage and survival function-based objectives can be an interesting research direction.Developing multiobjective models with different performance measures and solution approaches to obtain Pareto solutions can be another possible direction for further research [21,59,60,62].Integrating the strategic and tactical level ambulance location problems with operational level problems, such as routing and relocation, is necessary to develop a complete framework for EMS planning decisions.Our model assumes that all ambulances are homogeneous, but the EMS systems generally utilise multiple types of ambulances.Evaluating the impact of various types of EMS vehicles is also an important aspect as it requires different dispatch policies for individual patients.Design of experiments based approach can be applied to fine-tune the parameters of the proposed solution approach.Another possible direction of further research is developing an exact solution approach exploiting the problem structure to improve the solution quality.

Figure 3 .
Figure 3. Region of interest with potential locations for ambulance stations.

Figure 4 .
Figure 4. Comparison of the convergence behaviour of GA and MA.

Figure 5 .
Figure 5.Comparison of GA and MA based on running time and best fitness value achieved.

Figure 7 .
Figure 7. Variation in the number of stations located by model 1 compared to models 2 and 3.

Table 1 .
Comparison of selected relevant articles with the current work.

Table 2 .
Performance weight ( ℎ ) values for different types of calls.  Service rate at station  during period    Number of ambulances relocated from station  to station  during period  3.2.2.Objective function Set of demand zones,  ∈   Set of potential stations for ambulances,  ∈   Set of time periods,  ∈   Set of ranking of stations,  ∈  = {1, 2}, 1 = primary and 2 = secondary  Set of different types of calls, ℎ ∈  Parameters  ℎ Number of calls associated with call type ℎ from demand zone  during period   ℎ Performance weight related to demand zone  if covered by station  for call type ℎ  total Total number of ambulances available  max Maximum number of stations that can be located  max Maximum number of ambulances that can be assigned to a station   Number of calls per unit time received from demand zone  during period    Travel time for an ambulance at station  to reach demand zone    Service time to serve a call from demand zone  from an ambulance at station  Decision variables   1 if a station is located at location  during period , 0 otherwise   Number of ambulances allocated to station  during period    Proportion of demand from demand zone  covered by station  with rank  during period    1 indicates station  is the primary station for demand zone  during period , 0 otherwise   Number of calls received per unit time at station  during period  ∈ ∑︁ ∈ ∑︁ ∈ ∑︁ ∈ ∑︁ ℎ∈  ℎ  ℎ   .

Table 3 .
, which represents aggregate demand at station , is divided into  breakpoint values given by   = {︀  1 ,  2 , . . .,  || }︀ , where  1 and  || are the minimum and the maximum possible values for aggregate demand.Similarly, the denominator ∑︀ ∈ ∑︀ ∈       , which represents aggregate service time, is divided into  breakpoint values represented by   = {︀  1 ,  2 , . . .,  || }︀ .Then,   and   are defined as two SOS2 variables to approximate the numerator and denominator, respectively.Defining the variable   as SOS2 is implicitly equivalent to defining ∑︀ ∈   = 1, ∀ ∈ , ∀ ∈  , and at most, only two consecutive values of   can be non-zero for each station  in period  [68].Suppose for any station  and period  if the value of ∑︀ ∈ ∑︀ ∈     = Ê, then two consecutive values   ′  and  ( ′ +1) will be chosen such that   ′   ′  +  ( ′ +1)  ( ′ +1) ≥ Ê using constraint (3.16) to approximate the numerator of equation (3.11).Similarly, constraint (3.17) together with constraint (3.23) approximates the denominator of equation (3.11) using   .Equations (3.18) and (3.19) integrate both the SOS2 variables using the breakpoint variable   .Constraint (3.20) ensures that the breakpoint variable   = 1 if station  is selected and 0, otherwise.Finally, equation (3.21) calculates the value of   using the breakpoint value.Thus, the non-linear Summary of linearisation parameters and variables.Busy probability associated with breakpoints  and , when  number of ambulances are allocated at a station  SOS2 variable related to breakpoint  for service rate linearisation at station  during period   SOS2 variable related to breakpoint  for service rate linearisation at station  during period   Breakpoint variable associated with service rate linearisation at station  during period   SOS2 variable related to breakpoint  for arrival rate linearisation at station  during period   SOS2 variable related to breakpoint  for service rate linearisation at station  during period   Breakpoint variable associated with available probability linearisation during period   { 1 ,  2,...,   } is SOS2 ∀ ∈ , ∀ ∈  (3.29) { 1 ,  2,...,   } is SOS2 ∀ ∈ ,  ∈ .(3.30)For equation (3.12), we first divide the possible service rate and arrival rate values into  and  breakpoint values given by   = {︀  1 ,  2 , . . .,  || }︀ and   = {︀  1 ,  2 , . . .,  || }︀ , respectively.Next, we calculate the busy probability   of ambulances for each breakpoint  and  using equation (3.31), assuming  / /-loss queueing system.  and  is the number of ambulances.Then, constraints (3.25) and (3.26) introduce two SOS2 variables   and   for approximating arrival and service rates, respectively.Equations (3.27)-(3.29)combine the SOS2 variables   and   using   .The relationship of variable   with the number of ambulances allocated is defined in constraint (3.32).Then, constraint (3.13) is expressed in the form of constraint (3.33).The term   ensures that constraint (3.33) holds only when  ambulances are located at station  during period .

Table 4 .
Comparison of the proposed approach with exact approach and GA for small-size instances.

Table 5 .
Comparison of the proposed approach with exact approach and GA for medium-size instances.

Table 6 .
Comparison of the proposed approach with GA for large-size instances.

Table 7 .
Distribution of total calls across periods and different call types for a sample instance.

Table 9 .
Comparison of stations and ambulances located for large problem instances.