A probabilistic model applied to emergency service vehicle location

European Journal of Operational Research 196 (2009) 323–331

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

O.R. Applications

A probabilistic model applied to emergency service vehicle location P. Beraldi, M.E. Bruni * Dipartimento di Elettronica, Informatica, Sistemistica, Università della Calabria, Laboratorio di Ingegneria per le Decisioni in ambito Sanitario, Italy

a r t i c l e

i n f o

Article history: Received 9 May 2007 Accepted 25 February 2008 Available online 8 March 2008 Keywords: Stochastic programming Location Emergency medical service

a b s t r a c t This paper is concerned with the formulation and the solution of a probabilistic model for determining the optimal location of facilities in congested emergency systems. The inherent uncertainty which characterizes the decision process is handled by a new stochastic programming paradigm which embeds the probabilistic constraints within the traditional two-stage framework. The resulting model drops simplifying assumptions on servers independence allowing at the same time to handle the spatial dependence of demand calls. An exact solution method and different tailored heuristics are presented to efficiently solve the problem. Computational experience is reported with application to various networks. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction In recent years the design of emergency medical service (EMS) systems has attracted significant attention by the scientific community fostering several contributions. Providing an adequate service level is a struggle for many health care systems, but it is of special concern in the context of EMS. Despite the recent progress in accessibility, answer swiftness, availability of more advanced equipments, some pressing problems are worth to be further investigated: the funding shortfalls, the reduction of the high costs of preparedness, the workforce issues, to mention a few. Some authors have concentrated their efforts on the strategic planning, others on tactical support (including routing of ambulances for non-urgent calls), and on operational control (ambulance dispatch and relocation). Although it is possible to identify a number of improvement areas, successful works should proceed for subsequent refinements from the strategic level to the operational one. This paper focuses on the strategic level and aims at determining the optimal location and sizing of emergency stations so as to assure a given quality of service. Different performance measures can be considered for service level as, for example, the population covered by the EMS, the equity and accessibility of the service, the timeliness of the interventions. We measure the efficiency of an EMS by its ability to provide effective service when and where necessary. Complexity and uncertainty make the design and the management of efficient EMS systems very difficult and deserve particular attention by the system planners. In the scientific literature, the stochastic nature of the problem has been recognized by several authors according to two main approaches: * Corresponding author. Tel.: +39 0984494038. E-mail addresses: [email protected] (P. Beraldi), [email protected] (M.E. Bruni). 0377-2217/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2008.02.027

the stochastic programming paradigm, commonly through the use of probabilistic constraints, and the queueing framework. We shall focus on the first approach proposing a new stochastic programming model under probabilistic constraints for designing and planning an EMS system. The proposed model makes out from other models in the literature for various aspects. First of all, the model relaxes the assumption of servers independence. Server independence and system-wide server busy probability are simplificative assumptions used in Daskin’s maximum expected coverage location problem [16] and ReVelle and Hogan’s maximum availability location problem [41]. Even though some attempts have been made to ride out simplifying assumptions by embedding the hypercube model (see [6,20]) or the queueing theory [42] into mathematical programming models, the resulting formulations are of little practical use because of their complexity. Secondly, this paper tackles the EMS design by a novel stochastic programming paradigm which encapsulates the probabilistic constraints into the two-stage framework. The resulting model is particularly suitable to represent the problem under investigation allowing to address the system’s reliability jointly with other strategic issues as the service capacity planning, location selection, and distribution of services. Thirdly, although in the last years researchers have moved their attention from deterministic to probabilistic models, the assessment of the busy fraction, either by simulation or queueing models, is still a crucial issue. We cope with this issue by considering randomness in the emergency requests rather than in server availability. Clearly, the busy fraction is tightly related to the number of calls to be serviced, and therefore treating uncertainty on the demand seems to be a more direct way to handle the problem. The paper is organized as follows. In Section 2 we review the main literature relevant for our study. In Section 3 we formulate

324

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

our stochastic model for the location of emergency service facilities under uncertainty. In Section 4 we propose an exact specialized solution method and in Section 5 we present three heuristics able to determine an approximate solution in competitive computational times. Section 6 is devoted to the presentation and discussion of computational experiments. Conclusions are presented in Section 7. 2. Related literature The proposed model can be viewed as a variant of the classical capacitated facility location problem in a stochastic setting. In order to place our contribution in the right perspective, we briefly review two main streams of literature that may be of interest for comparison: the literature on EMS location problems, and on facility location problems under uncertainty. Clearly these two research areas are strongly related, having the location of emergency service stations inspired a significant amount of research in stochastic location theory. A variety of models has been proposed for EMS and ambulance location problems in stochastic setting. They differ in the way some important issues such as, for example, spatial and temporal coverage, congestion management, disaster preparedness, vehicle availability maximization, costs minimization, are handled. Typical recommendations provided by the different models are the numbers, locations, types and capacities of vehicle base stations. Only recently the issue of relocation and routing of the emergency vehicles has been addressed [21]. Early models dealt with the static and deterministic EMS location problem. The weakness of these formulations has been recognized in the early 1980s, when the issue of server congestion begun to be handled in a deterministic setting through backup coverage [17,25]. Probabilistic models were later introduced in order to cope with servers congestion in a stochastic setting. Most of the probabilistic EMS literature extends classical deterministic location problems additionally focusing on servers congestion. The aim is to guarantee an adequate service level accounting for server availability when demands arise. The models, wrapped around traditional covering models, address the reliability issue mainly through the use of chance constraints. In the seminal work of Chapman and White [1,12] probabilistic constraints were embedded in a classical location set covering problem in order to take into account the randomness of service availability. The model, although never implemented due to mathematical difficulties, stimulated a number of probabilistic models that we refer hereafter. The Maximum Expected Covering Location Problem (MEXCLP) developed by Daskin [15,16], simplifies the model of Chapman and White assuming independence and common busy probabilities among servers. The model maximizes the expected value of population coverage given that a fixed number of facilities have to be placed on a network. Daskin developed an exchange-based heuristic that approximates the solution of the problem for all values of the probability of a vehicle being busy. The independence assumption has been shown to be unrealistic by Batta et al. [6] who proposed some corrective terms to handle more realistic situations. The assumption of a system-wide busy probability has been relaxed by ReVelle and Hogan in [40,41] by replacing the system busy fraction with an area specific busy fraction. More recently, a site specific busy fraction (based on a workload evaluation at each service point instead of around each demand point) has been proposed in [9]. Rather than considering randomness in server availability, Ball and Lin [5] and Beraldi et al. [8] acknowledged the stochastic nature of the demand process. Even though the two models start from different considerations, they are similar in terms of recommendations that they

