- Email: [email protected]
The gravity multiple server location problem
Computers & Operations Research 38 (2011) 694–701
Contents lists available at ScienceDirect
Computers & Operations Research journal homepage: www.elsevier.com/locate/caor
The gravity multiple server location problem Tammy Drezner, Zvi Drezner Steven G. Mihaylo College of Business and Economics, California State University-Fullerton, Fullerton, CA 92834, United States
a r t i c l e in f o
a b s t r a c t
Available online 22 August 2010
Models for locating facilities and service providers to serve a set of demand points are proposed. The number of facilities is unknown, however, there is a given number of servers to be distributed among the facilities. Each facility acts as an M/M/k queuing system. The objective function is the minimization of the combined travel time and the waiting time at the facility for all customers. The distribution of demand among the facilities is governed by the gravity rule. Two models are proposed: a stationary one and an interactive one. In the stationary model it is assumed that customers do not consider the waiting time at the facility in their facility selection decision. In the interactive model we assume that customers know the expected waiting time at the facility and consider it in their facility selection decision. The interactive model is more complicated because the allocation of the demand among the facilities depends on the demand itself. The models are analyzed and three heuristic solution algorithms are proposed. The algorithms were tested on a set of problems with up to 1000 demand points and 20 servers. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Facility location Heuristics Multiple facilities
1. Introduction We consider models for locating facilities and allocating p servers among them so as to minimize customers’ travel time plus waiting time at the facilities. The number of facilities is a variable limited to a maximum of p facilities. We allow the location of one or more servers at each facility. Wang et al. [43] considered the problem of locating facilities that are modeled as M/M/1 queuing systems. Customers travel to the closest facility and a constraint is placed on the maximum expected waiting times at all facilities. The objective is to minimize the combined travel and waiting times. Berman et al. [5] considered a similar model with additional constraints on lost demand due to congestion and insufficient coverage. The objective is to minimize the number of facilities. Berman and Drezner [4] extended the model suggested by Wang et al. [43] by allowing more than one server to be located at each facility thus employing the M/M/k queuing system. Berman and Drezner [4] assume that customers patronize the closest facility. Aboolian et al. [1] investigated the Berman and Drezner [4] model applying the minimax objective rather than minimum. Typical to all of these models is the assumption that service is provided by the closest facility. This assumption is used in many location-allocation models such as the p-median, p-center or competitive models [27,28,24,36,6,31,18,29,40]. The proximity assumption is logical because it is in the user’s best interest to use Corresponding author. Tel.: + 1 714 657 278 2712.
E-mail address: [email protected] (Z. Drezner). 0305-0548/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2010.08.006
the services of the closest facility. This is indeed the best choice if the demand is allocated to facilities by a planner rather than by customers’ free choice and/or when facilities are equally attractive. Under the proximity assumption, customers use only one facility to the exclusion of all others, resulting in an ‘‘all or nothing’’ behavior. In reality there are many reasons why the proximity assumption, which leads to the ‘‘all or nothing’’ customer behavior, does not accurately describe customers’ behavior: (i) A small change in the location of facilities may discontinuously shift the whole demand at a demand point from one facility to another. Suppose that the closest facility is at a distance of 4.9 miles and the second closest one is at a distance of 5 miles. The whole demand is allocated to the first facility. However, if the first facility is relocated to a location at a distance of 5.1 miles from the demand point, the whole demand will now be allocated to the second facility. This does not reflect actual customer behavior. (ii) Different facilities have different levels of attractiveness and customers are willing to travel an extra distance to a farther and more attractive facility. This issue in competitive location models is addressed by Drezner [8,9] and in various utility (additive or multiplicative) models [35,11,16]. (iii) Even if the facilities are equally attractive, it is not clear that all customers at a demand point patronize the same facility. Customers do not necessarily measure the exact distance when deciding which facility to patronize for service. They select the facility that they perceive to be the closest to them.
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
Distances are likely perceived differently by different customers. Therefore, the selected facility to patronize will not be the same one for all customers at a demand point. Distance perception is discussed in Drezner and Eiselt [17]. (iv) Demand points are not mathematical points. Rather, they represent an area inhabited by customers. It is possible that the facility closest to customers residing on one side of the area is not the closest to customers residing on the other side of the same area. This is referred to as Type C error in the aggregation literature [30,23,22]. This issue can be illustrated by a Voronoi diagram based on the facilities [37,42]. A Voronoi diagram is a tessellation of the plane into polygons. Each polygon consists of all the points closest to one facility which must be located inside that polygon. The Voronoi diagram may partition a demand area into several parts, each closest to a different facility. Incorporating area of demand points in competitive location gravity-based models is addressed in Drezner and Drezner [12]. It is therefore more accurate and more realistic to employ the gravity rule proposed by Reilly [39] and later used by Huff [32,33] to explain customers’ spatial choice behavior and estimating buying power attracted to competing facilities. According to the gravity rule, the probability that a customer patronizes a facility is proportional to the facility’s attractiveness, and to a decreasing function of the distance to the facility, commonly known as the distance decay function. Huff [32,33] suggested a power of the distance as a decay function. Others suggested an exponent of the distance [44]. Drezner [10] showed empirically that an exponent of the distance more accurately describes customers’ behavior than a power of the distance. The gravity rule rectifies the four concerns raised above about the proximity rule. Some recent papers have re-analyzed well known models using the gravity rule rather than the proximity rule. For example, Drezner and Drezner [13] formulated and solved the hub location problem using the gravity rule. Drezner and Drezner [15] analyzed and solved p-median problems on a network and Drezner and Drezner [14] solved the p-median in the plane when allocation of customers among the facilities is governed by the gravity rule. The contribution of the current paper is two fold: The gravity rule is used for modeling, and two models are proposed, analyzed, and solved. The stationary model] following the one suggested by Berman and Drezner [4] but assuming that customers patronize a facility according to the gravity rule. The interactive model assuming that customers select the facility by incorporating the expected waiting time at the facility in the selection criteria. Customers may experience a long waiting time at a facility in a certain visit, thus selecting another facility in a following trip. The selected facility may be farther away but have shorter expected waiting time. This model is more difficult to analyze and solve because the attractiveness of the facility depends on the demand. Applications for such models include the location of automatic teller machines (ATMs) where at each facility one or more ATMs are installed; location of bank branches where a decision is required about the number of teller stations to be built at each branch; location of gas stations with the number of pumps in each gas station being a variable; location of post office branches, department of motor vehicle offices, or similar service establishments where the number of windows (servers) is also a variable. The paper is organized as follows: In the next section we formulate, analyze, and propose solution methods for the stationary model. In Section 3 we define and analyze the interactive model. In Section 4 we report computational experiments on a
695
set of test problems for both the stationary and the interactive models. We conclude in Section 5 and propose ideas for future research.
2. The stationary model The analysis of the stationary model follows the approach in Berman and Drezner [4]. We briefly summarize the changes required for implementing the approach in Berman and Drezner [4]. 2.1. Notation Let N wi
be a set of demand points of cardinality n be the mean demand generated at demand point iA N in a Poisson process M be a set of potential locations for facilities of cardinality m be the attractiveness of a facility located at potential Aj location j A M p be the total number of servers available for installation at the facilities, one or more servers at each facility S D M of unknown cardinality s rp be an unknown set where facilities are located be the distance between demand point i A N and dij potential location j A M measured in time units by any distance measure f(d) be the decay function of the distance governing the gravity rule Wk ðl, mÞ be the expected time spent in line at an M/M/k system consisting of k servers with arrival rate l and service rate
m Lk ðl, mÞ kj Z 1 K
¼ lWk ðl, mÞ be the average length of the line in an M/M/ k system of k servers be the unknown number of servers allocated to facility jAS be the unknown vector K ¼ fk1 , . . . ,ks g of the allocation of servers to the set of facilities S
2.2. Formulation The average demand lj ðSÞ for j A S is by the gravity rule:
lj ðSÞ ¼
X
wi P
iAN
Aj f ðdij Þ
ð1Þ
j A S Aj f ðdij Þ
where f ðdÞ ¼ ecd (or any other distance decay function). The objective function, to be minimized, is the combined travel time and expected time spent at the selected facility for all customers. Note that the total travel time TðSÞ to the facilities is independent of the servers allocation: TðSÞ ¼
X iAN
wi
X Aj f ðdij Þ dij P j A S Aj f ðdij Þ jAS
ð2Þ
The total waiting time at the facilities, WðS,KÞ, is X X WðS,KÞ ¼ lj ðSÞWkj lj ðSÞ, m ¼ Lkj lj ðSÞ, m jAS
ð3Þ
jAS
The objective function, to be minimized, is FðS,KÞ ¼ TðSÞ þWðS,KÞ
ð4Þ
696
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
2.3. Analysis of the stationary model
3. The interactive model
We first devise an algorithm that, for a given set S DM, finds the optimal allocation of p servers among the members of S. Each member in S is allocated at least one server.
Suppose that customers have information about the expected waiting time at each facility. Customers will base their selection on the combined travel and waiting time. This entails replacing dij in (1) with dij þ Wkj ðlj ðSÞ, mÞ when calculating the allocation of customers to facilities. If the decay function f ðdÞ ¼ ecd is used, it
2.3.1. Calculating the total time spent at a facility This algorithm is based on Pasternack and Drezner [38] and is also given in Berman and Drezner [4]. Consider an M/M/k facility with mean arrival rate l and mean service rate m. Calculate the sequence far g for r ¼ 1, . . . ,k: a1 ¼ 1;
ar þ 1 ¼ 1 þr
m l
ar
ð5Þ
then Lk ðl, mÞ ¼
l2
ðkmlÞ l þ ðkmlÞak
cWk ðlj ðSÞ, mÞ
means replacing ecdij with fe be interpreted as defining
l ðkmlÞ½l þ ðkmlÞak
2.3.2. Finding the best allocation of servers The best allocation of servers is found by the following algorithm [4]. 1. For each node j A S calculate the minimum of servers number necessary to service the demand kj ¼ int lj ðSÞ=m þ 1. P 2. If sj ¼ 1 kj 4 p, there is no feasible solution. Ps 3. If j ¼ 1 kj ¼ p, the solution is fkj g. P P 4. If sj ¼ 1 kj o p, repeat the following until sj ¼ 1 kj ¼ p. (a) Check the decrease in Lkj (6) by changing kj to kj þ 1 for j ¼1,y,s. (b) Increase the number of servers by one at the facility with the maximum decrease Lkj Lkj þ 1 . Note that changing kj to kj þ 1 is very simple once the vector fakj g is saved. Since akj is known, akj þ 1 is calculated by (5) and then Lkj þ 1 is calculated by (6). The value of akj for the facility with the increased number of servers is updated. Since the best allocation of servers to facilities is optimally found by the approach above, the function FðSÞ ¼ minfFðS,KÞg K
ð7Þ
is well defined and can be calculated for any set S D M.
2.4. Solving the stationary problem The analysis in Section 2.3.2 enables us to calculate the optimal value of the objective function for a given set S D M of cardinality sr p. The optimization problem is therefore reduced to finding the set S that minimizes FðSÞ (7). Berman and Drezner [4] constructed a descent algorithm, a tabu search algorithm, a simulated annealing algorithm, and a genetic algorithm. We constructed the first three algorithms to heuristically solve the stationary problem of finding the best subset S D M. The genetic algorithm developed by Berman and Drezner [4] required about five times the computer time than the other metaheuristics and therefore was not tested. The heuristic algorithms are detailed in the Appendix.
ð8Þ
j
where Aj(S) is the attractiveness of facility j as is applied in competitive facilities location models ([8,9,2]) but it is not fixed. Rather, it is a function of the allocation of demand which is determined by the set S. Equation (1) is then replaced by
lj ðSÞ ¼ ð6Þ
gecdij . This formula can
cWk ðlj ðSÞ, mÞ
Aj ðSÞ ¼ e
X iAN
Wk ðl, mÞ ¼
j
wi P
Aj ðSÞf ðdij Þ
ð9Þ
j A S Aj ðSÞf ðdij Þ
where Aj(S) is defined by (8). We note that using the number of servers at a facility as the attractiveness measure is not appropriate. It is possible to have longer lines at a facility with more servers because there may be higher demand at such a facility. The difficulty with solving such a model is that the ‘‘attractiveness’’ Aj(S) is a function of the allocation of demand and is not fixed. The same issue has risen in other papers that analyzed models where the allocation of demand depends on the demand itself [19–21]. Consider the problem of minimizing F(S,K) for a given S. If we are able to find minK fFðS,KÞg, then the heuristic approaches designed to find the best set S for solving the stationary problem can be employed for solving the interactive problem. This problem for a given set S DM can be represented as a traffic flow problem [7,41]. The following network is designed (see Fig. 1). One set of nodes is N and a second set of nodes is S. An additional sink node is also defined. Traffic (demand) wi is generated at each iA N and the destination of each traveler is at the sink. A link connects every demand point iA N with a facility j A S. The length of the link is dij and does not depend on the flow on the link. Each facility j A S is connected to the sink. The length of such a link is the waiting time at the facility, which depends on the traffic on the link (congestion). A customer selects a facility, travels to it, and then uses the link connecting the facility to the sink which is the waiting time at the facility. The waiting time depends on the volume of traffic Demand
dij
Facilities
Wk Sink
Fig. 1. The traffic flow representation.
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
(total demand) on the link connecting the facility and the sink. Such a traffic flow problem has an equilibrium [7,41].
3.1. Solving the traffic flow problem
697
general, using small values of y will always lead to convergence, but may not be the most efficient.
3.2. Solving the interactive model
We follow the widely used methods for solving traffic flow problems [41]. Consider a given set S D M. For a given attractiveness levels Aj(S), for example Aj ðSÞ ¼ 1 8j, the values of lj ðSÞ are calculated. The best solution K can be found as detailed in Section 2.3.2. The values of Lk and Wk change. Therefore, the attractiveness levels Aj(S) change, and the problem is re-solved, and so on. As is experienced in traffic flow problems, such an iterative algorithm may not converge to an equilibrium. We therefore dampen the oscillations in the values of Aj ðSÞ by changing Aj(S) by a proportion 0 o y r 1, for example y ¼ 0:1, of the change to the new one. This means that if Auj ðSÞ are the values calculated by (8), Aj(S) is replaced by Aj ðSÞ þ yðAuj ðSÞAj ðSÞÞ. This can be interpreted as assuming that only 10% of the customers are cognizant of the change in attractiveness and modify their behavior. In the network flow algorithm an incremental adjustment of the traffic flow rather than full adjustment leads to an equilibrium. In our experiments using y ¼ 0:1 always converged to an equilibrium. In
We propose to apply the same metaheuristic approaches detailed in the Appendix by applying the solution approach described in Section 3.1 for finding K to define F(S). Note that the solution approach in Section 3.1 finds an equilibrium solution which may be a local minimum.
3.3. An extension to the interactive model Suppose that only a proportion 0 r y r 1 of the customers are cognizant of the waiting time at the facilities and consider it in their facility selection. This entails replacing dij in (1) with dij þ yWkj ðlj ðSÞ, mÞ when calculating the allocation of customers to facilities. The attractiveness is thus Aj ðSÞ ¼ e
ycWk ðlj ðSÞ, mÞ j
Table 1 Optimal solutions for the stationary model for a given number of facilities. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 a
Not run.
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
% over best known (Table 2)
Run rime (min)
s ¼1
s ¼2
s¼3
s¼4
s ¼5
s ¼1
s¼ 2
s¼3
s¼4
s¼5
64.71 69.41 76.29 75.35 79.49 92.25 100.07 99.83 94.96 119.70 118.18 122.34 101.07 133.95 134.15 135.93 102.24 146.60 151.24 146.60 101.72 152.27 161.65 156.20 103.54 158.91 170.56 169.43 103.41 161.62 178.62 181.43 106.39 167.35 189.83 190.14 105.43 168.44 194.11 200.22
27.93 28.37 31.88 30.73 38.18 47.68 53.24 53.11 50.21 68.02 66.86 69.57 54.43 79.02 78.21 79.52 53.54 85.16 88.14 84.41 55.16 90.26 97.05 92.94 58.61 97.43 106.22 105.23 58.98 101.53 114.62 116.60 61.27 105.87 122.79 123.28 60.89 107.08 126.42 131.47
6.91 9.14 10.42 2.73 16.80 19.22 20.57 20.11 24.48 33.31 30.53 32.60 24.91 44.28 42.32 42.39 24.99 51.61 51.75 47.76 27.86 55.84 61.78 57.86 30.53 59.36 66.67 65.47 30.44 60.90 71.27 72.74 32.98 64.14 78.04 77.82 32.48 64.56 80.87 83.76
0.00 1.68 1.48 0.40 5.71 9.30 7.37 6.36 12.89 19.22 13.36 12.86 14.35 21.30 19.96 20.37 14.16 25.90 27.54 24.17
1.28 0.35 0.00 0.00 0.00 4.89 4.25 3.25 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.03 0.03 0.03 0.03 0.05 0.05 0.05 0.05 0.07 0.07 0.07 0.07 0.10 0.10 0.10 0.10 0.14 0.14 0.14 0.14
0.01 0.01 0.01 0.01 0.09 0.09 0.09 0.09 0.46 0.47 0.47 0.48 1.51 1.52 1.53 1.54 3.84 3.85 3.88 3.90 7.85 7.88 7.92 7.96 14.56 14.62 14.69 14.77 24.84 24.83 24.99 25.08 39.64 39.71 39.87 40.00 60.05 60.18 60.38 60.56
0.17 0.17 0.17 0.18 5.43 5.43 5.48 5.56 41.71 41.74 42.00 42.42 181.55 183.07 183.94 183.92 565.57 566.83 569.87 575.04
3.93 3.94 3.96 4.03 257.30 258.41 260.54 260.42 2987.53
a a a a a a a a a a a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a a a a a a a a a a a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
698
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
When y ¼ 0 the model is stationary. Otherwise, the methods proposed for the solution of the interactive model can be used for the solution of this problem with y 4 0.
demand resulting in problems which are easier to solve. The exponential decay function used is f ðdÞ ¼ e4d . 4.1. The stationary model results We attempted to require about the same total run times for the three algorithms tested. The descent algorithm was replicated 2000 times. Following extensive experiments the following parameters were selected. For the simulated annealing algorithm: the initial temperature T0 was set to 100; the number of iterations was set to N¼1000 np, and the temperature reduction factor was set to a ¼ 15=N. For the tabu search: the tabu tenure was randomly selected each iteration in the range ½0:1n,0:2n, and the number of iterations was set to 100 times the number of iterations in the descent phase. Each problem was solved ten times by both the simulated annealing algorithm and the tabu search. In some situations the planner wishes to have a given number of facilities. The number of servers is assumed to be p as in the original problem. By total enumeration of all possible sets S D M of cardinality 1, 2, or 3 the optimum for one, two, or three facilities was found. We also found the optimal solution for four and five
4. Computational experiments Programs were coded in Fortran using double precision arithmetic and compiled by the Compaq 6.6 Fortran compiler. The programs were run on a desk top Core 2 Duo PC, 2.3 GHz with 2 MB RAM. We randomly generated problems with n ¼ 100,200, . . . ,1000 points in a square of side 10. These points are both generating demand and potential locations for the facilities. Four values of p ¼5, 10, 15, 20 were tested for each value of n. Distances are assumed Euclidean. Each point generates demand wi A ð0,1Þ. P When a service rate of m ¼ ð1=pÞ ni¼ 1 wi is used, total demand is equal to the total service rate. The service rate at each server P was therefore calculated as m ¼ ð1:1=pÞ ni¼ 1 wi so it is 10% higher than total demand. This value of m defines difficult problems because there are many infeasible solutions. If one applies, for P example, m ¼ ð2=pÞ ni¼ 1 wi , capacity of the system is double the Table 2 Best solutions found for the stationary model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Best known
132.765 128.564 123.474 123.841 213.656 199.022 191.182 191.152 298.528 264.516 266.294 261.082 382.856 328.668 328.338 325.646 476.711 390.604 383.343 390.343 567.167 453.174 436.873 445.964 651.937 512.174 490.068 491.941 735.743 571.720 536.787 531.237 822.098 634.331 585.084 584.289 912.215 697.788 636.839 623.689
Descent
Sim. anneal.
Tabu search
a
b
a
b
a
b
23 0 49 2 0 0 7 0 6 0 0 0 3 0 0 0 15 0 0 0 8 0 0 0 61 0 0 0 22 0 0 0 100 0 0 0 30 0 0 0
0.00 0.35 0.00 0.00 0.18 2.67 0.00 1.25 0.00 2.98 0.63 1.61 0.00 3.20 1.75 2.78 0.00 6.50 6.44 2.00 0.00 1.29 9.35 3.08 0.00 1.77 10.95 3.47 0.00 0.96 12.98 6.98 0.00 1.37 13.36 7.39 0.00 0.58 15.25 10.06
9 4 5 2 3 2 1 4 10 1 0 1 7 1 1 0 10 3 1 1 4 2 1 0 10 3 1 1 8 3 1 1 10 3 1 1 10 2 1 1
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.12 0.00 0.00 0.00 0.00 0.84 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.77 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
10 0 3 2 2 2 7 0 10 0 1 0 7 0 0 1 7 0 0 0 9 1 0 1 10 0 0 0 9 1 0 0 10 0 0 0 10 1 0 0
0.00 0.35 0.00 0.00 0.00 0.00 0.00 0.15 0.00 2.06 0.00 1.03 0.00 0.61 0.45 0.00 0.00 0.30 2.94 0.68 0.00 0.00 2.83 0.00 0.00 0.24 2.68 2.45 0.00 0.00 2.00 5.00 0.00 0.59 1.59 4.04 0.00 0.00 0.92 2.69
8.15
Number of times best known solution found. Percentage of best found solution above best known solution.
3.28
3.25
0.09
2.60
0.84
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
facilities for some of the smaller problems. The results are summarized in Table 1. The best known solution for these problems is the best result obtained by the three heuristics reported in Table 2. We report the percentage of the optimal solution for a given number of facilities, s, above the best known solution for the particular problem, and the run time in minutes required for such a solution. If only one facility is allowed, the optimal value of the objective function can be more than three times the optimal solution (more than 200% above the best known solution). Note that for the problem with n ¼100 and p ¼5, s¼4 is optimal and for the problems with n¼ 200, 300 and p ¼5, s¼5 is optimal. The optimality of the best known solution for these three problems was thus proven. The performance of the descent, simulated annealing, and tabu search algorithms is reported in Tables 2 and 3. In Table 2 we report the best known solution, and the number of times (out of 2000 replications for the descent algorithm and 10 replications for the other two algorithms) that the best known solution was found, and the percentage of the best solution found by the corresponding algorithm above the best known solution. In Table 3 we report the percentage of the average solution for all
Table 3 Average solutions found for the stationary model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Sim. anneal.
Tabu search
a
Descent b
a
b
a
b
12.55 9.84 7.15 3.40 17.56 14.49 9.62 7.33 19.06 17.01 8.66 7.82 16.05 18.81 11.20 9.44 10.63 20.79 15.55 9.80 7.33 22.05 18.52 11.70 5.50 22.60 20.23 14.73 3.26 21.48 21.53 17.78 3.02 20.37 23.36 17.99 2.31 18.08 24.02 20.81
0.20 0.22 0.26 0.24 0.98 1.37 1.83 1.92 3.15 5.42 6.04 6.82 7.04 11.94 14.39 15.08 14.76 23.50 26.72 31.26 24.84 40.41 47.23 58.27 37.78 64.14 75.65 96.28 51.68 102.73 117.08 146.87 68.41 162.38 169.18 223.48 86.95 247.59 235.50 313.88
0.05 0.87 0.27 0.06 0.25 1.62 2.09 0.87 0.00 2.05 4.61 4.65 0.06 0.96 2.56 3.85 0.00 0.50 2.23 3.73 0.04 0.37 1.48 3.10 0.00 0.21 1.80 2.21 0.01 0.34 1.87 1.73 0.00 1.00 0.64 1.18 0.00 0.18 0.61 1.00
0.28 0.70 1.15 1.71 1.12 3.16 5.64 8.30 2.51 7.68 14.68 22.09 4.42 14.13 28.08 42.45 6.96 22.58 44.92 74.09 9.99 32.36 66.00 110.95 13.66 44.34 91.10 156.69 17.75 57.58 117.75 210.12 22.45 72.84 149.68 268.66 27.69 89.89 185.05 333.46
0.00 1.34 0.75 0.11 0.30 1.84 0.13 0.98 0.00 2.92 1.45 1.66 0.05 2.24 2.31 1.38 0.01 1.40 5.51 1.27 0.00 0.82 5.01 2.62 0.00 1.19 5.82 4.09 0.01 0.72 4.59 6.41 0.00 1.24 3.63 6.54 0.00 0.52 3.59 7.29
0.17 0.23 0.22 0.23 0.86 1.58 1.86 1.79 2.26 6.37 5.49 7.83 4.58 13.80 20.80 21.34 8.61 30.48 36.54 53.24 12.83 48.47 69.06 79.63 24.68 70.92 117.26 134.02 30.25 107.19 171.92 200.97 40.81 160.35 218.61 307.36 47.76 178.84 351.02 451.04
14.09
63.59
1.23
59.62
1.99
76.03
Percentage of average solution above best known solution. Total time in minutes.
699
replications above the best known solution for the corresponding algorithm and the total run time in minutes for all replications. By examining the results reported in Tables 2 and 3 we conclude that the simulated annealing algorithm performed much better than the other two algorithms for this set of problems and is the recommended algorithm for the stationary model. 4.2. The interactive model results We performed the same experiments for the interactive model. Since run time is longer, we experimented with n r500 problems and the descent algorithm was run only 1000 times. Following extensive experiments we also changed the number of iterations pffiffiffi for the simulated annealing algorithm to N ¼ 100np n. We did not repeat the given number of facilities experiment because the best K found by the iterative procedure is not proven optimal. Note that the s ¼1 results are identical to the results for the stationary model. We report in Tables 4 and 5 the performance measures that were reported in Tables 2 and 3 for the stationary model. By examining the results reported in Tables 4 and 5 we conclude that the tabu search algorithm performed better than the other two algorithms and is the recommended algorithm for the interactive model.
5. Conclusions Two models for locating facilities to serve a set of demand points were proposed. The number of facilities is unknown. However, there is a given number of servers to be distributed among the facilities. Each facility acts as an M/M/k queuing system. The objective function is the minimization of the combined travel time and the waiting time at the facility for all customers. The distribution of demand among the facilities is governed by the gravity rule. Two models were proposed. In the stationary model it is assumed that customers do not consider the expected waiting time at the facility in their selection decision. In Table 4 Best solutions found for the interactive model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Best known
129.618 125.115 121.609 120.376 208.688 194.731 189.612 187.333 296.184 259.033 260.194 253.655 380.856 323.206 319.310 315.332 474.793 385.589 374.061 374.990
Sim. anneal.
Tabu search
a
Descent b
a
b
a
b
49 0 127 39 106 0 220 17 172 4 0 0 27 3 0 0 130 7 0 0
0.00 0.25 0.00 0.00 0.00 0.80 0.00 0.00 0.00 0.00 0.21 0.28 0.00 0.00 1.52 0.54 0.00 0.00 2.23 0.93
10 1 9 8 10 0 1 4 10 3 0 0 10 1 0 0 10 4 0 0
0.00 0.00 0.00 0.00 0.00 0.27 0.00 0.00 0.00 0.00 0.64 0.08 0.00 0.00 0.08 0.24 0.00 0.00 0.21 0.09
10 5 10 10 10 1 7 10 10 8 2 3 10 6 2 2 10 9 1 1
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
45.05
0.34
4.05
0.08
6.35
Number of times best known solution found. Percentage of best found solution above best known solution.
0.00
700
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
Table 5 Average solutions found for the interactive model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Average a b
Descent
Sim. anneal.
Tabu search
a
b
a
b
a
b
5.01 4.48 1.99 2.36 7.04 7.07 2.44 3.27 3.82 7.78 4.59 4.30 2.60 7.89 6.75 5.33 1.88 6.06 10.04 5.29
1.49 3.13 8.28 11.80 6.75 15.87 36.19 55.27 19.06 49.06 79.06 129.13 37.14 102.70 157.37 230.19 59.45 197.68 260.96 393.20
0.00 0.23 0.01 0.01 0.00 0.76 0.06 0.15 0.00 0.12 1.42 0.29 0.00 0.17 0.64 0.91 0.00 0.09 0.49 0.38
1.73 5.93 15.78 31.85 8.32 19.61 44.54 74.66 23.06 53.03 99.55 161.84 47.16 102.30 186.01 294.88 76.30 173.22 301.94 502.47
0.00 0.13 0.00 0.00 0.00 0.48 0.02 0.00 0.00 0.05 0.73 0.16 0.00 0.09 0.41 0.50 0.00 0.03 0.28 0.74
1.66 3.59 10.08 14.04 7.32 17.04 31.33 44.70 20.39 44.96 52.63 84.62 38.15 107.49 121.23 173.69 66.23 228.54 244.87 301.53
5.00
92.69
0.29
111.21
0.18
80.71
Percentage of average solution above best known solution. Total time in minutes.
the interactive model we assume that customers know the expected waiting time at the facility and consider it in their selection decision. The interactive model is more complicated to analyze and solve because the allocation of the demand among the facilities depends on the demand itself. Three heuristic procedures for the solution of the two models were proposed: a descent algorithm, simulated annealing, and tabu search. The recommended approaches are simulated annealing for the stationary model and the tabu search for the interactive model. As future research we propose to investigate a host of related models such as an expansion of an existing system by opening new facilities and/or expanding the number of servers in existing facilities due to increase in demand, competitive models based on the framework proposed in this paper, and many variants of competitive location models as described in Berman et al. [3]. Note that the traffic flow algorithm can be employed for the interactive model based on the proximity assumption leading to another extension of the model in Berman and Drezner [4].
Appendix A The three heuristics are detailed in this appendix. A.1. The descent algorithm The descent algorithm is quite straightforward. A neighborhood of a set of sites S of cardinality s is defined by either removing a site from S (if s 4 1), or adding a site to S (if s o p), or adding a site to S and removing another one from S. The cardinality of such a neighborhood for 1 o s op is n þ sðnsÞ. 1. Start with a randomly generated S and calculate F(S). In our experiments we started with s¼2 because F(S) is feasible for any S of cardinality s¼1 and thus there is always a feasible solution in the neighborhood of the starting S.
2. Check all the subsets S0 in the neighborhood of S and calculate F(S0 ) for each. Find the minimum of these values Fmin and the subset Smin. 3. If Fmin oFðSÞ change S to Smin and go to Step 2. 4. Otherwise, the descent algorithm terminates with a solution S.
A.2. Simulated annealing The simulated annealing algorithm [34] simulates the cooling process of hot melted metals. The variant used in this paper depends on three parameters: the starting temperature T0, the number of iterations N, and the factor a o 1 by which the temperature is reduced every iteration. The neighborhood of S is defined in the descent algorithm. 1. Start with a randomly generated S and calculate F(S). Set the temperature T to T0. 2. A move in the neighborhood of S is randomly selected. Each move is selected with the same probability. The move defines S0 . 3. Calculate F(S0 ). 4. If FðSuÞ r FðSÞ change S to Su and go to Step 6. 5. Otherwise, (a) calculate d ¼ ðFðSuÞFðSÞÞ=T. (b) Change S to Su with probability ed . Otherwise, retain S. 6. Set T ¼ aT and go to Step 2 unless the predetermined number of iterations is reached. 7. The best encountered solution throughout the process is the resulting solution.
A.3. Tabu search A tabu search [25,26] proceeds from the point that the descent algorithm terminated. We establish a tabu list of ‘‘reverse’’ moves (adding back to S a site removed from S), and moves remain in the tabu list for a given tenure. Note that we just need to keep the last iteration number at which a site was removed from S for all the sites in M. 1. Apply the descent algorithm detailed above. Let S be the result of the descent algorithm. F* is set to F(S) and S* is set to S. Empty the tabu list. 2. Check all the subsets S0 in the neighborhood of S and calculate F(S0 ) for each. 3. Find the minimum of these values Fmin and the subset Smin. 4. If Fmin oF , change S and S to Smin , and go to Step 2. 5. Otherwise, (a) Find the minimum FðSuÞ for moves not in the tabu list (i.e., the list of sites that were removed from S during the tabu tenure). The minimal S0 , whether improving or not, is the S for the next iteration. (b) Add the facility removed from S (if any) to the tabu list. (c) Remove the sites whose tenure exceeds the tabu tenure from the tabu list. (d) If the stopping criterion is not met, go to Step 2. Otherwise, stop the search and S* is the resulting solution. We used a random tabu tenure each iteration as detailed in Section 4. To operationalize the tabu list, the last iteration when a site was removed from S is recorded for each j A M (a large negative number is recorded for sites that were never removed form S). A site is in the tabu list if and only if the difference between the current iteration number and its recorded last
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
iteration number is less than or equal to the tabu tenure. This is an easy way to handle varying tabu tenure values for each iteration. References [1] Aboolian R, Berman O, Drezner Z. The multiple server center location problem. Annals of Operations Research 2009;167:337–52. [2] Aboolian R, Berman O, Krass D. Efficient solution approaches for discrete multi-facility competitive interaction model. Annals of Operations Research 2009;167:297–306. [3] Berman O, Drezner T, Drezner Z, Krass D. Modeling competitive facility location problems: new approaches and results. In: Oskoorouchi M, editor. Tutorials in operations research. San Diego: INFORMS; 2009. p. 156–81. [4] Berman O, Drezner Z. The multiple server location problem. Journal of the Operational Research Society 2007;58:91–9. [5] Berman O, Krass D, Wang J. Locating service facilities to reduce lost demand. IIE Transactions 2006;38:933–46. [6] Current J, Daskin M, Schilling D. Discrete network location models. In: Drezner Z, Hamacher HW, editors. Facility location: applications and theory. Berlin: Springer-Verlag; 2002. p. 81–118. [7] Dafermos S. Traffic equilibrium and variational inequalities. Transportation Science 1980;14:42–54. [8] Drezner T. Locating a single new facility among existing unequally attractive facilities. Journal of Regional Science 1994;34:237–52. [9] Drezner T. Competitive facility location in the plane. In: Drezner Z, editor. Facility location: a survey of applications and methods. New York, NY: Springer; 1995. p. 285–300. [10] Drezner T. Derived attractiveness of shopping malls. IMA Journal of Management Mathematics 2006;17:349–58. [11] Drezner T, Drezner Z. Competitive facilities: market share and location with random utility. Journal of Regional Science 1996;36:1–15. [12] Drezner T, Drezner Z. Replacing discrete demand with continuous demand in a competitive facility location problem. Naval Research Logistics 1997;44:81–95. [13] Drezner T, Drezner Z. A note on applying the gravity rule to the airline hub problem. Journal of Regional Science 2001;41:67–73. [14] Drezner T, Drezner Z. Multiple facilities location in the plane using the gravity model. Geographical Analysis 2006;38:391–406. [15] Drezner T, Drezner Z. The gravity p-median model. European Journal of Operational Research 2007;179:1239–51. [16] Drezner T, Drezner Z, Eiselt HA. Consistent and inconsistent rules in competitive facility choice. Journal of the Operational Research Society 1996;47:1494–503. [17] Drezner T, Eiselt HA. Consumers in competitive location models. In: Drezner Z, Hamacher HW, editors. Facility location: applications and theory. Berlin: Springer-Verlag; 2002. p. 151–78. [18] Drezner Z. Competitive location strategies for two facilities. Regional Science and Urban Economics 1982;12:485–93. [19] Drezner Z, Wesolowsky GO. Location-allocation on a line with demanddependent costs. European Journal of Operational Research 1996;90:444–50. [20] Drezner Z, Wesolowsky GO. Allocation of demand when cost is demanddependent. Computers and Operations Research 1999;26:1–15.
701
[21] Drezner Z, Wesolowsky GO. Allocation of discrete demand with changing costs. Computers and Operations Research 1999;26:1335–49. [22] Francis RL, Lowe TJ, Rayco MB, Tamir A. Aggregation error for location models: survey and analysis. Annals of Operations Research 2009;167: 171–208. [23] Francis RL, Lowe TJ, Tamir A. Aggregation error bounds for a class of location models. Operations Research 2000;48:294–307. [24] Ghosh A, Rushton G. Spatial analysis and location-allocation models. New York, NY: Van Nostrand Reinhold Company; 1987. [25] Glover F. Future paths for integer programming and links to artificial intelligence. Computers and Operations Research 1986;13:533–49. [26] Glover F, Laguna M. Tabu search. Boston, Massachusetts: Kluwer Academic Publishers; 1997. [27] Hakimi SL. Optimum locations of switching centres and the absolute centres and medians of a graph. Operations Research 1964;12:450–9. [28] Hakimi SL. Optimum distribution of switching centers in a communication network and some related graph theoretic problems. Operations Research 1965;13:462–75. [29] Hakimi SL. On locating new facilities in a competitive environment. European Journal of Operational Research 1983;12:29–35. [30] Hillsman EL, Rhoda R. Errors in measuring distances from populations to service centers. Annals of Regional Science 1978;12:74–88. [31] Hotelling H. Stability in competition. Economic Journal 1929;39:41–57. [32] Huff DL. Defining and estimating a trade area. Journal of Marketing 1964;28:34–8. [33] Huff DL. A programmed solution for approximating an optimum retail location. Land Economics 1966;42:293–303. [34] Kirkpatrick S, Gelat CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [35] Leonardi G, Tadei R. Random utility demand models and service location. Regional Science and Urban Economics 1984;14:399–431. [36] Minieka E. The centers and medians of a graph. Operations Research 1977;25:641–50. [37] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts and applications of Voronoi diagrams. John Wiley: Wiley Series in Probability and Statistics; 2000. [38] Pasternack BA, Drezner Z. A note on calculating steady state results for an M/ M/k queuing system when the ratio of the arrival rate to the service rate is large. Journal of Applied Mathematics & Decision Science 1998;2:133–5. [39] Reilly WJ. The law of retail gravitation. New York, NY: Knickerbocker Press; 1931. [40] ReVelle C. The maximum capture or sphere of influence problem: hotelling revisited on a network. Journal of Regional Science 1986;26:343–57. [41] Sheffi Y. Urban transportation networks: equilibrium analysis with mathematical programming methods. Englewood Cliffs, NJ: Prentice-Hall; 1984. [42] Suzuki A, Okabe A. Using Voronoi diagrams. In: Drezner Z, editor. Facility location: a survey of application and methods. New York: Springer; 1995. p. 103–18. [43] Wang Q, Batta R, Rump C. Algorithms for a facility location problems with stochastic customer demand and immobile servers. Annals of Operations Research 2002;111:17–34. [44] Wilson AG. Retailers’ profits and consumers’ welfare in a spatial interaction shopping mode. In: Masser I, editor. Theory and practice in regional science. London: Pion; 1976. p. 42–59.
Contents lists available at ScienceDirect
Computers & Operations Research journal homepage: www.elsevier.com/locate/caor
The gravity multiple server location problem Tammy Drezner, Zvi Drezner Steven G. Mihaylo College of Business and Economics, California State University-Fullerton, Fullerton, CA 92834, United States
a r t i c l e in f o
a b s t r a c t
Available online 22 August 2010
Models for locating facilities and service providers to serve a set of demand points are proposed. The number of facilities is unknown, however, there is a given number of servers to be distributed among the facilities. Each facility acts as an M/M/k queuing system. The objective function is the minimization of the combined travel time and the waiting time at the facility for all customers. The distribution of demand among the facilities is governed by the gravity rule. Two models are proposed: a stationary one and an interactive one. In the stationary model it is assumed that customers do not consider the waiting time at the facility in their facility selection decision. In the interactive model we assume that customers know the expected waiting time at the facility and consider it in their facility selection decision. The interactive model is more complicated because the allocation of the demand among the facilities depends on the demand itself. The models are analyzed and three heuristic solution algorithms are proposed. The algorithms were tested on a set of problems with up to 1000 demand points and 20 servers. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Facility location Heuristics Multiple facilities
1. Introduction We consider models for locating facilities and allocating p servers among them so as to minimize customers’ travel time plus waiting time at the facilities. The number of facilities is a variable limited to a maximum of p facilities. We allow the location of one or more servers at each facility. Wang et al. [43] considered the problem of locating facilities that are modeled as M/M/1 queuing systems. Customers travel to the closest facility and a constraint is placed on the maximum expected waiting times at all facilities. The objective is to minimize the combined travel and waiting times. Berman et al. [5] considered a similar model with additional constraints on lost demand due to congestion and insufficient coverage. The objective is to minimize the number of facilities. Berman and Drezner [4] extended the model suggested by Wang et al. [43] by allowing more than one server to be located at each facility thus employing the M/M/k queuing system. Berman and Drezner [4] assume that customers patronize the closest facility. Aboolian et al. [1] investigated the Berman and Drezner [4] model applying the minimax objective rather than minimum. Typical to all of these models is the assumption that service is provided by the closest facility. This assumption is used in many location-allocation models such as the p-median, p-center or competitive models [27,28,24,36,6,31,18,29,40]. The proximity assumption is logical because it is in the user’s best interest to use Corresponding author. Tel.: + 1 714 657 278 2712.
E-mail address: [email protected] (Z. Drezner). 0305-0548/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2010.08.006
the services of the closest facility. This is indeed the best choice if the demand is allocated to facilities by a planner rather than by customers’ free choice and/or when facilities are equally attractive. Under the proximity assumption, customers use only one facility to the exclusion of all others, resulting in an ‘‘all or nothing’’ behavior. In reality there are many reasons why the proximity assumption, which leads to the ‘‘all or nothing’’ customer behavior, does not accurately describe customers’ behavior: (i) A small change in the location of facilities may discontinuously shift the whole demand at a demand point from one facility to another. Suppose that the closest facility is at a distance of 4.9 miles and the second closest one is at a distance of 5 miles. The whole demand is allocated to the first facility. However, if the first facility is relocated to a location at a distance of 5.1 miles from the demand point, the whole demand will now be allocated to the second facility. This does not reflect actual customer behavior. (ii) Different facilities have different levels of attractiveness and customers are willing to travel an extra distance to a farther and more attractive facility. This issue in competitive location models is addressed by Drezner [8,9] and in various utility (additive or multiplicative) models [35,11,16]. (iii) Even if the facilities are equally attractive, it is not clear that all customers at a demand point patronize the same facility. Customers do not necessarily measure the exact distance when deciding which facility to patronize for service. They select the facility that they perceive to be the closest to them.
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
Distances are likely perceived differently by different customers. Therefore, the selected facility to patronize will not be the same one for all customers at a demand point. Distance perception is discussed in Drezner and Eiselt [17]. (iv) Demand points are not mathematical points. Rather, they represent an area inhabited by customers. It is possible that the facility closest to customers residing on one side of the area is not the closest to customers residing on the other side of the same area. This is referred to as Type C error in the aggregation literature [30,23,22]. This issue can be illustrated by a Voronoi diagram based on the facilities [37,42]. A Voronoi diagram is a tessellation of the plane into polygons. Each polygon consists of all the points closest to one facility which must be located inside that polygon. The Voronoi diagram may partition a demand area into several parts, each closest to a different facility. Incorporating area of demand points in competitive location gravity-based models is addressed in Drezner and Drezner [12]. It is therefore more accurate and more realistic to employ the gravity rule proposed by Reilly [39] and later used by Huff [32,33] to explain customers’ spatial choice behavior and estimating buying power attracted to competing facilities. According to the gravity rule, the probability that a customer patronizes a facility is proportional to the facility’s attractiveness, and to a decreasing function of the distance to the facility, commonly known as the distance decay function. Huff [32,33] suggested a power of the distance as a decay function. Others suggested an exponent of the distance [44]. Drezner [10] showed empirically that an exponent of the distance more accurately describes customers’ behavior than a power of the distance. The gravity rule rectifies the four concerns raised above about the proximity rule. Some recent papers have re-analyzed well known models using the gravity rule rather than the proximity rule. For example, Drezner and Drezner [13] formulated and solved the hub location problem using the gravity rule. Drezner and Drezner [15] analyzed and solved p-median problems on a network and Drezner and Drezner [14] solved the p-median in the plane when allocation of customers among the facilities is governed by the gravity rule. The contribution of the current paper is two fold: The gravity rule is used for modeling, and two models are proposed, analyzed, and solved. The stationary model] following the one suggested by Berman and Drezner [4] but assuming that customers patronize a facility according to the gravity rule. The interactive model assuming that customers select the facility by incorporating the expected waiting time at the facility in the selection criteria. Customers may experience a long waiting time at a facility in a certain visit, thus selecting another facility in a following trip. The selected facility may be farther away but have shorter expected waiting time. This model is more difficult to analyze and solve because the attractiveness of the facility depends on the demand. Applications for such models include the location of automatic teller machines (ATMs) where at each facility one or more ATMs are installed; location of bank branches where a decision is required about the number of teller stations to be built at each branch; location of gas stations with the number of pumps in each gas station being a variable; location of post office branches, department of motor vehicle offices, or similar service establishments where the number of windows (servers) is also a variable. The paper is organized as follows: In the next section we formulate, analyze, and propose solution methods for the stationary model. In Section 3 we define and analyze the interactive model. In Section 4 we report computational experiments on a
695
set of test problems for both the stationary and the interactive models. We conclude in Section 5 and propose ideas for future research.
2. The stationary model The analysis of the stationary model follows the approach in Berman and Drezner [4]. We briefly summarize the changes required for implementing the approach in Berman and Drezner [4]. 2.1. Notation Let N wi
be a set of demand points of cardinality n be the mean demand generated at demand point iA N in a Poisson process M be a set of potential locations for facilities of cardinality m be the attractiveness of a facility located at potential Aj location j A M p be the total number of servers available for installation at the facilities, one or more servers at each facility S D M of unknown cardinality s rp be an unknown set where facilities are located be the distance between demand point i A N and dij potential location j A M measured in time units by any distance measure f(d) be the decay function of the distance governing the gravity rule Wk ðl, mÞ be the expected time spent in line at an M/M/k system consisting of k servers with arrival rate l and service rate
m Lk ðl, mÞ kj Z 1 K
¼ lWk ðl, mÞ be the average length of the line in an M/M/ k system of k servers be the unknown number of servers allocated to facility jAS be the unknown vector K ¼ fk1 , . . . ,ks g of the allocation of servers to the set of facilities S
2.2. Formulation The average demand lj ðSÞ for j A S is by the gravity rule:
lj ðSÞ ¼
X
wi P
iAN
Aj f ðdij Þ
ð1Þ
j A S Aj f ðdij Þ
where f ðdÞ ¼ ecd (or any other distance decay function). The objective function, to be minimized, is the combined travel time and expected time spent at the selected facility for all customers. Note that the total travel time TðSÞ to the facilities is independent of the servers allocation: TðSÞ ¼
X iAN
wi
X Aj f ðdij Þ dij P j A S Aj f ðdij Þ jAS
ð2Þ
The total waiting time at the facilities, WðS,KÞ, is X X WðS,KÞ ¼ lj ðSÞWkj lj ðSÞ, m ¼ Lkj lj ðSÞ, m jAS
ð3Þ
jAS
The objective function, to be minimized, is FðS,KÞ ¼ TðSÞ þWðS,KÞ
ð4Þ
696
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
2.3. Analysis of the stationary model
3. The interactive model
We first devise an algorithm that, for a given set S DM, finds the optimal allocation of p servers among the members of S. Each member in S is allocated at least one server.
Suppose that customers have information about the expected waiting time at each facility. Customers will base their selection on the combined travel and waiting time. This entails replacing dij in (1) with dij þ Wkj ðlj ðSÞ, mÞ when calculating the allocation of customers to facilities. If the decay function f ðdÞ ¼ ecd is used, it
2.3.1. Calculating the total time spent at a facility This algorithm is based on Pasternack and Drezner [38] and is also given in Berman and Drezner [4]. Consider an M/M/k facility with mean arrival rate l and mean service rate m. Calculate the sequence far g for r ¼ 1, . . . ,k: a1 ¼ 1;
ar þ 1 ¼ 1 þr
m l
ar
ð5Þ
then Lk ðl, mÞ ¼
l2
ðkmlÞ l þ ðkmlÞak
cWk ðlj ðSÞ, mÞ
means replacing ecdij with fe be interpreted as defining
l ðkmlÞ½l þ ðkmlÞak
2.3.2. Finding the best allocation of servers The best allocation of servers is found by the following algorithm [4]. 1. For each node j A S calculate the minimum of servers number necessary to service the demand kj ¼ int lj ðSÞ=m þ 1. P 2. If sj ¼ 1 kj 4 p, there is no feasible solution. Ps 3. If j ¼ 1 kj ¼ p, the solution is fkj g. P P 4. If sj ¼ 1 kj o p, repeat the following until sj ¼ 1 kj ¼ p. (a) Check the decrease in Lkj (6) by changing kj to kj þ 1 for j ¼1,y,s. (b) Increase the number of servers by one at the facility with the maximum decrease Lkj Lkj þ 1 . Note that changing kj to kj þ 1 is very simple once the vector fakj g is saved. Since akj is known, akj þ 1 is calculated by (5) and then Lkj þ 1 is calculated by (6). The value of akj for the facility with the increased number of servers is updated. Since the best allocation of servers to facilities is optimally found by the approach above, the function FðSÞ ¼ minfFðS,KÞg K
ð7Þ
is well defined and can be calculated for any set S D M.
2.4. Solving the stationary problem The analysis in Section 2.3.2 enables us to calculate the optimal value of the objective function for a given set S D M of cardinality sr p. The optimization problem is therefore reduced to finding the set S that minimizes FðSÞ (7). Berman and Drezner [4] constructed a descent algorithm, a tabu search algorithm, a simulated annealing algorithm, and a genetic algorithm. We constructed the first three algorithms to heuristically solve the stationary problem of finding the best subset S D M. The genetic algorithm developed by Berman and Drezner [4] required about five times the computer time than the other metaheuristics and therefore was not tested. The heuristic algorithms are detailed in the Appendix.
ð8Þ
j
where Aj(S) is the attractiveness of facility j as is applied in competitive facilities location models ([8,9,2]) but it is not fixed. Rather, it is a function of the allocation of demand which is determined by the set S. Equation (1) is then replaced by
lj ðSÞ ¼ ð6Þ
gecdij . This formula can
cWk ðlj ðSÞ, mÞ
Aj ðSÞ ¼ e
X iAN
Wk ðl, mÞ ¼
j
wi P
Aj ðSÞf ðdij Þ
ð9Þ
j A S Aj ðSÞf ðdij Þ
where Aj(S) is defined by (8). We note that using the number of servers at a facility as the attractiveness measure is not appropriate. It is possible to have longer lines at a facility with more servers because there may be higher demand at such a facility. The difficulty with solving such a model is that the ‘‘attractiveness’’ Aj(S) is a function of the allocation of demand and is not fixed. The same issue has risen in other papers that analyzed models where the allocation of demand depends on the demand itself [19–21]. Consider the problem of minimizing F(S,K) for a given S. If we are able to find minK fFðS,KÞg, then the heuristic approaches designed to find the best set S for solving the stationary problem can be employed for solving the interactive problem. This problem for a given set S DM can be represented as a traffic flow problem [7,41]. The following network is designed (see Fig. 1). One set of nodes is N and a second set of nodes is S. An additional sink node is also defined. Traffic (demand) wi is generated at each iA N and the destination of each traveler is at the sink. A link connects every demand point iA N with a facility j A S. The length of the link is dij and does not depend on the flow on the link. Each facility j A S is connected to the sink. The length of such a link is the waiting time at the facility, which depends on the traffic on the link (congestion). A customer selects a facility, travels to it, and then uses the link connecting the facility to the sink which is the waiting time at the facility. The waiting time depends on the volume of traffic Demand
dij
Facilities
Wk Sink
Fig. 1. The traffic flow representation.
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
(total demand) on the link connecting the facility and the sink. Such a traffic flow problem has an equilibrium [7,41].
3.1. Solving the traffic flow problem
697
general, using small values of y will always lead to convergence, but may not be the most efficient.
3.2. Solving the interactive model
We follow the widely used methods for solving traffic flow problems [41]. Consider a given set S D M. For a given attractiveness levels Aj(S), for example Aj ðSÞ ¼ 1 8j, the values of lj ðSÞ are calculated. The best solution K can be found as detailed in Section 2.3.2. The values of Lk and Wk change. Therefore, the attractiveness levels Aj(S) change, and the problem is re-solved, and so on. As is experienced in traffic flow problems, such an iterative algorithm may not converge to an equilibrium. We therefore dampen the oscillations in the values of Aj ðSÞ by changing Aj(S) by a proportion 0 o y r 1, for example y ¼ 0:1, of the change to the new one. This means that if Auj ðSÞ are the values calculated by (8), Aj(S) is replaced by Aj ðSÞ þ yðAuj ðSÞAj ðSÞÞ. This can be interpreted as assuming that only 10% of the customers are cognizant of the change in attractiveness and modify their behavior. In the network flow algorithm an incremental adjustment of the traffic flow rather than full adjustment leads to an equilibrium. In our experiments using y ¼ 0:1 always converged to an equilibrium. In
We propose to apply the same metaheuristic approaches detailed in the Appendix by applying the solution approach described in Section 3.1 for finding K to define F(S). Note that the solution approach in Section 3.1 finds an equilibrium solution which may be a local minimum.
3.3. An extension to the interactive model Suppose that only a proportion 0 r y r 1 of the customers are cognizant of the waiting time at the facilities and consider it in their facility selection. This entails replacing dij in (1) with dij þ yWkj ðlj ðSÞ, mÞ when calculating the allocation of customers to facilities. The attractiveness is thus Aj ðSÞ ¼ e
ycWk ðlj ðSÞ, mÞ j
Table 1 Optimal solutions for the stationary model for a given number of facilities. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 a
Not run.
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
% over best known (Table 2)
Run rime (min)
s ¼1
s ¼2
s¼3
s¼4
s ¼5
s ¼1
s¼ 2
s¼3
s¼4
s¼5
64.71 69.41 76.29 75.35 79.49 92.25 100.07 99.83 94.96 119.70 118.18 122.34 101.07 133.95 134.15 135.93 102.24 146.60 151.24 146.60 101.72 152.27 161.65 156.20 103.54 158.91 170.56 169.43 103.41 161.62 178.62 181.43 106.39 167.35 189.83 190.14 105.43 168.44 194.11 200.22
27.93 28.37 31.88 30.73 38.18 47.68 53.24 53.11 50.21 68.02 66.86 69.57 54.43 79.02 78.21 79.52 53.54 85.16 88.14 84.41 55.16 90.26 97.05 92.94 58.61 97.43 106.22 105.23 58.98 101.53 114.62 116.60 61.27 105.87 122.79 123.28 60.89 107.08 126.42 131.47
6.91 9.14 10.42 2.73 16.80 19.22 20.57 20.11 24.48 33.31 30.53 32.60 24.91 44.28 42.32 42.39 24.99 51.61 51.75 47.76 27.86 55.84 61.78 57.86 30.53 59.36 66.67 65.47 30.44 60.90 71.27 72.74 32.98 64.14 78.04 77.82 32.48 64.56 80.87 83.76
0.00 1.68 1.48 0.40 5.71 9.30 7.37 6.36 12.89 19.22 13.36 12.86 14.35 21.30 19.96 20.37 14.16 25.90 27.54 24.17
1.28 0.35 0.00 0.00 0.00 4.89 4.25 3.25 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.03 0.03 0.03 0.03 0.05 0.05 0.05 0.05 0.07 0.07 0.07 0.07 0.10 0.10 0.10 0.10 0.14 0.14 0.14 0.14
0.01 0.01 0.01 0.01 0.09 0.09 0.09 0.09 0.46 0.47 0.47 0.48 1.51 1.52 1.53 1.54 3.84 3.85 3.88 3.90 7.85 7.88 7.92 7.96 14.56 14.62 14.69 14.77 24.84 24.83 24.99 25.08 39.64 39.71 39.87 40.00 60.05 60.18 60.38 60.56
0.17 0.17 0.17 0.18 5.43 5.43 5.48 5.56 41.71 41.74 42.00 42.42 181.55 183.07 183.94 183.92 565.57 566.83 569.87 575.04
3.93 3.94 3.96 4.03 257.30 258.41 260.54 260.42 2987.53
a a a a a a a a a a a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a a a a a a a a a a a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
698
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
When y ¼ 0 the model is stationary. Otherwise, the methods proposed for the solution of the interactive model can be used for the solution of this problem with y 4 0.
demand resulting in problems which are easier to solve. The exponential decay function used is f ðdÞ ¼ e4d . 4.1. The stationary model results We attempted to require about the same total run times for the three algorithms tested. The descent algorithm was replicated 2000 times. Following extensive experiments the following parameters were selected. For the simulated annealing algorithm: the initial temperature T0 was set to 100; the number of iterations was set to N¼1000 np, and the temperature reduction factor was set to a ¼ 15=N. For the tabu search: the tabu tenure was randomly selected each iteration in the range ½0:1n,0:2n, and the number of iterations was set to 100 times the number of iterations in the descent phase. Each problem was solved ten times by both the simulated annealing algorithm and the tabu search. In some situations the planner wishes to have a given number of facilities. The number of servers is assumed to be p as in the original problem. By total enumeration of all possible sets S D M of cardinality 1, 2, or 3 the optimum for one, two, or three facilities was found. We also found the optimal solution for four and five
4. Computational experiments Programs were coded in Fortran using double precision arithmetic and compiled by the Compaq 6.6 Fortran compiler. The programs were run on a desk top Core 2 Duo PC, 2.3 GHz with 2 MB RAM. We randomly generated problems with n ¼ 100,200, . . . ,1000 points in a square of side 10. These points are both generating demand and potential locations for the facilities. Four values of p ¼5, 10, 15, 20 were tested for each value of n. Distances are assumed Euclidean. Each point generates demand wi A ð0,1Þ. P When a service rate of m ¼ ð1=pÞ ni¼ 1 wi is used, total demand is equal to the total service rate. The service rate at each server P was therefore calculated as m ¼ ð1:1=pÞ ni¼ 1 wi so it is 10% higher than total demand. This value of m defines difficult problems because there are many infeasible solutions. If one applies, for P example, m ¼ ð2=pÞ ni¼ 1 wi , capacity of the system is double the Table 2 Best solutions found for the stationary model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Best known
132.765 128.564 123.474 123.841 213.656 199.022 191.182 191.152 298.528 264.516 266.294 261.082 382.856 328.668 328.338 325.646 476.711 390.604 383.343 390.343 567.167 453.174 436.873 445.964 651.937 512.174 490.068 491.941 735.743 571.720 536.787 531.237 822.098 634.331 585.084 584.289 912.215 697.788 636.839 623.689
Descent
Sim. anneal.
Tabu search
a
b
a
b
a
b
23 0 49 2 0 0 7 0 6 0 0 0 3 0 0 0 15 0 0 0 8 0 0 0 61 0 0 0 22 0 0 0 100 0 0 0 30 0 0 0
0.00 0.35 0.00 0.00 0.18 2.67 0.00 1.25 0.00 2.98 0.63 1.61 0.00 3.20 1.75 2.78 0.00 6.50 6.44 2.00 0.00 1.29 9.35 3.08 0.00 1.77 10.95 3.47 0.00 0.96 12.98 6.98 0.00 1.37 13.36 7.39 0.00 0.58 15.25 10.06
9 4 5 2 3 2 1 4 10 1 0 1 7 1 1 0 10 3 1 1 4 2 1 0 10 3 1 1 8 3 1 1 10 3 1 1 10 2 1 1
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.12 0.00 0.00 0.00 0.00 0.84 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.77 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
10 0 3 2 2 2 7 0 10 0 1 0 7 0 0 1 7 0 0 0 9 1 0 1 10 0 0 0 9 1 0 0 10 0 0 0 10 1 0 0
0.00 0.35 0.00 0.00 0.00 0.00 0.00 0.15 0.00 2.06 0.00 1.03 0.00 0.61 0.45 0.00 0.00 0.30 2.94 0.68 0.00 0.00 2.83 0.00 0.00 0.24 2.68 2.45 0.00 0.00 2.00 5.00 0.00 0.59 1.59 4.04 0.00 0.00 0.92 2.69
8.15
Number of times best known solution found. Percentage of best found solution above best known solution.
3.28
3.25
0.09
2.60
0.84
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
facilities for some of the smaller problems. The results are summarized in Table 1. The best known solution for these problems is the best result obtained by the three heuristics reported in Table 2. We report the percentage of the optimal solution for a given number of facilities, s, above the best known solution for the particular problem, and the run time in minutes required for such a solution. If only one facility is allowed, the optimal value of the objective function can be more than three times the optimal solution (more than 200% above the best known solution). Note that for the problem with n ¼100 and p ¼5, s¼4 is optimal and for the problems with n¼ 200, 300 and p ¼5, s¼5 is optimal. The optimality of the best known solution for these three problems was thus proven. The performance of the descent, simulated annealing, and tabu search algorithms is reported in Tables 2 and 3. In Table 2 we report the best known solution, and the number of times (out of 2000 replications for the descent algorithm and 10 replications for the other two algorithms) that the best known solution was found, and the percentage of the best solution found by the corresponding algorithm above the best known solution. In Table 3 we report the percentage of the average solution for all
Table 3 Average solutions found for the stationary model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Sim. anneal.
Tabu search
a
Descent b
a
b
a
b
12.55 9.84 7.15 3.40 17.56 14.49 9.62 7.33 19.06 17.01 8.66 7.82 16.05 18.81 11.20 9.44 10.63 20.79 15.55 9.80 7.33 22.05 18.52 11.70 5.50 22.60 20.23 14.73 3.26 21.48 21.53 17.78 3.02 20.37 23.36 17.99 2.31 18.08 24.02 20.81
0.20 0.22 0.26 0.24 0.98 1.37 1.83 1.92 3.15 5.42 6.04 6.82 7.04 11.94 14.39 15.08 14.76 23.50 26.72 31.26 24.84 40.41 47.23 58.27 37.78 64.14 75.65 96.28 51.68 102.73 117.08 146.87 68.41 162.38 169.18 223.48 86.95 247.59 235.50 313.88
0.05 0.87 0.27 0.06 0.25 1.62 2.09 0.87 0.00 2.05 4.61 4.65 0.06 0.96 2.56 3.85 0.00 0.50 2.23 3.73 0.04 0.37 1.48 3.10 0.00 0.21 1.80 2.21 0.01 0.34 1.87 1.73 0.00 1.00 0.64 1.18 0.00 0.18 0.61 1.00
0.28 0.70 1.15 1.71 1.12 3.16 5.64 8.30 2.51 7.68 14.68 22.09 4.42 14.13 28.08 42.45 6.96 22.58 44.92 74.09 9.99 32.36 66.00 110.95 13.66 44.34 91.10 156.69 17.75 57.58 117.75 210.12 22.45 72.84 149.68 268.66 27.69 89.89 185.05 333.46
0.00 1.34 0.75 0.11 0.30 1.84 0.13 0.98 0.00 2.92 1.45 1.66 0.05 2.24 2.31 1.38 0.01 1.40 5.51 1.27 0.00 0.82 5.01 2.62 0.00 1.19 5.82 4.09 0.01 0.72 4.59 6.41 0.00 1.24 3.63 6.54 0.00 0.52 3.59 7.29
0.17 0.23 0.22 0.23 0.86 1.58 1.86 1.79 2.26 6.37 5.49 7.83 4.58 13.80 20.80 21.34 8.61 30.48 36.54 53.24 12.83 48.47 69.06 79.63 24.68 70.92 117.26 134.02 30.25 107.19 171.92 200.97 40.81 160.35 218.61 307.36 47.76 178.84 351.02 451.04
14.09
63.59
1.23
59.62
1.99
76.03
Percentage of average solution above best known solution. Total time in minutes.
699
replications above the best known solution for the corresponding algorithm and the total run time in minutes for all replications. By examining the results reported in Tables 2 and 3 we conclude that the simulated annealing algorithm performed much better than the other two algorithms for this set of problems and is the recommended algorithm for the stationary model. 4.2. The interactive model results We performed the same experiments for the interactive model. Since run time is longer, we experimented with n r500 problems and the descent algorithm was run only 1000 times. Following extensive experiments we also changed the number of iterations pffiffiffi for the simulated annealing algorithm to N ¼ 100np n. We did not repeat the given number of facilities experiment because the best K found by the iterative procedure is not proven optimal. Note that the s ¼1 results are identical to the results for the stationary model. We report in Tables 4 and 5 the performance measures that were reported in Tables 2 and 3 for the stationary model. By examining the results reported in Tables 4 and 5 we conclude that the tabu search algorithm performed better than the other two algorithms and is the recommended algorithm for the interactive model.
5. Conclusions Two models for locating facilities to serve a set of demand points were proposed. The number of facilities is unknown. However, there is a given number of servers to be distributed among the facilities. Each facility acts as an M/M/k queuing system. The objective function is the minimization of the combined travel time and the waiting time at the facility for all customers. The distribution of demand among the facilities is governed by the gravity rule. Two models were proposed. In the stationary model it is assumed that customers do not consider the expected waiting time at the facility in their selection decision. In Table 4 Best solutions found for the interactive model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500 Average a b
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Best known
129.618 125.115 121.609 120.376 208.688 194.731 189.612 187.333 296.184 259.033 260.194 253.655 380.856 323.206 319.310 315.332 474.793 385.589 374.061 374.990
Sim. anneal.
Tabu search
a
Descent b
a
b
a
b
49 0 127 39 106 0 220 17 172 4 0 0 27 3 0 0 130 7 0 0
0.00 0.25 0.00 0.00 0.00 0.80 0.00 0.00 0.00 0.00 0.21 0.28 0.00 0.00 1.52 0.54 0.00 0.00 2.23 0.93
10 1 9 8 10 0 1 4 10 3 0 0 10 1 0 0 10 4 0 0
0.00 0.00 0.00 0.00 0.00 0.27 0.00 0.00 0.00 0.00 0.64 0.08 0.00 0.00 0.08 0.24 0.00 0.00 0.21 0.09
10 5 10 10 10 1 7 10 10 8 2 3 10 6 2 2 10 9 1 1
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
45.05
0.34
4.05
0.08
6.35
Number of times best known solution found. Percentage of best found solution above best known solution.
0.00
700
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
Table 5 Average solutions found for the interactive model by various algorithms. n
p
100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400 500 500 500 500
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
Average a b
Descent
Sim. anneal.
Tabu search
a
b
a
b
a
b
5.01 4.48 1.99 2.36 7.04 7.07 2.44 3.27 3.82 7.78 4.59 4.30 2.60 7.89 6.75 5.33 1.88 6.06 10.04 5.29
1.49 3.13 8.28 11.80 6.75 15.87 36.19 55.27 19.06 49.06 79.06 129.13 37.14 102.70 157.37 230.19 59.45 197.68 260.96 393.20
0.00 0.23 0.01 0.01 0.00 0.76 0.06 0.15 0.00 0.12 1.42 0.29 0.00 0.17 0.64 0.91 0.00 0.09 0.49 0.38
1.73 5.93 15.78 31.85 8.32 19.61 44.54 74.66 23.06 53.03 99.55 161.84 47.16 102.30 186.01 294.88 76.30 173.22 301.94 502.47
0.00 0.13 0.00 0.00 0.00 0.48 0.02 0.00 0.00 0.05 0.73 0.16 0.00 0.09 0.41 0.50 0.00 0.03 0.28 0.74
1.66 3.59 10.08 14.04 7.32 17.04 31.33 44.70 20.39 44.96 52.63 84.62 38.15 107.49 121.23 173.69 66.23 228.54 244.87 301.53
5.00
92.69
0.29
111.21
0.18
80.71
Percentage of average solution above best known solution. Total time in minutes.
the interactive model we assume that customers know the expected waiting time at the facility and consider it in their selection decision. The interactive model is more complicated to analyze and solve because the allocation of the demand among the facilities depends on the demand itself. Three heuristic procedures for the solution of the two models were proposed: a descent algorithm, simulated annealing, and tabu search. The recommended approaches are simulated annealing for the stationary model and the tabu search for the interactive model. As future research we propose to investigate a host of related models such as an expansion of an existing system by opening new facilities and/or expanding the number of servers in existing facilities due to increase in demand, competitive models based on the framework proposed in this paper, and many variants of competitive location models as described in Berman et al. [3]. Note that the traffic flow algorithm can be employed for the interactive model based on the proximity assumption leading to another extension of the model in Berman and Drezner [4].
Appendix A The three heuristics are detailed in this appendix. A.1. The descent algorithm The descent algorithm is quite straightforward. A neighborhood of a set of sites S of cardinality s is defined by either removing a site from S (if s 4 1), or adding a site to S (if s o p), or adding a site to S and removing another one from S. The cardinality of such a neighborhood for 1 o s op is n þ sðnsÞ. 1. Start with a randomly generated S and calculate F(S). In our experiments we started with s¼2 because F(S) is feasible for any S of cardinality s¼1 and thus there is always a feasible solution in the neighborhood of the starting S.
2. Check all the subsets S0 in the neighborhood of S and calculate F(S0 ) for each. Find the minimum of these values Fmin and the subset Smin. 3. If Fmin oFðSÞ change S to Smin and go to Step 2. 4. Otherwise, the descent algorithm terminates with a solution S.
A.2. Simulated annealing The simulated annealing algorithm [34] simulates the cooling process of hot melted metals. The variant used in this paper depends on three parameters: the starting temperature T0, the number of iterations N, and the factor a o 1 by which the temperature is reduced every iteration. The neighborhood of S is defined in the descent algorithm. 1. Start with a randomly generated S and calculate F(S). Set the temperature T to T0. 2. A move in the neighborhood of S is randomly selected. Each move is selected with the same probability. The move defines S0 . 3. Calculate F(S0 ). 4. If FðSuÞ r FðSÞ change S to Su and go to Step 6. 5. Otherwise, (a) calculate d ¼ ðFðSuÞFðSÞÞ=T. (b) Change S to Su with probability ed . Otherwise, retain S. 6. Set T ¼ aT and go to Step 2 unless the predetermined number of iterations is reached. 7. The best encountered solution throughout the process is the resulting solution.
A.3. Tabu search A tabu search [25,26] proceeds from the point that the descent algorithm terminated. We establish a tabu list of ‘‘reverse’’ moves (adding back to S a site removed from S), and moves remain in the tabu list for a given tenure. Note that we just need to keep the last iteration number at which a site was removed from S for all the sites in M. 1. Apply the descent algorithm detailed above. Let S be the result of the descent algorithm. F* is set to F(S) and S* is set to S. Empty the tabu list. 2. Check all the subsets S0 in the neighborhood of S and calculate F(S0 ) for each. 3. Find the minimum of these values Fmin and the subset Smin. 4. If Fmin oF , change S and S to Smin , and go to Step 2. 5. Otherwise, (a) Find the minimum FðSuÞ for moves not in the tabu list (i.e., the list of sites that were removed from S during the tabu tenure). The minimal S0 , whether improving or not, is the S for the next iteration. (b) Add the facility removed from S (if any) to the tabu list. (c) Remove the sites whose tenure exceeds the tabu tenure from the tabu list. (d) If the stopping criterion is not met, go to Step 2. Otherwise, stop the search and S* is the resulting solution. We used a random tabu tenure each iteration as detailed in Section 4. To operationalize the tabu list, the last iteration when a site was removed from S is recorded for each j A M (a large negative number is recorded for sites that were never removed form S). A site is in the tabu list if and only if the difference between the current iteration number and its recorded last
T. Drezner, Z. Drezner / Computers & Operations Research 38 (2011) 694–701
iteration number is less than or equal to the tabu tenure. This is an easy way to handle varying tabu tenure values for each iteration. References [1] Aboolian R, Berman O, Drezner Z. The multiple server center location problem. Annals of Operations Research 2009;167:337–52. [2] Aboolian R, Berman O, Krass D. Efficient solution approaches for discrete multi-facility competitive interaction model. Annals of Operations Research 2009;167:297–306. [3] Berman O, Drezner T, Drezner Z, Krass D. Modeling competitive facility location problems: new approaches and results. In: Oskoorouchi M, editor. Tutorials in operations research. San Diego: INFORMS; 2009. p. 156–81. [4] Berman O, Drezner Z. The multiple server location problem. Journal of the Operational Research Society 2007;58:91–9. [5] Berman O, Krass D, Wang J. Locating service facilities to reduce lost demand. IIE Transactions 2006;38:933–46. [6] Current J, Daskin M, Schilling D. Discrete network location models. In: Drezner Z, Hamacher HW, editors. Facility location: applications and theory. Berlin: Springer-Verlag; 2002. p. 81–118. [7] Dafermos S. Traffic equilibrium and variational inequalities. Transportation Science 1980;14:42–54. [8] Drezner T. Locating a single new facility among existing unequally attractive facilities. Journal of Regional Science 1994;34:237–52. [9] Drezner T. Competitive facility location in the plane. In: Drezner Z, editor. Facility location: a survey of applications and methods. New York, NY: Springer; 1995. p. 285–300. [10] Drezner T. Derived attractiveness of shopping malls. IMA Journal of Management Mathematics 2006;17:349–58. [11] Drezner T, Drezner Z. Competitive facilities: market share and location with random utility. Journal of Regional Science 1996;36:1–15. [12] Drezner T, Drezner Z. Replacing discrete demand with continuous demand in a competitive facility location problem. Naval Research Logistics 1997;44:81–95. [13] Drezner T, Drezner Z. A note on applying the gravity rule to the airline hub problem. Journal of Regional Science 2001;41:67–73. [14] Drezner T, Drezner Z. Multiple facilities location in the plane using the gravity model. Geographical Analysis 2006;38:391–406. [15] Drezner T, Drezner Z. The gravity p-median model. European Journal of Operational Research 2007;179:1239–51. [16] Drezner T, Drezner Z, Eiselt HA. Consistent and inconsistent rules in competitive facility choice. Journal of the Operational Research Society 1996;47:1494–503. [17] Drezner T, Eiselt HA. Consumers in competitive location models. In: Drezner Z, Hamacher HW, editors. Facility location: applications and theory. Berlin: Springer-Verlag; 2002. p. 151–78. [18] Drezner Z. Competitive location strategies for two facilities. Regional Science and Urban Economics 1982;12:485–93. [19] Drezner Z, Wesolowsky GO. Location-allocation on a line with demanddependent costs. European Journal of Operational Research 1996;90:444–50. [20] Drezner Z, Wesolowsky GO. Allocation of demand when cost is demanddependent. Computers and Operations Research 1999;26:1–15.
701
[21] Drezner Z, Wesolowsky GO. Allocation of discrete demand with changing costs. Computers and Operations Research 1999;26:1335–49. [22] Francis RL, Lowe TJ, Rayco MB, Tamir A. Aggregation error for location models: survey and analysis. Annals of Operations Research 2009;167: 171–208. [23] Francis RL, Lowe TJ, Tamir A. Aggregation error bounds for a class of location models. Operations Research 2000;48:294–307. [24] Ghosh A, Rushton G. Spatial analysis and location-allocation models. New York, NY: Van Nostrand Reinhold Company; 1987. [25] Glover F. Future paths for integer programming and links to artificial intelligence. Computers and Operations Research 1986;13:533–49. [26] Glover F, Laguna M. Tabu search. Boston, Massachusetts: Kluwer Academic Publishers; 1997. [27] Hakimi SL. Optimum locations of switching centres and the absolute centres and medians of a graph. Operations Research 1964;12:450–9. [28] Hakimi SL. Optimum distribution of switching centers in a communication network and some related graph theoretic problems. Operations Research 1965;13:462–75. [29] Hakimi SL. On locating new facilities in a competitive environment. European Journal of Operational Research 1983;12:29–35. [30] Hillsman EL, Rhoda R. Errors in measuring distances from populations to service centers. Annals of Regional Science 1978;12:74–88. [31] Hotelling H. Stability in competition. Economic Journal 1929;39:41–57. [32] Huff DL. Defining and estimating a trade area. Journal of Marketing 1964;28:34–8. [33] Huff DL. A programmed solution for approximating an optimum retail location. Land Economics 1966;42:293–303. [34] Kirkpatrick S, Gelat CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [35] Leonardi G, Tadei R. Random utility demand models and service location. Regional Science and Urban Economics 1984;14:399–431. [36] Minieka E. The centers and medians of a graph. Operations Research 1977;25:641–50. [37] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts and applications of Voronoi diagrams. John Wiley: Wiley Series in Probability and Statistics; 2000. [38] Pasternack BA, Drezner Z. A note on calculating steady state results for an M/ M/k queuing system when the ratio of the arrival rate to the service rate is large. Journal of Applied Mathematics & Decision Science 1998;2:133–5. [39] Reilly WJ. The law of retail gravitation. New York, NY: Knickerbocker Press; 1931. [40] ReVelle C. The maximum capture or sphere of influence problem: hotelling revisited on a network. Journal of Regional Science 1986;26:343–57. [41] Sheffi Y. Urban transportation networks: equilibrium analysis with mathematical programming methods. Englewood Cliffs, NJ: Prentice-Hall; 1984. [42] Suzuki A, Okabe A. Using Voronoi diagrams. In: Drezner Z, editor. Facility location: a survey of application and methods. New York: Springer; 1995. p. 103–18. [43] Wang Q, Batta R, Rump C. Algorithms for a facility location problems with stochastic customer demand and immobile servers. Annals of Operations Research 2002;111:17–34. [44] Wilson AG. Retailers’ profits and consumers’ welfare in a spatial interaction shopping mode. In: Masser I, editor. Theory and practice in regional science. London: Pion; 1976. p. 42–59.