JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998, pp. 401–424 PROBABILISTIC, MAXIMAL COVERING LOCATION– ALLOCATION MODELS FOR CONGESTED SYSTEMS* Vladimir Marianov Department of Electrical Engineering, The Catholic University of Chile, Casilla 6177, Correo 22, Santiago, Chile Daniel Serra Department of Economics and Business, Universitat Pompeu Fabra, Trias Fargas, 25-27, Barcelona 08005, Spain, E-mail: daniel.serra@econ.upf.es ABSTRACT. When dealing with the design of service networks, such as health and emergency medical services, banking or distributed ticket-selling services, the location of service centers has a strong influence on the congestion at each of them, and, consequently, on the quality of service. In this paper, several probabilistic maximal covering location– allocation models with constrained waiting time for queue length are presented to consider service congestion. The first model considers the location of a given number of single-server centers such that the maximum population is served within a standard distance, and nobody stands in line for longer than a given time or with more than a predetermined number of other users. Several maximal coverage models are then formulated with one or more servers per service center. A new heuristic is developed to solve the models and tested in a 30-node network. 1. INTRODUCTION When dealing with the design of service networks, the location of service centers has a strong influence on the amount of congestion at each of them, and, consequently, on the quality of service. Examples of these types of systems include primary health care center services, banks, and distribution centers. Primary health care centers have one physician attending general cases who refers those cases that cannot be cured to a specialist. These centers must be located in such a way that they may be reached from any demand point within a reasonably short time. Then because waiting time is an important determinant of the perception of service quality, once a patient has arrived at the center, his or her waiting time should be as short as possible. Other examples are banks, which have branches distributed over a geographical area, where clients stand *This research has been possible thanks to grants by NORTEL—External Research, FONDECYT (the Chilean State Scientific Research Fund) and the Fundación Banco Bilbao Vizcaya, Spain. Received October 1996; revised March 1997; accepted April 1997. © Blackwell Publishers 1998. Blackwell Publishers, 350 Main Street, Malden, MA 02148, USA and 108 Cowley Road, Oxford, OX4 1JF, UK. 401 402 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 in lines, waiting to be served, or distribution centers, where trucks arrive to deliver their load. The presence of more trucks than servers at a center can easily make it congested. In this paper, several probabilistic, maximal covering, location–allocation models with constrained waiting time for queue length are presented to consider service congestion. The first model addresses the issue of the location of a given number of single-server centers such that the maximum population is served within a standard distance, and nobody stands in line for longer than a given time, or with more than a predetermined number of other clients. We then formulate several maximal coverage models, with one or more servers per service center. A new heuristic is developed to solve the models, and is tested in the case of a 30-node network. 2. LITERATURE REVIEW Existing literature has presented several models that are explicitly intended for the design of spatial queueing systems. However in general, these models are oriented to emergency systems in which servers travel to the site of the emergency, rather than systems in which servers are fixed. Some models assume a single server in the region under study. For example, Berman, Larson, and Chiu (1985) develop a heuristic algorithm to locate optimally one server on a congested network and formulate a Stochastic Queue Median. Batta (1988) considers the situation in which there may be a selective rejection of calls by the dispatcher. The model is presented together with a greedy heuristic procedure for location of the server. Batta, Larson, and Odoni (1988) present a model and an algorithm for locating one server when there are calls of different priority. Batta (1989) presents a model to study the effect of using expected service-time-dependent queueing disciplines on the optimal location of a single server. In some cases, models assuming single servers are used as building blocks for algorithms that locate more than one server. Based on the one-server location algorithm of Berman, Larson, and Chiu (1985), Berman, Larson, and Parkan (1987) develop two heuristics for locating p servers on a congested network. At each step a procedure called Mean Time Calibration is used (Larson and Odoni, 1981) which in turn includes solving the exact hypercube. This creates an algorithm suitable for systems with only a few servers. After the Mean Time Calibration is performed, the first heuristic uses the 1-median to improve the location of each one of the servers, and the second one uses the Stochastic Queue Median. All these models are nonlinear. Some authors have also considered the problem of districting, or determination of optimal service territories. Berman and Larson (1985) solve the problem of districting for a two-server network in the presence of queueing. Given the locations of two servers on a congested network, a nonlinear model and a heuristic algorithm are developed to determine the optimal service territories of each server. Each district behaves as a M/G/1 system, with FIFO © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 403 queues. In this case optimality refers to the minimum average response time to a random customer. There is no interaction between districts. Berman and Mandowsky (1986) combined the two-server districting algorithm with the Stochastic Queue Median, to develop a general location–districting iterative algorithm for two units, n-nodes and m-server networks. Initially the best districting is found for a first location choice, followed by a relocation of both servers given fixed districts. This procedure is iterated as needed. In the case of m servers, at each step two servers are optimally located and their optimal service areas found while all the remaining units stay at a fixed location. Again, there is no interaction between districts. In general, a fairly large computational effort is required if these models are to be used for the location of servers. All the models are nonlinear, heuristic, and their objective is to minimize expected response time of servers that travel to the site of an emergency. Also, all models use approximations in order to model the system. On the other hand, optimization models for location, derived from the Location Set Covering Problem (LSCP) by Toregas et al. (1974), the p-median problem (Hakimi, 1964; ReVelle and Swain, 1970), and the Maximal Covering Location Problem (MCLP) Church and ReVelle (1974), even though they are NP-hard and therefore cannot be solved in general to optimality, can be solved to optimality for reasonable size problems when linear approximations that result from relaxing the integer constraints are used. Optimization models typically include a capacity constraint to deal with congestion. This forces the demand for service at each center to be smaller than its maximum capacity. Under these constraints, the maximum capacity is usually set by the number of servers at the center, or by some estimate of the maximum number of users that the center can serve at the same time. To estimate the demand for service at any center, two different measures are often used. The first is the total population allocated to the center, multiplied by some experimental or practical factor; this gives an estimate of the expected number of simultaneous requests for service, or required workload, of that center. The second is an average of the historical rate of requests originating from the population allocated to the center, if a record is available. This demand is then constrained to be smaller than the maximum capacity of the center. Other models that deal with congestion examine the behavior of the potential users of the service (see Brandeau et al., 1995). Because the utilization of a facility depends on the decisions of all the users in the aggregate, user-choice models have at their core an equilibrium problem. The demand at each node is an independently distributed random variable and consumers patronize a facility with a given probability. The optimal location of the facilities is found given the equilibrium behavior of the users. In the traditional optimization models where demand is implicitly assumed to be constant over time and equal to an average (which is always strictly constrained to be smaller than the maximum capacity of a center) there is a contradiction between this strong, rigorous constraint, and the fact that it is © Blackwell Publishers 1998. 404 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 applied to an average that by definition, is exceeded 50 percent of the time. In this paper, we propose models based on the more realistic assumption that the number of requests for service is not constant over time, but is, instead, a stochastic process. This stochasticity of the demand is explicitly taken into account in order to derive a capacity-like constraint. Instead of imposing an upper bound on the demand equal to the capacity of the center, this forces a lower bound on the quality of the service at the center, particularly the waiting time or the number of people in line, awaiting service. We present several optimization models. Rather than minimizing the average response time, as other authors do, we restrict either the response time or the queue length to be smaller than a predetermined figure. This formulation allows us to obtain linear objectives and linear constraints. Therefore, although the problem is NP-hard, we are able to solve the location problem to optimality, if desired, using linear-programing relaxation and branch and bound (LP+BB). Although the LP+BB procedure may be burdensome in terms of computing time. In contrast to other location models in the literature, it does not require any approximations of the queueing model. A distinct contribution of the formulations presented here (as opposed to models that utilize averages, or models that minimize some quality indicator) is that a particular value of service quality of the system (as given by queue lengths or waiting times) is explicitly embedded in the optimization model. Therefore, by using these models when designing a system, the system can be fine-tuned, allowing the designer to trade off investment and operating cost versus service quality. The models assume that consumers patronize the nearest service and no considerations are made on the user’s choice with respect to the level of congestion. In addition, an infinite queue length is assumed. This assumption is standard for most queueing models, including situations where there is a (relatively large) finite upper bound on the permissible number of customers. The first model focuses on the location of a given number of single-server centers so that the maximum population is served within a standard distance, and nobody stands in line for longer than a specified time, or with more than a predetermined maximum number of other clients. We then formulate several maximal coverage models, with one or more servers per service center. 3. THE QUEUEING, MAXIMAL COVERING, LOCATION–ALLOCATION MODEL (QM-CLAM) The maximal covering, location–allocation model can be stated as: Locate p centers and allocate users to them so as to maximize covered population, where coverage is defined as: (i) covered population is allocated to a center within a time or distance standard from their home location, and (ii) if a user is covered on arrival at the center, he or she will wait in a line with no more than b other people, with a probability of at least α. © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 405 or Locate p centers and allocate users to them so as to maximize covered population, where coverage is defined as: (i) covered population is allocated to a center within a time or distance standard from their home location, and (ii) if a user is covered on arrival at the center, he or she will be served within time τ of arrival at the center, with a probability of at least α. The Church and ReVelle (1974) maximal covering, location problem (MCLP) model cannot be modified to accommodate the congestion constraint. The MCLP model does not contain allocation variables so it is not possible to allocate servers to demands covered by them. Without allocation variables it is not possible for the model to aggregate all the calls for service arriving at a server and determine when congestion occurs. Thus, we rewrite the maximal covering model as a p-median-like model, using location-allocation variables. The formulation of the model is (1) Max ∑a x i ij i, j (2) (3) s. t. xij ≤ y j ∀i, j ∑x ∀i ij ≤1 j (4) (4a) (5) P(center j has ≤ b people in queue) ≥ α P(waiting time at center j ≤ τ) ≥ α ∀j ∑y i ∀j =p i yj , xij = 0,1, j ∈ Ni ∀i,j Zero–one variable yj is one if a center is located at node j, and zero otherwise. Zero–one variable xij is one if users located at demand node i are allocated to a center located at node j, and zero otherwise. If the set Ni is defined as the set of candidate locations that are within a distance standard from node i, or the set of candidate locations which can be reached from node i within a certain standard time, variables xij are only defined for the pair of subscripts (i, j) such that j ∈ Ni. This is because we are forcing coverage of every demand from a center located within the distance standard, thus it is unnecessary to define variables that will never equal one. The parameter ai is the total population at demand node i. Equation (1) maximizes the population allocated to a center. Note that there is no constraint forcing every demand node to be covered. Constraint (2) states that it is not possible to allocate a demand node i to a node j unless there is a center at the last one. Equation (3) forces each demand node i to be allocated to only one service center j. Equation (5) sets the number of centers to be located, whereas Equation (4) forces every center to have a maximum of b people in line with a probability of at least α. This constraint ensures that on arrival at the center, a user usually will find a line no longer than b. Equation (4a) explicitly © Blackwell Publishers 1998. 406 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 makes the total time spent by a user at the center shorter than or equal to τ with probability of at least α, thus assuring every user timely service. The usual capacity constraints would be less restrictive. (4b) expected number of customers in center j ≤ b (4c) expected time spent at center j ≤ τ ∀j ∀j The first constraint, Equation (4b), forces the expected (average) number of users at center j (the one being attended plus those in line) to be less than or equal to b. If this constraint holds, on arrival at the center users will find queues shorter than b 50 percent of the time [rather than 100 × α percent, as in Constraint (4)]. Finally, Constraint (4c) ensures the average time spent in the system is less than or equal to τ. In order to write Constraint (4) in a tractable form, we make the customary assumption that requests for service at each demand node i appear according to a Poisson process with intensity fi. Each center serves a set of demand nodes, therefore the requests for service at that center are the union of the requests for service of the nodes in the set. Thus they can be described as a stochastic process, equal to the sum of several Poisson processes. It is easy to show that this stochastic process can also be a Poisson process, with an intensity λj equal to the sum of the intensities of the processes at the nodes served by the center. This set of nodes is unknown until the mathematical programing problem is solved. However, we can use variables xij in order to rewrite the parameter λj as λj = ∑f x i ij i ∈l Using this definition, if a particular variable xij is one, meaning that node i is allocated to center j, the corresponding intensity fi will be included in the computation of λj. We also assume an exponentially distributed service time, with an average rate of µj where µj ≥ λj to maintain equilibrium in the system. This is a reasonable assumption because only the service time at the center is considered and not the travel time of clients. If we assume steady state, we can use the results of an M/M/1 queueing system for each center and its allocated users. If we define the state k of the system as k users in the system (either being attended or in line), the state-transition diagram of the system is the one shown in Figure 1. Here, state k corresponds to k users at the center, that is, state k = 0 corresponds to the center being idle, state k = 1 to one user being attended at the center, state k = 2 to two users at the center: one of them getting attention and one in the queue, and so on. We want the probability of a user being in a line with no more than b other people, to be at least equal to α. If we represent as pk the steady-state probability of being in state k, this requirement is written as (6) © Blackwell Publishers 1998. p0 + p1 +....+ pb+1 ≥ α MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 407 FIGURE 1: State Transition Diagram, M/M/1 Queuing System. Writing and solving the steady-state balance equations of the M/M/1 system, we get the following expression for the steady state probabilities (Wolff, 1989) pk = (1 – ρj)ρjk where ρj = λj /µj. Hence, Equation (6) becomes: (1 – ρj) + (1 – ρj)ρj + (1 – ρj)ρj2 + ... + (1 – ρj)ρjb+1 ≥ α or b+ 1 d 1− ρ j i∑ ρ k j ≥α k= 0 which is equivalent to d 1− ρ j i 1 − ρ bj + 2 ≥α 1− ρ j or 1 − ρ bj + 2 ≥ α ρ bj + 2 ≤ 1 − α ρ j ≤ b+ 2 1 − α Then because ρj = λj /µj, (7) λ j ≤ µ j b+ 2 1 − α Equation (7) is equivalent to Constraint (4). Using the relationship between the intensity at the center and the intensities at the demand nodes, Constraint (4) is rewritten as (8) ∑f x i ij ≤ µ j b+ 2 1 − α i ∈l Equation 8 is a linear, deterministic equivalent of Constraint (4). If we choose to use Constraint (4a) instead of Constraint (4), that is P[time spent at center j ≤ τ] ≥ α © Blackwell Publishers 1998. ∀j 408 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 we may use the probability distribution function of the waiting time in a M/M/1 queue, w, which has the following expression (Larson and Odoni, 1981) fw w j = µ j − λ j e d i d i − µ j −λ j w j d i to derive its cumulative distribution (9) bg P w j ≤ τ = Fw τ = 1 − e d i − µ j −λ j τ d i The probability in Equation (9) is set at greater or equal to α 1− e e − µ j −λ j τ d i − µ j −λ j τ d i ≥α ∀j ≤ 1− α ∀j b − µ j − λ j τ ≤ ln 1 − α d i λj ≤µj + 1 ln 1 − α τ b g ∀j g ∀j Rewriting the parameter λj as a function of the variables, we get ∑f x (10) i i ∈l ij ≤ µ j + 1 ln 1 − α τ b g ∀j Equation (10) is the linear, deterministic equivalent of Constraint (4a). The final formulation of the model is as follows Max ∑a x i ij i, j s. t. xij ≤ y j ∀i, j ∑x ∀i ≤1 ij j ∑y i =p i ∑f x i ij ≤ µ j b+ 2 1 − α or i ∈l i i ∈l xij = 0, 1, j ∈ N i ∀i, j y j = 0, 1 ∀i, j © Blackwell Publishers 1998. ∑f x ij ≤ µ j + 1 ln 1 − α τ b g ∀j MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 409 If the problem is solved using its linear relaxation, in these last equations the right-hand side may be multiplied by variable yj in order to improve the integer characteristics of the variables at the solution. 4. THE QUEUEING, MAXIMAL COVERING LOCATION–ALLOCATION MODEL WITH CO-LOCATION OF m SERVERS PER CENTER This queueing maximal covering location model can be stated as: Locate p service centers, each with m servers (where m may depend on the location, that is, it could be an mj) and allocate users to them so as to maximize covered population, where coverage is defined as: (i) covered population is allocated to a center within a time or distance standard from its home location, and (ii) if a user is covered, on arrival at the center, he or she will wait in line with no more than b other people, with a probability of at least α. The formulation of the model is the same as the previous one. We develop this model using only the constraint on queue length. Elsewhere, we formulate a model with an upper probabilistic bound on the waiting time. Again, we assume an exponentially distributed service time for each center, with a service rate of µj. If the number of servers at the center is mj, the inequality µj ≥ mjλj must hold or the system does not reach an equilibrium. By virtue of the above assumptions, the results for an M/M/m queueing system can be used for each center and its allocated users. The state-transition diagram of the system is shown in Figure 2. In this figure, state k corresponds to k users at the center with each user being attended to. That is, state k = 0 corresponds to the center being idle, state k = 1 to one user being attended at the center, state k = 2 to two users at the center, both of them receiving attention, and so on up to state m in which all m users in the system are being attended. In state m + 1, however, m users are being attended while another stands in line, state m + 2 represents m users in service and two in line, and so on. We want the probability of a user finding a queue with no more than b other people, to at least equal α. This requirement is written as p0 + p1 + .... + pm+b ≥ α that is, the probability of the queue being shorter than or equal to b users at the arrival of the next request, is greater than α. Also, since p0 + p1 +....+p∞ = 1, (11) pm+b+1 + pm+b+2 + ... + p∞ ≤ 1 – α which means that the probability of the queue being longer than b is smaller than 1 – α. Note that the special case b = 0 does not necessarily mean that the user finds one server available; it may happen that all m servers at the center are busy, but there are no users in the queue. In this case, the arriving customer must wait until one of the servers becomes idle. If immediate attention is desired, that is, at least one server free with probability α, then p0 + p1 + .... + pm–1 must be made greater or equal to α. © Blackwell Publishers 1998. 410 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 FIGURE 2: State Transition Diagram, M/M/m Queuing System. Solving the steady-state balance equations of the M/M/m system, we get the following expression for the steady-state probabilities (Wolff, 1989) pk = p0 ρk/k! pk = p0 ρk/m! mk–m k≤m k>m I F G ρ ρ J + =G GGH FGH 1 − mρ IJK m! ∑ j ! JJJK m− 1 m p0 −1 j j =0 Where ρ = λ/µ. Although these parameters are specific to each server center, we will not use any subscript for the moment. Using these expressions for the steady-state probabilities, Equation (11) becomes k =∞ p0 m m ρ ≤ 1− α m! m k= m + b+ 1 ∑ FG IJ H K or p0 m m ρ m! m FG IJ H K m+ b+ 1 ∞ ρ ∑ FGH m IJK k ≤ 1− α k= 0 Since ρ/m ≤ 1, the summation converges, and can be written in a simpler form. If we also replace the expression for p0, we obtain O LM ρ ρ P MM F ρ I + ∑ j ! PP PQ MN GH 1 − m JK m! m m− 1 −1 j m + b+ 1 F FG ρ IJ m GG H m K m! G ρ GH 1 − m m j =0 I JJ ≤ 1 − α JJ K After some algebraic manipulation, this equation becomes (12) m− 1 bm − kgm! m ∑ k=0 © Blackwell Publishers 1998. k! b 1 ≥ 1 ρ m + b + 1− k 1 − α MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 411 Because ρ = λ/µ, and λ is a function of the variables xij, Equation (12) can be rewritten as a function of variables xij, thus becoming the deterministic equivalent of Equation (4). It is intuitively easy to see that for any fixed value of α, the value of the left-hand side of Equation (12) can be made large enough to make the equation hold by making ρ small enough because its exponent is always positive. The value of variable ρ is decreased by manipulating variables xij, (making as many of them equal to zero as needed). Furthermore, for any value of α there must exist a value ρα which makes Equation (12) hold as an equality, as well as a range of values of ρ such that Equation (12) holds as a strict inequality. Although it is the deterministic equivalent of Equation (4), Equation (12) cannot be used in a linear model because of its nonlinearity. However, we show next that its left-hand side (LHS) is strictly decreasing with increasing ρ, and we later use this characteristic to find a linear equivalent. LEMMA: The left-hand side (LHS) of Equation (12) strictly decreases with ρ Proof The derivative of the LHS with respect to ρ is (13) ∂LHS = ∂ρ m−1 ∑ b k=0 g bm − kkg!m! m − m + b + 1− k b 1 ρ m + b+ 2 − k This derivative is strictly negative because (m + b + 1 – k) is strictly positive, as are the remaining factors of each term of the summation. The entire summation is also strictly negative, as is the whole derivative. Thus, the LHS of Equation (12) is a strictly decreasing function of ρ, which also means that it strictly increases when ρ decreases. THEOREM: The left-hand side of Equation (12) is strictly decreasing with ρ, so if ρα is the value of ρ which makes Equation (12) hold as an equality, the equation ρ ≤ ρα , is a linear equivalent of Equation (12) when used as a constraint for a linear programing formulation Proof Let ρα be the value of ρ which makes Equation (12) hold as an equality. Since the LHS strictly increases when ρ decreases, for any value of ρ ≤ ρα Equation (12) also holds. In other words, the inequality ρ ≤ ρα is a sufficient condition for Equation (12) to hold. Furthermore, by virtue of the strictly increasing nature of the LHS of Equation (12), for any ρ > ρα Equation (12) cannot hold. Thus, ρ ≤ ρα is also a necessary condition for this equation to hold. Because ρ ≤ ρα is a necessary and sufficient condition for Equation (12) to hold, we can use it in place of the Equation. Thus, we have Equation (14) the new deterministic, linear equivalent to Equation (12). (14) © Blackwell Publishers 1998. ρ ≤ ρα 412 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 Once the value of α is given, ρα can be found by using any numeric root-finding technique (for example, Newton methods). This procedure must be repeated for each service center j, and a value ραj found for each one, obtaining the set of equations ρj ≤ ραj ∀j Since ρj = λj / µj λj ≤ µj ραj ∀j Recalling that λj is a function of variables xij ∑f x (15) i ij ≤ µ j ρ αj ∀j i ∈l is the set of linear, deterministic equivalents of Constraint (4). Remember also that ρα j is a function of the number of servers at center j (mj), which can be specified for each center before solving the model. The final model is Max ∑a x i ij i, j s. t. xij ≤ y j ∑x ≤1 ij ∀i j ∑y i =p i ∑f x i ij ≤ µ j ρ αj ∀j i ∈l xi j = 0,1, j ∈ Ni yj = 0,1 ∀i, j ∀i, j In this model, variable xi j is one if mj servers are located at center j, and zero otherwise. 5. COMPUTATIONAL EXPERIENCE: BRANCH-AND-BOUND AND HEURISTIC METHODS We now assume that the service centers are primary health care facilities. The servers are physicians, and there are one or more physicians at each center. The average service time is set at 20 minutes. We use the 30-node network, shown in Figure 3. Each demand center is also a potential server location, and the distances are Euclidean. We further assume that centers serve only patients located within a distance of 1.5 miles. The populations and coordinates of the nodes of the network are shown in Table 1. Assuming one server per center, we computed the value of the right-hand side of Equation (8) (limit values of arrival © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 413 FIGURE 3: 30-Node Test Network. rates, λ), given an average service time 1/µ of 20 minutes, different values of α, and different values of queue length b. In order to compare both criteria of quality of service (queue length and total response time), we show the total response time w (time in queue plus time in service that the clients will not exceed with TABLE 1: 30-node Test Network Node X Y Popn Node X Y Popn 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 3.2 2.9 2.7 2.9 3.2 2.6 2.4 3.0 2.9 2.9 3.3 1.7 3.4 2.5 2.1 3.1 3.2 3.6 2.9 2.9 2.5 3.3 3.5 2.7 2.1 2.8 5.3 3.0 6.0 2.8 710 620 560 390 350 210 200 190 170 170 160 150 140 120 120 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 3.0 1.9 1.7 2.2 2.5 2.9 2.4 1.7 6.0 1.9 1.0 3.4 1.2 1.9 2.7 5.1 4.7 3.3 4.0 1.4 1.2 4.8 4.2 2.6 2.1 3.2 5.6 4.7 3.8 4.1 110 100 100 90 90 90 80 80 80 80 70 60 60 60 60 © Blackwell Publishers 1998. 414 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 probability α), that would result in the same values of λ if Equation (10) is used instead of Equation (8). The values are shown in Table 2. Queue length b = 0 means that, at the time a new user arrives, there are either no other users in the system, or there is one user being attended by the server. It is interesting to note that, for a particular value of λ, total response time appears intuitively to be long compared to queue length. This is because when computing a queue length of, for example b = 5, occupation states k = 0 to k = 6 (no users in the system, up to one user being attended plus five users in line) with their associated probabilities p0 to p6, have to be considered. However, when dealing with response time, the probabilities of long waiting times are aggregated over all possible states k. Branch-and-Bound; Constrained Waiting Time; One Server per Center The commercial integer programming software CPLEX 3.0, on a cluster of eight DEC 3000–700 AXP computers was used to solve the QM–CLAM model. The branching was limited to 20,000 in each case. With this limit, the run time was between 10 and 40 minutes per run. The daily call rate for each demand node was set at 0.006 times the population. Table 3 shows the results for α = 85 percent, 90 percent, and 95 percent, and different values of the waiting time, τ. In this table, S is the number of servers to be located. A “y” in the columns “L” indicate when the branching limit of 20,000 was reached. The columns “Obj” indicate the value of the objective and the columns “Locns” the locations of the centers. Note that for short waiting times, a difference of one minute can dramatically change the coverage. This can be seen for example in the cases τ = 48 and 49, α = 90 percent. For τ = 49 min, nine servers are enough for complete coverage of the population. However, if the required waiting time is reduced by one minute, complete coverage becomes impossible, no matter how many servers are deployed. The limit of 20,000 branches was reached in many of the runs. In these cases the solution is not necessarily optimal. Note also, that there are some cases in which it is not possible to achieve total coverage of the population, no matter how many servers are located in the system, so some demand nodes are not covered by the service. This is because the call rate of these nodes exceeds the service capability of a single server. TABLE 2: Limit Values of Arrival Rates α = 0.99: α = 0.9: α = 0.5: b λ [1/min] w [min] b λ [1/min] w [min] b λ [1/min] 0 1 2 3 4 0.005 0.0108 0.0158 0.0199 0.0232 102.33 117.35 134.70 153.10 171.90 0 1 2 3 4 0.0158 0.0232 0.0281 0.0315 0.0341 67.35 85.94 105.22 124.79 144.50 0 1 2 3 4 0.0354 0.0397 0.042 0.0435 0.0445 © Blackwell Publishers 1998. w [min] 47.5 67.5 87.5 107.0 127.0 α = 85 percent τ = 40 S Obj L Locns τ = 41 Obj L τ = 42 Locns Obj L – – – – τ = 49 Locns Obj L – – – – – – – – τ = 52 Locns Obj L Locns – – – – – – – – – 9 4140 n 1,3,6,7,13, 16,23,24,30 8 4140 n 1,2,6,14,19, 24,25,30 5470 n 1,6,7,8,9, 19,22,24 7 4060 y 2,4,7,17, 22,25,30 5390 y 1,3,4,5, 10,22,25 – – – – – – 6 3500 y 1,7,10, 13,15,22 5110 y 2,4,6,7,17, 22 5470 n 5,6,7,15,22, 24 – – – – – – 5 2940 y 5,6,7,16,19 4310 y 4,7,10,22,30 5390 y 1,7,9,15,22 – – – – – – 4 2370 y 6,9,19,22 3500 y 6,22,29,30 4520 y 7,10,15,22 5470 n 10,15,22,24 5470 n 10,15,22,24 3 1830 y 4,19,22 2670 y 2,19,29 3400 y 9,10,17 5390 n 6,15,22 5400 n 6,22,24 2 1220 y 9,19 1780 y 4,7 2300 y 6,29 5210 n 9,19 5320 n 9,22 – – α = 90 percent τ = 48 S Obj L Locns τ = 49 Obj L Locns τ = 50 τ = 60 Obj L Locns Obj L MARIANOV & SERRA: LOCATION-ALLOCATION MODELS © Blackwell Publishers 1998. TABLE 3: CPLEX Results for One-Server per Center, Constrained Waiting Time τ = 70 Locns Obj L Locns 3580 n 5,6,8,10,14, 23,24,26,30 5470 n 1,2,3,6,7,18, 22,24,30 – – – – – – – – – 8 3500 y 4,6,9,10,14,1 5,23,30 5390 y 3,7,9,11,13, 15,22,30 – – – – – – – – – 7 2960 y 4,6,7,10,12,1 8,27 4880 y 1,4,6,7,8,9, 22 5470 n 1,10,11,19,2 2,24,29 – – – – – – 6 2710 y 5,6,7,17,18,2 2 4240 y 1,7,15,19, 22,30 5390 y 6,7,9,15,22, 30 – – – – – – 415 9 416 © Blackwell Publishers 1998. α = 90 percent (continued) τ = 48 S Obj L Locns τ = 49 Obj L τ = 50 Locns L τ = 60 Locns 4580 y 3,6,7,22,29 Obj L – – τ = 70 Locns – Obj L – – Locns 5 2300 y 3,6,14,15,17 4 1900 y 9,10,14,17 2750 y 8,17,19,25 3720 y 10,11,15,19 5470 n 9,22,24,29 5470 n 15,20,22,24 3 1430 y 15,19,22 2140 y 9,17,25 2810 y 1,5,25 5390 n 6,15,22 5400 n 9,22,24 1430 y 1,7 1880 y 9,22 5210 n 9,19 5320 n 9,22 2 3530 y 9,10,17,19,22 Obj – α = 95 percent τ = 62 S L τ = 63 τ = 64 Locns Obj L Locns Obj L 11 3580 n 1,2,3,9,10,12, 13,16,24,25,30 – – – – – 10 3520 y 1,2,3,9,10, 12,13,16,24,25 – – – – – 9 3340 y 3,8,9,10,11,12, 13,25,27 4140 n 2,5,13,15,17,2 2,25,29 8 2970 y 3,5,10,12,13, 14,20,25 7 τ = 74 Locns Obj L – – – – – – 5470 n 1,2,6,9,13,1 5,22,24,29 – – 4060 y 2,6,7,8,14, 15,17,29 5390 y 1,2,6,9,13, 22,25,29 – 2670 y 5,6,7,10,12, 14,25 3940 y 1,3,6,7,15, 22,25 5180 y 6,7,8,9,17, 22,30 6 2280 y 2,4,9,12,14,25 3300 y 4,7,10,15, 19,22 5 1900 y 14,15,16,23,25 2730 y 10,12,15, 22,25 4 1570 y 9,13,14,19 τ = 84 Locns Obj L Locns – – – – – – – – – – – – – – – 4350 y 2,7,9,10,15, 19 – – – – – – 3680 y 4,9,10,19,29 – – – – – – 2270 y 7,10,22,30 2890 y 9.10.29,30 5470 n 7,10,22,24 5470 n 6,15,22,24 3 1720 y 22,25,30 2260 y 5,9,26 5390 n 7,10,22 5400 n 6,22,24 2 1160 y 7,16 1520 y 20,23 5210 y 2,6 5320 n 6,22 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 Obj MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 417 Given its capacity-like form when the probabilistic constraint was binding, the run time and the number of branches was much higher than in the case in which the distance limit was binding. Allocation to nearest server was not required, so in many cases demand was not allocated to the nearest server. In reality, users often go to a server which gives better service in terms of waiting time instead of choosing a server because of its closeness. Figure 4 shows an example of this. If allocation to nearest server is desired, additional constraints must be included in the models. Branch-and-Bound; Constrained Queue Length; One or More Servers per Center The queue lengths were set at zero, one, and two users. To show a more congested case the call rates for one server per center were set at 0.015 times the node population and for three and five servers per center the call rates were set at 0.042 times the node population. These latter cases represent a highly congested system. The results for one server are shown in Table 4 for three values of α: 95 percent, 90 percent, and 85 percent. The results for three and five servers TABLE 4: CPLEX Results for One Server per Center, Cconstrained Waiting Time α = 95 percent S 7 6 5 4 3 2 b=0 b=1 b=2 Obj L Locns Obj L Locns Obj L Locns 5470 5390 5260 4170 3170 2140 n y y y y y 1,2,5,6,15,22,24 2,10,11,20,22,23 1,6,7,22,30 7,15,17,25 9,10,25 6,19 – – 5470 5390 5270 3520 – – n y y y – – 2,9,15,22,24 1,6,15,22 9,15,17 10,18 – – – 5470 5390 4540 – – – n n y – – – 6,15,22,24 7,9,22 8,19 α = 90 percent S 5 4 3 2 b=0 b=1 b=2 Obj L Locns Obj L Locns Obj L Locns 5470 5390 4480 3020 n n y y 3,7,10,22,24 6,7,22,25 7,10,22 10,17 – 5470 5390 4440 – n n y – 6,15,22,24 10,19,22 2,7 – 5470 5390 5210 – n n n – 10,15,22,24 6,15,22 9,19 α = 85 percent b=0 S 4 3 2 b=1 b=2 Obj L Locns Obj L Locns Obj L Locns 5470 5390 3700 n n y 10,22,24,29 6,7,22 1,19 5470 5390 5100 n n n 6,15,22,24 9,15,22 6,19 5470 5390 5210 n n n 10,15,22,24 6,15,22 9,19 © Blackwell Publishers 1998. 418 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 are shown in Table 5, for α = 95 percent. Again, the branching limit was reached in many cases, especially for three and five servers per center, resulting in solutions that are not necessarily optimal. Several simple heuristics were developed for obtaining fast solutions of the models. These heuristics were all based on satisfying the time-limit or queuelength restriction. We describe the best one, which follows: Step 1: Make a list of the candidate nodes j; ordered by decreasing call rate fj, if they are also demand nodes, ordered at random if the candidate nodes are not demand nodes. Call this list Dj. Then for each candidate node compute the right-hand side of Equation (8) or Equation (10), depending on the model being utilized. This right-hand side is the limit value of average calls per time unit (λlim, j) that can be accepted by a center located at this node. If this limit value is exceeded, the constraint on queue length or waiting time limitdoes not hold. TABLE 5: CPLEX Results for m Servers per Center. Constrained Waiting Time Three servers per Center, α = 95 percent S B=0 B=1 Obj L Locns 10 4760 n 9 4760 n 8 4680 y 7 4530 y 1,2,3,5,6,11,14, 15,18,24 3,4,5,6,7,10,13, 22,24 3,5,7,9,15,19,22,2 3 1,6,7,9,19,22,29 6 3860 y 5 4 3 2 3230 2490 1960 1320 y y y y b=2 Obj L Locns 5470 n 5390 y 6,9,17,22,25,30 4620 y 6,9,15,16,19 9,10,17,22 4,5,20 15,28 3920 3210 2460 1640 y y y y 2,7,8,9,15, 19,22,24 4,6,8,9,15, 22,30 6,9,10,17, 19,30 6,7,10,13,19 3,9,15,19 15,19,25 19,23 Obj L Locns 5470 n 5390 y 4620 3630 2810 1880 y y y y 10,13,15, 22,24,29,30 6,9,11,22, 29,30 6,7,8,22,29 2,10,15,19 15,19,25 9,22 Five Servers per Center, α = 95 percent b=0 S 6 5 4 3 2 b=1 b=2 Obj L Locns Obj L Locns Obj L Locns 5470 5410 5320 4010 2860 n y y y y 6,10,14,15,24,29 9,10,22,24,29 7,9,22,30 7,11,23 2,15 5470 5390 4570 3080 n n y y 6,10,15,22,24 6,10,15,22 7,8,19 9,22 5470 5390 5070 3400 n y y y 1,9,15,22,24 6,7,15,22 6,15,22 3,13 © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 419 For each candidate node j, make a list of all the demand nodes i within the distance standard, ordered by increasing distance to the candidate node. Call this list Dji. Note that the same demand node can be in several lists Dji. Step 2: For each candidate node, make the current total incoming call rate equal to zero (λinc, j). Step 3: Starting with the first candidate node j on the list Dj, add to its incoming call rate λinc, j, the call rate fi of the first node on the list Dji. Then, add the call rate fi of the second node on the list Dji, then the third, and so on, until the point where adding any extra demand node would exceed the limit value of calls λlim,j. Temporarily allocate all these demand nodes to a hypothetical server at node j. Step 4: Repeat Step 3 for all nodes in list Dj. Note that the same demand node could be temporarily allocated to several candidate locations. Step 5: Locate a server on the node with the highest λinc,j. Take all demand nodes allocated to it out of all the lists Dji of all potential centers. Allocate them definitively to the located server. Step 6: Take each of the servers already located, de-allocate the demands that were allocated to it, and move it to all possible unused candidate locations (repeating each time steps 2 to 5). If some location gives a better objective, keep the server at that location. Repeat for all servers already located. Repeat step 5 until there are no further changes that improve the solution. Step 7: Repeat steps 2 to 6 until all available servers are located. Step 8: Repeat procedure of step 6 for all servers. Two versions of the heuristic were run for each value of the parameters. The first, with step 8, but without step 6; the second, with steps 6 and 8. For all cases, the heuristic took no longer than 2.5 seconds on a 486SX computer running at 100 MHz. In most cases the second version performed better, and we present only the results of the second heuristic. Heuristics; Constrained Waiting Time; One Server per Center Table 6 shows the results obtained by the heuristic. The figures shown correspond to the average percentage difference between the heuristic solution and the CPLEX solution, for different numbers of servers. Also shown are the maximum percentage differences between the two solution methods. It is interesting to note that, for α = 95 percent, τ = 62 m, the solution given by the heuristic is better than the CPLEX solution. This happens in many cases, and is due to the branching limit imposed on CPLEX. Because of this branching © Blackwell Publishers 1998. 420 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 TABLE 6: Results of the Heuristics—Constrained Waiting Time α = 85 percent τ = 40 m 4.17 (7.64)a α = 90 percent τ = 48 m 3.98 (7.43)a α = 95 percent τ = 62 m –0.42 (0.90)a τ = 41 m τ = 42 m τ = 49 m τ = 52 m τ = 62 m 5.04 (9.39) 3.10 (6.12) 0.00 (0.00) 2.06 (13.16) 5.49 (13.16) τ = 49 m τ = 50 m τ = 60 m τ = 70 m 1.21 (3.69) 3.15 (5.34) 1.66 (4.99) 3.16 (4.51) τ = 63 m τ = 64 m τ = 74 m τ = 84 m 1.62 (5.33) 0.78 (5.60) 7.48 (12.86) 2.74 (5.83) aThe figures in parenthesis correspond to maximum values of percentage of difference between heuristic and CPLEX. The figures on top are the average percentages of difference between both solution methods. limit, the solution found by CPLEX is not optimal, and the heuristic finds a better solution than the integer programing package. For one server, constrained waiting time, the run time for both heuristics never exceeded 0.11 seconds. Due to the structure of the heuristics, the solution is robust in many cases, in the sense that when more than k servers are being located, the locations of the first k servers correspond to the best solution to the problem of locating k servers. In other words, the solution (locations) to the problem of locating k servers is a subset of the solutions to the problems of locating more than k servers. However, when CPLEX is used locations may be very different for different numbers of servers, possibly because different CPLEX runs choose different alternate optima. Because of the way the heuristics fill the service capacity of each potential center, in most cases the resulting “districting” consists of attention regions that do not overlap. This differs from the regions found by CPLEX. For example, Figure 5 shows the heuristic solution for τ = 48, α = 90 percent. The CPLEX solution for the same case, displayed in Figure 4, shows a center located on a demand node (17) which is assigned to a different center (14). Also, it shows attention regions that overlap, as in the case of demand nodes 7, 18, and 19. Heuristics; Constrained Queue Length; One or More Servers per Center Table 7 shows the results of the heuristic, compared to the solutions obtained with CPLEX. Column “Ohe” displays the value of the objective and the column labelled “percent” shows the percentage of difference between CPLEX and the heuristic. Table 8 shows the results for three and five servers located at each center, for α = 95 percent. © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 421 FIGURE 4. CPLEX solution for τ = 48, α = 90%. The 4 servers are located at nodes 9,10,14 and 17. The arrows show the allocations. 14 27 12 16 28 22 17 23 30 19 29 3 26 18 7 8 2 4 9 15 1 13 5 11 6 24 10 25 20 21 FIGURE 5. Heuristic solution for τ = 48, α = 90%. The 4 servers are located at nodes 16,18,20 and 23. The demand at these is self-assigned. The arrows show the allocations. © Blackwell Publishers 1998. 422 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 TABLE 7: Results, Constrained Queue Length, One Server per Center α = 95 percent S b=0 b=1 b=2 % Ohe Locns % Ohe Locns % 7 0 5470 6,11,18,1, 22,20,24 0 5470 2,4,17,20,24,2 6,27 – – – 6 0 5390 6,11,18,1, 22,20 1.10 5410 2,4,17,20,24,2 6 – – – 5 0.95 5210 6,11,18,1, 22 2.38 5340 2,4,17,20,24 4 –1.68 4240 1,4,18,19 2.41 5260 2,4,17,20 1.28 5400 2,9,22,21 3 –0.31 3180 1,15,19 3.60 5080 2,4,17 2.97 5230 2,9,22, 2 0.94 2120 1, 6 0.28 3510 2,4 1.54 4470 2,9 % Locns % 0 Ohe 5470 Locns 2,9,22,21,24 α = 90 percent S b=0 b=1 % Ohe Locns 6 0 5470 5 1.46 5390 7,18,1,22, 10,24 7,18,1,22, 10 4 6.49 5040 3 2 2.68 0.67 4360 3000 – Ohe – b=2 – Ohe Locns – – – – – – 0 5470 15,3,22,25,24 7,18,1,22 1.46 5390 15,3,22,25 0 5470 15,30,12,20 7,18,1 7,18 6.12 1.35 5060 4380 15,3,22 15,3 0 0 5390 5210 15,30,12 15,30 α = 85% percent b=0 S b=1 b=2 % Ohe Locns % Ohe Locns % Ohe Locns 5 0 5470 7,10,19,14, 24 0 5470 2,7,22,20,24 0 5470 4,19,10,14,24 4 3 1.46 3.34 5390 5210 7,10,19,14 7,10,19 1.46 3.34 5390 5210 2,7,22,20 2,7,22 1.46 3.34 5390 5210 4,19,10,14 4,19,10 2 0 3700 7,10 9.21 4630 2,7 8.25 4780 4,19 6. CONCLUSIONS New models for locating service centers in a congested situation have been presented. These models, being linear, explicitly include a constraint on service quality, specifically the waiting time or queue length at each center. Heuristic solutions are presented for the models and compared to solutions obtained with a commercial integer programming package (CPLEX). When the run time (or branching) is limited, heuristic solutions appear to be very close to the solutions obtained by using CPLEX. In some cases, heuristic solutions are better. These heuristics are very simple, and may be improved by the introduction of randomselection steps or random restarts. © Blackwell Publishers 1998. MARIANOV & SERRA: LOCATION-ALLOCATION MODELS 423 TABLE 8: Results, Constrained Queue Length, Three and Five Servers per Center Three Servers per Center, α = 95 percent b=0 S b=1 % Ohe Locns 10 0 4760 11,18,20,17, 2,3,4,8,14,24 9 1.68 4680 8 3.85 4500 7 4.86 4310 6 2.85 5 4 3 2 b=2 % Ohe Locns % Ohe 11,18,20,17, 2,3,4,8,14 0 5470 1,2,5,12,152 0,24,27,29 0 5470 11,18,20,17, 2,3,4,8 11,18,20,17, 2,3,4 1.0 5410 1.10 5410 1.1 5330 1,2,5,12,15, 20,24,29 3,15,20,2,4,1, 22 2.56 5330 10,13,15,22, 24,29,30 3750 11,18,20,17, 2,3 –1.95 4710 4.45 5150 6,9,11,22, 29,30 1.24 3190 11,18,20,17,2 –2.04 4000 3,15,20,2,4 3.68 4450 6,7,8,22,29 –3.21 2570 11,18,20,17 –1.56 3260 3,15,20,2 0.83 3600 2,10,15,19 1.53 2.27 1930 1290 11,18,20 11,18 0.41 1.2 3,15,20 20,15 2.85 1.60 2730 1850 15,19,25 9,22 Locns % Ohe Locns 2440 1620 3,15,20,2,4,1 Locns Five Servers per Center, α = 95 percent b=0 S b=1 b=2 % Ohe Locns % Ohe 6 0 5470 8,30,5,25, 12,24 0 5470 5 0.36 5390 1,8,9,19,22 1.46 5390 10,23,2,15,16 0 5470 1,9,15,22,24 4 6.01 5000 8,30,5,25 1.67 5300 10,23,2,15 0 5390 6,7,15,22 3 0.25 4000 8,10,18 –0.88 4610 10,23,15 2.37 4950 6,15,22 2 7.30 2650 19,25 0 10,23 0.29 3390 3,13 3080 REFERENCES Batta, Rajan. 1988. “Single Server Queueing–Location Models with Rejection,” Transportation Science, 22, 209–216. ———. 1989. “A Queueing–Location Model with Expected Service Time–dependent Queueing Disciplines,” European Journal of Operational Research, 39, 192–205. Batta, Rajan, Richard Larson, and Amedeo Odoni. 1988. “A Single-Server Priority Queueing– Location Model,” Networks, 8, 87–103. Berman, Oded and Richard Larson. 1985. “Optimal 2–Facility Network Districting in the Presence of Queueing,” Transportation Science, 19, 261–277. Berman, Oded, Richard Larson, and Samuel Chiu. 1985. “Optimal Server Location on a Network Operating as a M/G/1 Queue,” Operations Research, 12, 746–771. Berman, Oded, Richard Larson, and Celik Parkan. 1987. “The Stochastic Queue p-Median Location Problem,” Transportation Science, 21, 207–216. Berman, Oded and R. Mandowsky. 1986. “Location–Allocation on Congested Networks,” European Journal of Operational Research, 26, 238–250. © Blackwell Publishers 1998. 424 JOURNAL OF REGIONAL SCIENCE, VOL. 38, NO. 3, 1998 Brandeau, Margaret, Samuel Chiu, Shiv Khumar, and Thomas Grossman Jr. 1995. “Location with Market Externalities,” in Zvi Drezner (ed.), Facility Location : A Survey of Applications and Methods. New York: Springer-Verlag. Church, Richard and Charles ReVelle. 1974. “The Maximal Covering Location Problem,” Papers of the Regional Science Association, 32, 101–118. Hakimi, S. Louis. 1964. “Optimal Locations of Switching Centers and the Absolute Centers and Medians of a Graph,” Operations Research, 12, 450–459. Larson, Richard and Amedeo Odoni. 1981. Urban Operations Research. Englewood Cliffs, New Jersey: Prentice–Hall. ReVelle, Charles, and Ralph Swain. 1970. “Central Facilities Location,” Geographical Analysis, 2, 30–42. Wolff, Ronald. 1989. Stochastic Modeling and the Theory of Queues. Englewood Cliffs, New Jersey: Prentice–Hall. © Blackwell Publishers 1998.