provide to the system planner. In particular, in [8] both the location and the dimensioning problems are faced considering the random requests separately. Starting from information on the cumulative demand, an upper bound on the failure probability of the system is derived in [5]. The case of large scale emergency has been recently addressed in [26,27]. Given the probability of occurrence of an emergency situation, three different models are proposed and several heuristic strategies are compared for the solution of the resulting facility location models. A detailed review of probabilistic and deterministic location problems can be found in Owen and Daskin [39], Brotcorne and Laporte [11,35] and Galvão et al. [20]. We should point out that most of the probabilistic models mentioned so far, capture the inherent uncertainty of the EMS functioning, but often do not consider other important aspects as the presence of different type of servers, the different priority of the emergency calls, the presence of conflicting objective functions. We should mention here the fuzzy-multiobjective covering vehicle location model by Araz et al. [4] in which three important objective functions are simultaneously considered within a goal programming model. The model considers the possibility of imprecisions in aspirations levels through the fuzzy set theory. Another goal programming model is proposed in [3] where the maximization of the expected demand covered (which represents one of the two objective functions in the model) is calculated by considering rather originally, the probability of arriving to the demand point within a specified target time. Notwithstanding, the definition of a comprehensive model, able to faithfully represents the EMS functioning, is far to be achieved. Probability distributions of stochastic parameters have been combined with results of queuing theory to examine additional aspects of EMS problems. In particular, the emergency service system is modelled as a queuing system making different assumptions on the way the servers and the queue operate. Since our focus is on proactive models rather than descriptive ones, we refer the interested readers to [22,32,33,36,44] and to the references therein. The second important stream of the literature, similar in motivation to our work, concerns facility location problems under uncertainty. For extensive reviews on facility location under uncertainty we refer the readers to [10,34,45]. Some remarks are in order. Uncertainty has been incorporated in different ways in location models: the use of scenario planning techniques is an example. We only cite some contributions on the class of twostage scenario based models because of their relationship with our model, disregarding the conspicuous research area devoted to the robustness approach of decision making in location analysis [30]. A two-stage stochastic programming problem with uncertainty in demand, selling prices, and production and transportation costs has been proposed by Laporte et al. [31]. Eppen et al. presented in [19] a two-stage formulation of the capacity expansion problem. Stochastic facility location problems have been also addressed in [18] where emphasis is given to the proposed parallel interior point algorithm. In this paper, we propose a novel approach to address uncertainty based on the introduction of probabilistic constraints in the traditional two-stage framework. Unlike previous works in the EMS area, our model presents a reliability approach new in the stochastic programming literature and relaxes the independence assumption among servers. We remark that the dependence assumption has posed severe computational difficulty in the past. The only attempt to incorporate the dependence assumption in an optimization model has been made by Ball and Lin [5] although with some weaknesses. This issue will be explored in the following section.

325

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

3. Problem statement Uncertainty represents a critical issue in the design of reliable systems. This is particularly true for EMS systems where most of the key elements to consider have stochastic nature (i.e. emergency calls, travelling times). We deal with uncertainty by means of a novel stochastic programming paradigm which embeds the probabilistic constraints within the traditional two-stage framework. This choice has been motivated by the inherent two-stage nature of the EMS design problem. In effect, the location of the facilities and the definition of the corresponding capacities, can be viewed as first-stage strategic decisions. Once the uncertainty has been resolved, tactical decisions concerning allocation of customers to facilities can be taken. In a classical two-stage formulation the system reliability is not a critical issue, in that the system is designed so to hedge against all the possible future circumstances that may occur. In other words, the reliability level imposed is of 100% and this choice might lead to the design of an over-dimensioned system often not practicable for economic reasons. The probabilistic constraints allow to relax the stochastic constraints and provide the decision-maker with a tool for evaluating different alternatives in terms of cost-reliability trade-off. According to the literature, an EMS is said to be reliable if, when a service demand arises, there is at least one vehicle ready to provide service. The ability to provide adequate coverage is tied to the balance between servers available and incoming requests. In our model we assume that the main source of uncertainty comes from the emergency calls process. Congestion occurs when the service calls arising from an area during a given time interval exceed the number of available vehicles. Let us denote by I the set of demand points indexed by i and by J the set of potential locations, indexed by j, where service facilities may be located. A candidate location j can provide service to a demand point i only if the average travelling distance T ij , between i and j, is within a given threshold value T. In this case we say that i can be covered by j. According to the de-facto national standard, we can define the set of locations j that can cover the demand point i as: Covi ¼ fj 2 JjT ij 6 Tg or equivalently the set of demand point i that can be covered by the location j: Demj ¼ fi 2 IjT ij 6 Tg. Let zj denote a 0–1 variable indicating whether a facility is located at j and xj denote the integer variable indicating the number of vehicles to be housed at station j. The maximum number of servers that can be located at j is denoted by U j . The cost parameters fj and cj denote the fixed cost for opening j and the per-unit capacity cost, respectively. Let yij denote the 0–1 variable indicating the assignment of i to j and qij the corresponding cost. We assume that uncertainty affecting the service requests is represented by a random vector dðxÞ defined on a given probability space X equipped with a r-algebra F and with a probability measure P. The two-stage EMS model can be formulated as follows: min x;z

J X ðcj xj þ fj zj Þ þ Ex ½Q ðx; z; xÞ

ð1Þ

j¼1

xj 6 U j zj zj 2 f0; 1g xj integer

j ¼ 1; . . . J;

ð3Þ

j ¼ 1; . . . J:

where Q ðx; z; xÞ ¼ min y

X

ð2Þ

j ¼ 1; . . . J;

ð4Þ

J I X X i¼1

di ðxÞyij ðxÞ 6 xj

qij yij ðxÞ

ð5Þ

j ¼ 1; . . . J;

ð6Þ

j¼1

Assuming that dðxÞ follows a discrete distribution with finite support and denoting by dis the sth-realization of the demand point i under scenario s ¼ 1; . . . ; S and by ps the corresponding probability, model (1)–(9) results in a classical two-stage problem with integer first and second stage variables. The novelty of our model stems from the replacement of the stochastic constraints with their probabilistic counterparts allowing the system manager to evaluate different solutions by varying the reliability level. From a mathematical standpoint, we get the following model: min

J J S X I X X X ðcj xj þ fj zj Þ þ ps qij yijs s¼1

j¼1

xj 6 U j zj

j¼1

ð10Þ

i¼1

j ¼ 1; . . . J;

ð11Þ

zj 2 f0; 1g

j ¼ 1; . . . J;

ð12Þ

xj integer

j ¼ 1; . . . J;

ð13Þ

yijs 6 xj

i ¼ 1; . . . I; j ¼ 1; . . . J; s ¼ 1; . . . S;

yijs 2 f0; 1g 0 P B i2Demj P@ P

i ¼ 1; . . . I; j ¼ 1; . . . J; s ¼ 1; . . . S;

dis yijs 6 xj

j ¼ 1; . . . J; s ¼ 1; . . . S;

1

yijs P 1

i ¼ 1; . . . I; s ¼ 1; . . . S;

C A P a:

ð14Þ ð15Þ

ð16Þ

j2Covi

where yijs denote the assignment variable of demand point i to service facility located at j under scenario s, and a 2 ð0; 1 represents the reliability level. It is evident that a reliability level of 1 results in the classical two-stage model. We point out that the constraint 0 1 X @ P dis yijs 6 xj ; j ¼ 1; . . . ; J; s ¼ 1; . . . ; SA P a i2Demj

can be viewed as an extension of the reliability constraint proposed by Ball and Lin [5]: Y P½DðjÞ P K j  6 1  a; i ¼ 1; . . . ; I; j2CovðiÞ

where K j is the number of vehicles housed at station j and DðjÞ represents the number of calls generated by all the demand points covered by station j, that is i 2 Demj . The upper bound is established relying on the fact that, thanks to the independence assumption among demand points, DðjÞ forms a Poisson process with an arrival rate equal to the sum of the arrival rates of all i 2 Demj . We notice that the computation of DðjÞ is valid only if the demand calls arising from every circular coverage area Demj do not overlap. When considered in a real context such an hypothesis could be acceptable only as a rough approximation. Moreover, the probabilistic constraints are imposed individually on each demand point. Although Ball and Lin claim that in arriving at the probabilistic constraint they did not assume independence among servers, it is not difficult to see that the two assumptions work in such a way that one strengthens the other: the relative independence of the servers is imposed by means of the hypothesis of empty intersection among sets of demand points covered by j 2 Covi . As evident, the absence of simplifying assumptions makes the modelling process more adherent to real-world situations having the capability of providing a more accurate description of the system. Under this respect, our model extends the work of Ball and Lin in three ways:

i2Demj

X

yij ðxÞ P 1

i ¼ 1; . . . I;

ð7Þ

j2Covi

yij ðxÞ 6 xj

i ¼ 1; . . . I; j ¼ 1; . . . J;

yij ðxÞ 2 f0; 1g

i ¼ 1; . . . I; j ¼ 1; . . . J:

ð8Þ ð9Þ

 rather than considering all the demand calls that arise in Demj for each j, we consider, by the assignment variables y, only those calls which originated in Demj are assigned to the location j;  the non-independence among servers is explicitly modelled and tackled with joint probabilistic constraints;

326

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

 the model does not require the independence assumption among demand points. Even though in normal operating condition this hypothesis is well justified, it is not satisfied in the case of large scale emergency situations. We would like to remark that the joint probabilistic framework [37] is more general than the separate chance constrained formulation [13], as imposing a reliability level on each individual point does not ensure to attain a reliability level a for the entire geographical area. Moreover, the introduction of probabilistic constraints into a two-stage stochastic integer programming problem represents an original contribution [23,24,43] regardless the specific application considered. 3.1. Deterministic equivalent reformulation Problem (10)–(16) is clearly non-convex. In order to derive a deterministic equivalent formulation, we introduce big-M constraints: min

J J S X I X X X ðcj xj þ fj zj Þ þ ps qij yijs s¼1

j¼1

xj 6 Uzj ;

j¼1

i¼1

j ¼ 1; . . . ; J;

zj 2 f0; 1g;

j ¼ 1; . . . ; J;

xj integer; j ¼ 1; . . . ; J; X dis yijs 6 xj þ ds M; j ¼ 1; . . . ; J; s ¼ 1; . . . ; S; i2Demj

X

yijs þ ds M P 1;

i ¼ 1; . . . ; I; s ¼ 1; . . . ; S;

ð17Þ

j2Covi

yijs 6 xj ;

i ¼ 1; . . . ; I; j ¼ 1; . . . ; J; s ¼ 1; . . . ; S;

yijs 2 f0; 1g; S X

i ¼ 1; . . . ; I; j ¼ 1; . . . ; J; s ¼ 1; . . . ; S;

ds ps 6 ð1  aÞ

s¼1

ds 2 f0; 1g;

s ¼ 1; . . . ; S:

Here ds are the big-M parameters that make the inequalities redundant when ds ¼ 1, whereas the fulfillment of the probabilistic constraints is driven by the knapsack constraint. In order to make the presentation clearer, it is useful to rewrite problem (17) in a compact form. Let us denote by dr a feasible solution of the knapsack constraint and by D the subset of scenarios for which drs ¼ 0; s 2 D. The problem can be trivially reformulated as follows: min hðDÞ X ps P a; s2D

ð18Þ

D # f1; . . . ; Sg; where hðDÞ is problem (17) without the knapsack constraint. The proposed model is not easy to solve. Factors that contribute to the difficulty of the problem are the large number of constraints and the presence of integer and binary variables (whose number is related to the number of scenarios). Furthermore, the capacitated facility location problem even in its deterministic counterpart, is known to be NP-complete. The following section is devoted to the presentation of solution methods (both exact and heuristic) designed so to exploit the specific problem structure. 4. Exact solution approach A straightforward approach for solving problem (18) consists in exploring all the feasible solutions of the knapsack constraint evaluating for each of them the objective function value of an integer problem with a large number of variables and constraints. It is evident that an explicit enumeration of all the feasible solutions would result computationally overwhelming even for a limited

number of scenarios. A significant reduction would be obtained by adopting a dominance criterion. Given two feasible solutions d1 and d2 of the knapsack constraint we say that d1 dominates d2 if d1 P d2 , where inequality is understood componentwise. The next subsection introduces an exact method based on the concept of dominance. 4.1. Branch and Bound method We propose a Branch and Bound scheme that generalizes the method introduced in [7], where the branching is done on the scenarios to consider in the solution, rather than on continuous variables as usual. Exploring the tree, we explicitly avoid the generation of dominated solutions which cannot provide a better objective function value. To get a flavor of our Branch and Bound scheme, we refer to a simple example with only 4 scenarios with associated probabilities f0; 408; 0; 275; 0; 173; 0; 144g and a ¼ 0; 8. The Branch and Bound tree for this example is depicted in Fig. 1. At root node we fix all the d variables to one and set the upper bound UB ¼ inf. We recall that setting all the d variables to one implies the violation of all the probabilistic constraints. At the second level of the tree we obtain jSj disjoint subsets of scenarios. Exploring the tree according to a depth first strategy, we examine node 1 obtaining from the relaxation of the partial problem a lower bound LB ¼ 2625; 771. If the lower bound of any of the nodes of the tree is greater than the current upper bound, that node is fathomed, otherwise, starting form the node, we generate jS  lastordðds ¼ 0Þj children, where lastordðds ¼ 0Þ is the position of the last d set to 0 in the node. In our example the first time that a node is fathomed is at node 3 where the probability level is satisfied (with probability 0; 855). In this case, the relaxed problem yields an integer solution and the upper bound is updated accordingly. Since node 4 fulfils the knapsack constraint, it is fathomed and included in a list L of nodes to be solved as integer. Node 5 cannot be fathomed, and thus children nodes should be generated from it. However, we observe that it is not worth generating nodes that will violate the knapsack constraint. In fact in this case lastordðds ¼ 0Þ ¼ 3 and only one child ð0100Þ could be generated with probability 0:725 6 0; 8. Node 6 is fathomed since no children are allowed to be generated from this node (jS  lastordðds ¼ 0Þj ¼ 0). Similar comments can be made for the other nodes in the tree. At the end of the tree exploration we come up with the list L of subproblems to be solved as integer. This list contains a set of solutions that are proven to be superior solutions. In fact they are all not dominated and not worse than the upper bound found in the search tree. By solving the problems in the aforementioned list, the global optimum of the problem is obtained. We recall that the continuous relaxation of each problem in the list has already been solved during the exploration of the Branch and Bound tree. Thus, we can sort the nodes to be solved to integrality according to the best bound rule. The basic scheme of the proposed algorithm is reported below: Step 0 (Initialization) Fix all the d variables to one in problem (17) and put the problem in the list C of nodes to be examined. Set UB ¼ inf, L ¼ ;. Step 1 (Stop Criterion) If the list C is empty go to Step 4, otherwise select a problem p from the list and go to Step 2. Step 2 (Solution) Solve the linear relaxation of problem p and let zp be the corresponding objective function. If zp P UB then fathom the node. Otherwise check the feasibility of the PS knapsack constraint. If s¼1 ds ps > ð1  aÞ go to Step 3, otherwise check if the solution is integer. If it is integer and better than the UB update the incumbent value. If the solution is not integer but satisfies the knapsack constraint, put problem p into L.

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

327

ROOT 1111 UB inf

1 0111 LB 2625,771

2 0011 LB 2705,379

5 0101 LB 2637,036

3 0001 LB 2706,339 INTEGER

7 1011 LB 2698,337

8 1101 LB 2526,025

9 1110 LB 2565

6 0110 LB 2636,042

4 0010 LB 2705,839

Fig. 1. Example of Branch and Bound tree.

Step 3 (Branching) Let lastordðds ¼ 0Þ the position of the last d set to 0 in the node. Create jS  lastordðds ¼ 0Þj branches and add the generated subproblems to C. Go to step 1. Step 4 (Termination) Solve all the problems in L starting form the problem with lowest bound value and update the incumbent value whenever needed.

is the minimization of the number of objects, and the sum of the weights of the subset of scenarios considered must be greater than the reliability level a. The problem can be formulated as minjDj D X

ps P a;

s2D

The Branch and Bound clearly terminates after a finite number of steps, since the number of generated nodes is finite. We observe that the number of problems to solve depends on the number of feasible solutions of the knapsack constraint. Even if the proposed approach avoids an explicit enumeration of all the feasible solutions, the number of explored nodes in the search tree can be very high. This is the main motivation behind the design of our heuristic approaches. 5. Heuristic solution approaches We have designed three heuristic strategies: the first two are constructive approaches, the third one finds a pool of non-dominated solutions which are proven to be the best among all possible solutions. We observe that all heuristics focus more on the solution of the probabilistic constraints than on the solution of the location–allocation problem reflecting the characteristic bilevel structure of the problem under study. Let D ¼ fs 2 f1; . . . ; Sgjds ¼ 0g be the set of scenarios included in the optimal solution of problem (17). The solution of the problem can be viewed as a bilevel hierarchical process. In fact, if D were known, the optimal solution of the problem would consist in solving problem (17) considering only the constraints relative to scenarios s 2 D . Unfortunately, detecting D is not a trivial task. The heuristics are aimed at finding an approximate subset of scenarios to be used as a proxy of the optimal subset D . 5.1. The Mincard heuristic The first heuristic is a greedy-adding-like strategy. By using the probability associated with each scenario, we compute the smallest number of scenarios fulfilling the probability requirement. Thus, we solve a knapsack problem where each object corresponding to a scenario has a weight equal to its probability, the objective

where j  j is the cardinality of the subset considered. We denote the optimal solution of the problem above by Dmincard . Once Dmincard has been determined, we proceed by fixing in problem (17) the d vector and minimizing the problem with respect to the x; z; y variables. In spite of its simplicity, this heuristic has been shown to perform quite well. 5.2. The Approx heuristic The second heuristic introduces a more sophisticated criterion for the choice of the scenarios to include in the solution. The strategy is straightforward: we look for a subset of scenarios of minimal weight. The weights associated to the scenarios are computed separately for each scenario using the variable splitting technique [38]. This splitting scheme introduces copies xj1 ; . . . ; xjS and zj1 ; . . . ; zjS of the first-stage variable x; z and adds simple copy constraints which have a specific meaning in the context of stochastic programming. They represent the non-anticipativity principle, which states that the first-stage decisions should not depend on the scenario that will occur in the second stage. We note that these constraints link the blocks relative to different scenarios and make the problem non-separable. Ignoring the non-anticipativity constraints, we obtain jSj different subproblems that can be solved independently. Thanks to this splitting scheme we can restate problem (18) as follows: min x;z;y;d

X

^ hðsÞð1  ds Þ

s¼1;...;S

zjs ¼ zjðsþ1Þ ;

j ¼ 1; . . . ; J; s ¼ 1; . . . ; S  1;

xjs ¼ xjðsþ1Þ ; j ¼ 1; . . . ; J; s ¼ 1; . . . ; S  1; X ð1  ds Þps P a; s¼1;...;S

ds 2 f0; 1g;

s ¼ 1; . . . ; S:

ð19Þ

328

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

This reformulation makes use of binary variables ð1  dÞ because the scenario s is satisfied (i.e. included in the solution) only if ds ¼ 0. Here ! J J I X X X ^ hðsÞ ¼ p ðcj xjs þ fj zjs Þ þ q y s

ij ijs

j¼1

j¼1

i¼1

xjs 6 Uzjs j ¼ 1; . . . J; X dis yijs 6 xjs j ¼ 1; . . . J;

Table 1 Characteristic of the test problems Problem name

jIj

jJj

Number of integer variables (FS+SS)

Rows(FS+SS)

Test1 Test2 Test3 Test4 Test5

50 100 100 100 150

25 25 50 100 150

50 þ 1250jSj þ jSj 50 þ 2500jSj þ jSj 100 þ 5000jSj þ jSj 200 þ 10000jSj þ jSj 300 þ 22500jSj þ jSj

26 þ 1325jSj 26 þ 2625jSj 51 þ 5150jSj 101 þ 10200jSj 151 þ 22800jSj

i2Demj

X

yijs P 1

i ¼ 1; . . . I;

j2Covi

yijs 6 xjs

i ¼ 1; . . . I; j ¼ 1; . . . J;

yijs 2 f0; 1g

i ¼ 1; . . . I; j ¼ 1; . . . J;

zjs 2 f0; 1g

j ¼ 1; . . . J:

Table 2 Branch and Bound results

^ for each s we Minimizing the continuous relaxation of problem hðsÞ obtain jSj values vs used to solve the following knapsack problem: min

S X

Problem

Scenarios

a

Test1

10

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

vs bs 20

s¼1 S X

bs ps P a;

ð20Þ 30

s¼1

bs 2 f0; 1g: The optimal solution of problem (20) detects the subset of scenarios Dapprox ¼ fsjbs ¼ 1g that we expect will provide a good solution. Once determined Dapprox , we fix the value of the dsjs2Dapprox ¼ 0; and dsjsRDapprox ¼ 1; in problem (17) and we minimize the problem with respect to the x; z; y variables leading to a feasible solution of problem (17) with objective value s. Further insights on the quality of the bound obtained can be drawn by a deeper characterization of the optimal solution of problem (20) in comparison with the solution of problem (18). We observe that the difference between D and Dapprox depends on the relaxation of both integrality and non-anticipativity xj ) where ^ xj and ^zj are constraints. The difference (zjs  ^zj ), (xjs  ^ non-anticipative solutions, is referred to as non-anticipativity gap, whereas the difference in the value of objective function considering or disregarding the integrality constraints is commonly referred to as integrality gap. Thus, the quality of the bound obtained depends on the sensitivity of the problem to the variation of the coefficients in the objective function. Sensitivity analysis procedures could be used to determine the tolerance for the optimality of D . Unfortunately, estimating the sensitivity of our problem is as difficult as solving the original problem. Nonetheless, we can argue from the analysis of results that problem (18) is not highly sensitive and the solution of problem (20) is often equal to D .

**

Test2

10

20

**

30

40

Test3

10

20

30

**

40

Test4

10

5.3. The Comb heuristic The third heuristic is based on the application of the Branch and Bound scheme described in Section 4 initialized with a feasible solution. In effect, the Branch and Bound algorithm could be terminated before solving the integer problems (associated with the non-dominated solutions) having a lower bound value not worse than the current upper bound. Unlike classical truncated Branch and Bound, in this case the cardinality of the list L provides a worst-case measure of the number of solutions that can improve the current objective function and, therefore, is akin to an ordinal quality heuristic, where the cardinality of the list L counts the number of solutions which might be better than the current solution. This information becomes relevant when, as usually happens, the cardinality of the list L is very low compared with the number of solutions of the knapsack constraint.

40

20

**

30

40

Test5

10

20

**

30

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

B&B time

– – – –

264 166

– –

6 99 94 1806 12 110

–

20 19 – – 351 3300

–

327 133 – – –

– 1693 18 180 97 22

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

13 13 10 8314 2919 66 831 41,639 5223 – 71,839 50,266

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

180 228 313 670 1839 3086

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

CPLEX time

5 5 5 86 21 7 1692 429 190

– 5789 18,23 4578 1789 1833 423 316 313 8534 2839 1086 – 48,99 58,23

2010 1800 – – – – 23 24 26 – – 36 1028 44,375 – 80,314 – – 516 670 600 1018 2115 3470 45,086 – – – – – 1297 – – 16,018 9968 – 50,086 – –

329

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

5.4. Considerations We observe that the solution of problem (17), once determined the value of the d vector, can be very time consuming. When this is the case, a simple inner heuristic can be performed by exploiting the solution of the relaxed problem. The basic idea is to consider the so called split customers. A customer i a defined as split customer in scenario s, if the associated variables yijs have at least one fractional value for a j 2 J. If the assignment variables are spontaneously integer, we fix them also in the integer solution, other-

wise we assign customers with the following rule: we assign each split customer in a single scenario to the nearest open station not already used up to capacity. If such station does not exist, we open a new station able to serve the highest number of split customers not assigned yet. This strategy prevents the overassignment of demand points to emergency stations and limits the increase of the first stage cost. We observe that the continuous relaxation of the full problem is completely meaningless in that the binary variables driving the scenario choice have a purely logical function. On the other hand,

Table 3 Heuristics results Problem

Scenarios

a

MINCARD MINCARD error

Test1

10

20

30

**

40

Test2

10

20

**

30

40

Test3

10

20

30

**

40

Test4

10

20

**

30

40

Test5

10

20

**

30

APPROX MINCARD time

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

0 0 0 1.57 0.016 0 15.42 0.009 0 8.52 0 0.004

2 15 16 4 4 4 139 39 166 47 41 16

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

38.02 0 0 0.02 0 0 0 0.03 0 0.86 0 0

17 2 2 12 15 22 22 12 22 34 30 33

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

2 4.3 4.5 3.4 0.08 0.1 8 0.2 8.2 5.1 2.2 6.2

APPROX error 0 0 0 1.56 0.016 0 15.39 0 0 0.7 0 0.004

COMB APPROX time

COMB time

Num. Best Sol.

7 18 18 22 22 21 135 309 68 58 4 93

2 2 2 79 19 2 1396 21 21 2010 155 36

1 1 1 1 1 1 10 8 3 10 2 3

0 0 0 0 0 0 0 0 0 0 0 0

14 28 28 58 60 60 127 58 60 190 188 178

4 12 12 58 12 20 1650 1250 18 56 57 2

1 6 6 6 0 1 6 1 0 1 1 1

5 6 5 33 41 13 34 79 76 347 486 760

0 4.2 2.4 7 0.07 3.4 2.3 0.1 5.5 1.3 2.1 4.5

7 6 6 164 151 15 198 373 400 988 132 400

8 8 7 794 209 31 171 30,857 3400 22,460 13,857 34,000

2 2 2 1 4 1 0 2 1 6 4 2

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

8.6 3.1 3.1 13.4 7.8 4.5 33 6.9 7.2 8.1 15.2 4.2

24 27 61 102 102 100 349 399 460 1147 486 446

6.1 2.9 1.2 2 2 1 2.6 2.1 5 5.3 2.1 1.4

35 94 95 134 151 137 698 773 700 2988 332 211

180 128 313 444 609 3086 3440 3857 9800 32,870 9957 3400

0 1 0 1 4 0 5 2 3 1 2 4

0.95 0.975 0.99 0.95 0.975 0.99 0.95 0.975 0.99

4.8 4.3 2.9 7.6 8 3.4 14.5 3.2 4.2

77 67 60 140 118 150 338 790 1967

3.1 1.2 1.8 5.6 2 2.7 1.9 6 7.8

99 97 90 238 219 175 670 1873 00

320 200 227 7940 2090 800 31,710 28,070 34,000

1 1 1 1 1 2 5 2 1

% error = 100(Optimal Sol.  Heuristic Sol.)/Optimal Sol.

330

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

we have experimentally observed that relaxing only partially the integrality restrictions on the variables (preserving as a byproduct the logical value of the model) does not produce any substantial reduction in the solution time. Henceforth, the inner heuristic we propose offers a compromise strategy leading to solutions of acceptable quality. 6. Computational experience This section is devoted to the presentation and discussion of the computational experience carried out to assess and validate our model. The proposed methods were coded in AIMMS 3.6 [2] and CPLEX 9.1 [14] was used as LP solver. The computational experience has been carried out on a AMD Athlon processor at 1.79 GHz. We have considered different topologies used as test problems for the Two-Stage Capacitated Facility Location problem with single sourcing [28,29] and two public data sets available from [46]. Starting from this basis, we have generated different instances by varying the number of scenarios and the reliability level. We have considered capacitated instances to better simulate a real-world situation, setting the maximum number of vehicles to house at each station equal to 4. As far as costs are concerned, we have considered randomly generated values for the variable costs cj . In particular, costs were drawn from an uniform distribution ranging between 20 and 40. The assignment costs qij have been assumed to be proportional to the distance between i 2 I and j 2 J for adherence with the closest vehicle dispatch rule. When not available in the data sets, fixed costs have been generated randomly from a uniform distribution ranging between 50 and 5000. In order to generate the scenarios, we have designed the following procedure: first we have generated the demand rate randomly for each demand point (rather than considering an equal arrival rate as common in the literature), consequently, for each demand point hourly calls were randomly generated accordingly to a Poisson distribution. Table 1 shows the size of the problems parametrized with respect to the number of scenarios. The computational experiments have been carried out with the twofold aim of assessing the performance of the Branch and Bound method and of evaluating the heuristic approaches in terms of time savings and quality. Table 2 reports the solution times (in seconds) of the Branch and Bound method and the CPLEX solver. Both evaluations have been carried out by initializing the solvers with the upper bound value found by the greedy-adding-like heuristic. As far as the CPLEX time is concerned, we decided to stop the solver if its solution time is at least four times higher than the Branch and Bound time. In that case, the column CPLEX time shows one hyphen. On the contrary if the solution time of our method doubles CPLEX time we have reported in the column B&B time one hyphen. For all the other instances, we have reported both the time spent by CPLEX and by our Branch and Bound algorithm. Table 3 reports the computational results for the three heuristic approaches. In the following we shall use the acronym MINCARD and APPROX for the first two heuristics, respectively. In Table 3, the columns labelled APPROX time and MINCARD time indicate the CPU times whereas the columns APPROX error and MINCARD error report the quality of the solution obtained. As far as the third heuristic is concerned, we have reported the cardinality of the list L which contains the solutions that are be proven to be the best. The column reporting this information is labelled Num. Best Sol. and COMB time is the time spent to find the result. From the analysis of the results we can draw the following considerations:

OBJECTIVEFUNCTION 2.680 2.660 2.640 2.620 2.600 2.580 2.560 0.900

0.950

0.975 RELIABILITY

0.990

1

Fig. 2. Cost-reliability trade-off.

 The Branch and Bound scheme is efficient in most (87%) of the considered instances. However, we observe that its performance is negatively influenced by two main factors. The first one is related to the quality of the lower bounds generated during the exploration of the Branch and Bound tree. We have experimentally observed that the continuous relaxation of problems offers in most of the instances lower bounds of poor quality. The second issue regards the solutions of the problems. Clearly, the efficiency of the Branch and Bound would benefit of the integer solutions generated along the exploration of the Branch and Bound tree. Unfortunately, the solutions associated to the nodes of the tree are almost always fractional.  The Branch and Bound algorithm behaves poorly for those instances (highlighted with ** in Tables 2 and 3) in which the number of problems to be solved as integer is high (see column Num. Best Sol.). We can observe that, especially for low reliability levels, the cardinality of the list L may be high enough to slow down the algorithm in Step 4. We remark, though, that for all these instances the APPROX heuristic is able to find good near optimal solutions in a small fraction (2.41% on average) of the time spent by CPLEX to find the optimal solution.  The APPROX heuristic takes on average more time than the MINCARD, but provides better quality upper bounds. This clearly depends on the fact that very often the APPROX heuristic detects the optimal solution D .  The CPU time of the heuristics is affected by the number of scenarios considered and by the different values of a. However, the inherent difficulty of the problem seems to have an impact on the behaviour of the proposed heuristics. To overcome this difficulty we have successfully applied the inner heuristic to the most time consuming problems (Test3, Test4 and Test5). It is worth observing that when the inner heuristic is applied the quality of bounds deteriorate, due to the use of heuristics in sequence. Nevertheless, the error remains bounded between 0% and 7.8% for the APPROX heuristic, while for the easier MAXCARD heuristic the error amplification becomes more significant.  The value of the objective function is related to the imposed reliability value. As one expects, the higher is a, the higher is the objective function value. Nevertheless, at least for the tested problems, good savings in the cost function value may be obtained even keeping high reliability levels. Fig. 2 shows the objective function value as function of a for Test1 with 10 scenarios.

7. Conclusions In this paper the general and complex problem of designing and planning emergency medical services has been modelled as a twostage stochastic programming problem with probabilistic constraints. The choice of this modelling framework can be considered very satisfactory especially considering the advantage of combining two powerful and flexible paradigms that usually are separate.

P. Beraldi, M.E. Bruni / European Journal of Operational Research 196 (2009) 323–331

An exact solution method and three heuristics have been developed to solve the problem. Extensive numerical experiments have been carried out to assess the performance of proposed methods and evaluate the quality of the heuristic solutions. Although the exact approach has been shown to be very effective in most of the instances, for some problems this method may take more time than a standard general purpose software. The reason behind this is that the continuous relaxation for this problem is very poor. Future research will focus on the study of tighter relaxations (for example, Lagrangian relaxations) combined with decomposition strategies. This is beyond the scope of this paper and is the subject of ongoing research. Finally, we note that the proposed approach can be of interest for modelling and solving more general facility location problems under uncertainty. In this respect we would like to emphasize that the EMS problem has been presented to illustrate the potential of the proposed novel framework. References [1] A. Adel, J.A. White, Probabilistic formulation of the emergency service location problem, The Journal of the Operational Research Society 29 (12) (1978) 1167– 1179. [2] AIMMS 3.6 Paragon Decision Technology B.V., The Netherlands. [3] O.I. Alsalloum, G.K. Rand, Extensions to emergency vehicle location models, Computers and Operations Research 33 (2006) 2725–2743. [4] C. Araz, H. Selim, I. Ozkarahan, A fuzzy multi-objective covering-based vehicle location model for emergency services, Computers and Operations Research 34 (2007) 705–726. [5] M. Ball, F. Lin, A reliability model applied to emergency service vehicle location, Operations Research 41 (1) (1993) 18–36. [6] R. Batta, J.M. Dolan, N.N. Krishnamurthy, The maximal expected covering location problem: Revisited, Transportation Science 23 (1989) 277–287. [7] P. Beraldi, M.E. Bruni, An exact approach for solving integer problems under probabilistic constraints with random technology matrix, Annals of Operations Research, in press. [8] P. Beraldi, M.E. Bruni, D. Conforti, Designing robust emergency medical service via stochastic programming, European Journal of Operational Research 158 (2004) 183–193. [9] F. Borras, J.T. Pastor, The ex-post evaluation of the minimum local reliability level: An enhanced probabilistic location set covering model, Annals of Operations Research 111 (2002) 51–74. [10] O. Berman, D. Krass, Facility location problems with stochastic demands and congestion, in: Z. Drezner, H.W. Hamacher, Applications and Theory: Facility Location, Springer-Verlag, Berlin, 2001. [11] L. Brotcorne, G. Laporte, F. Semet, Ambulance location and relocation models, European Journal of Operational Research 147 (2003) 451–463. [12] S. Chapman, J. White, Probabilistic formulations of emergency service facilities location problems, in: Paper presented at the 1974 ORSA/TIMS Conference, San Juan, Puerto Rico. [13] A. Charnes, W.W. Cooper, G.H. Symonds, Cost horizons and certainty equivalents: An approach to stochastic programming of heating oil, Management Science 4 (1958) 235–263. [14] ILOG CPLEX 9.1 User’s Manual ILOG Inc., CPLEX Division, Mountain View, USA. [15] M.S. Daskin, Application of an expected covering model to EMS system design, Decision Science 13 (1982) 416–439. [16] M. Daskin, A maximal expected covering location model: Formulation, properties, and heuristic solution, Transportation Science 17 (1983) 48–69. [17] M.S. Daskin, E.H. Stern, A hierarchical objective set covering model for emergency medical service vehicle deployment, Transportation Science 15 (1981) 137–152. [18] A. De Silve, D. Abramson, A parallel interior point method and its application to facility location problems, Computational Optimization and Applications 9 (1998) 249–273.

331

[19] G.D. Eppen, R.K. Martin, L. Schrage, A scenario approach to capacity planning, Operations Research 37 (4) (1989) 517–526. [20] R.D. Galvão, F.Y. Chiyoshia, R. Morabito, Towards unified formulations and extensions of two classical probabilistic location models, Computers and Operations Research 32 (2005) 15–33. [21] M. Gendreau, G. Laporte, F. Semet, The maximal expected coverage relocation problem for emergency vehicles, Journal of the Operational Research Society 57 (2006) 22–28. [22] J.B. Goldberg, Operations research models for the deployment of emergency services vehicles, EMS Management Journal 1 (1) (2004) 20–39. [23] A. Gupta, C.D. Maranas, Managing demand uncertainty in supply chain planning, Computers and Chemical Engineering 27 (8–9) (2003) 1219–1227. [24] A. Gupta, C.D. Maranas, A two-stage modeling and solution framework for multisite midterm planning under demand uncertainty, Industrial and Engineering Chemistry Research 39 (2000) 3799–3813. [25] K. Hogan, C.S. Revelle, Concepts and applications of backup coverage, Management Science 34 (1986) 1434–1444. [26] H. Jia, F. Odonez, M. Dessouky, A modeling framework for facility location of medical services for large-scale emergencies, IEE Transaction 39 (2007) 41– 55. [27] H. Jia, F. Odonez, M. Dessouky, Solution approaches for facility location of medical supplies for large scale emergencies, Computers & Industrial Engineering 52 (2007) 257–276. [28] A. Klose, An LP-based heuristic for two-stage capacitated facility location problems, Journal of Operational Research Society 50 (1998) 157–166. [29] A. Klose, A Lagrangean relax-and-cut approach for the two-stage capacitated facility location problem, European Journal of Operational Research 126 (2000) 408–421. [30] P. Kouvelis, G. Yu, Robust Discrete Optimization and its Applications, Kluwer, Dordrecht, 1996. [31] G. Laporte, F.V. Louveaux, L. van Hamme, Exact solution to a location problem with stochastic demands, Transportation Science 28 (2) (1994) 95– 103. [32] R.C. Larson, A hypercube queueing model for facility location and redistricting in urban emergency service, Computers and Operations Research 1 (1974) 67– 95. [33] R.C. Larson, Approximating performance of urban emergency service systems, Operations Research 23 (5) (1975) 845–868. [34] F.V. Louveaux, Discrete stochastic location models, Annals of Operations Research 6 (1986) 87–94. [35] V. Marianov, S. ReVelle, Siting emergency services, in: Z. Drezner (Ed.), Facility Location: A Survey of Applications and Methods, Springer, Heidelberg, 1995, pp. 199–223. [36] V. Marianov, C. ReVelle, The queuing maximal availability location problem: A model for the siting of emergency vehicles, European Journal of Operational Research 93 (1996) 110–112. [37] L.B. Miller, H. Wagner, Chance-constrained programming with joint constraints, Operations Research 13 (1965) 930–945. [38] M. Nasberg, K.O. Jonstern, P.A. Smeds, Variable Splitting – A New Lagrangian Relaxation Approach to Some Mathematical Programming Problems, Technical Report, LITH-MATH-R-85-04, Linkoping University, 1985. [39] S.H. Owen, M.S. Daskin, Strategic facility location: A review, European Journal of Operational Research 111 (1998) 42347. [40] C.S. ReVelle, K. Hogan, A reliability-constrained siting model with local estimates of busy fractions, Environment and Planning B: Planning and Design 15 (1988) 143–152. [41] C.S. ReVelle, K. Hogan, The maximum availability location problem, Transportation Science 23 (1989) 192–199. [42] C.S. ReVelle, K. Hogan, The maximum reliability location problem and a reliable p-center problem: Derivatives of the probabilistic set covering model, Annals of Operations Research 18 (1989) 155–174. [43] A. Ruszczynski, Probabilistic programming with discrete distributions and precedence constrained knapsack polyhedra, Mathematical Programming 93 (2002) 195–215. [44] F. Silva, D. Serra, Locating emergency services with different priorities: The priority queuing covering location problem, Journal of the Operational Research Society (2007) 1–10. [45] L.V. Snyder, Facility location under uncertainty: A review, IIE Transactions 38 (7) (2006) 537–554. [46] .