- Email: [email protected]
The minimum equitable radius location problem with continuous demand
Available online at www.sciencedirect.com
European Journal of Operational Research 195 (2009) 17–30 www.elsevier.com/locate/ejor
Continuous Optimization
The minimum equitable radius location problem with continuous demand Atsuo Suzuki a, Zvi Drezner b,* b
a Department of Mathematical Sciences, Nanzan University, 27 Seirei, Seto, Aichi 489-0863, Japan College of Business and Economics, California State University-Fullerton, Fullerton, CA 92834, United States
Received 6 September 2007; accepted 10 January 2008 Available online 26 January 2008
Abstract We analyze the location of p facilities satisfying continuous area demand. Three objectives are considered: (i) the p-center objective (to minimize the maximum distance between all points in the area and their closest facility), (ii) equalizing the load service by the facilities, and (iii) the minimum equitable radius – minimizing the maximum radius from each point to its closest facility subject to the constraint that each facility services the same load. The paper offers three contributions: (i) a new problem – the minimum equitable radius is presented and solved by an efficient algorithm, (ii) an improved and efficient algorithm is developed for the solution of the p-center problem, and (iii) an improved algorithm for the equitable load problem is developed. Extensive computational experiments demonstrated the superiority of the new solution algorithms. Ó 2008 Elsevier B.V. All rights reserved. Keywords: Location; p-Center; Equitable loads; Voronoi diagrams
1. Introduction Recently, problems of locating p facilities to satisfy continuous demand in an area were investigated. Ohya et al. (1984) formulated and solved the p-median problem in an area by applying an iterative procedure on the Voronoi diagram of the facilities (Okabe et al., 2000; Suzuki and Okabe, 1995). A similar approach was utilized by Okabe and Suzuki (1987) to formulate and solve the competitive location problem in an area and by Suzuki and Drezner (1996) to formulate and solve the p-center problem in an area. Baron et al. (2007) suggested an equitable load model (EL) for uniform continuous demand in an area. Their model seeks the location of p facilities in the area, each point in the area is serviced by the closest facility as long as it is within a given coverage distance r from the facility. The objective is equalizing the loads serviced by each facility which is formulated as minimizing the maximum load serviced by the facilities. If the value of r is too small, there may be no feasible solution. For a large enough value of r, the optimal solution is equal in areas for all facilities. They solved the EL problem by an iterative procedure based on Voronoi diagrams. The question of finding the minimum possible r among all optimal solutions to the equitable load problem, is a practical and interesting problem. We term this problem ‘‘the minimum equitable radius” (MER) problem. The EL solution is usually not unique. Therefore, the maximum distance to the nearest facility may vary from one optimal solution to another. The MER problem is to find the minimum value among those maximum distances. To illustrate the MER problem consider the location of four facilities in a square. If the four facilities are placed equally spaced on the periphery of a circle centered at the center of the square (the radius of the circle is irrelevant), they divide the *
Corresponding author. Tel.: +1 714 278 2712; fax: +1 714 278 5940. E-mail address: [email protected] (Z. Drezner).
0377-2217/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2008.01.022
18
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
square into four equal areas (see Fig. 1). This is because the perpendicular bisector between two facilities must pass through the center of the square, and the four areas are thus congruent. Each angle h and any radius lead to the locations of four facilities which are optimal solutions to the EL problem. Therefore, the number of EL solutions is infinite. It is easy to pffiffi show that the radius enclosing each shape is 2 cos 2ph . This radius is obtained when the facilities are located at the centers ð4 Þ of the lines connecting the intersection points of the perpendicular bisectors with the sides of the square. The optimal solution of the MER problem (the minimum possible radius), as proven in this paper, is obtained for h ¼ p4 when the facilities are located on the diagonals at the centers of the lines connecting the center of the square and its vertices. The square is divided into four equal small squares. The p-center problem is investigated in other environments as well. In the plane, it is investigated by Drezner (1984) and Vijay (1985) for Euclidean distances and by Drezner (1987) for rectilinear distances. In a network environment the p-center problem is extensively investigated (for example, Handler, 1990), see, Daskin (1995) and Current et al. (2002) for a review. The EL problem with no coverage radius constraints was investigated by Drezner and Drezner (2006) for planar problems and by Berman et al. (in press) in the network environment. We are not aware of any paper suggesting the MER problem. There are many situations in which equal loads of the facilities is desirable. When cell phone towers are established in an area, the busiest tower determines the efficiency of the system thus equalizing loads are desirable. This applies to any design of p identical M/M/1 facilities where minimizing the maximum waiting time at the facilities is equivalent to minimizing the maximum arrival rate at the facilities which is proportional to the number of customers getting service at the facility. When carving market territories to be assigned to marketing representatives, the objective is to assign to each territory an equitable market potential. Another example is the problem of designing machines with similar capacities in a production facility. The throughput of the system depends on the machine with the maximum load. Voter districting is another application. The objective is that each district will have about the same number of voters. The maximum distance to the closest facility is a measure of good service used in numerous p-center applications. The MER model combines the two objectives. We seek the minimum possible radius among all optimal solutions to the EL model thus achieving both the EL and the p-center objectives. In every application described above for the EL model, it is desirable to minimize the maximum distance to the closest facility in addition to equalizing the loads as much as possible. For example, it is desirable that voters in a district are close to the voting place. It is advantageous for customers to be close to the facilities. In this paper, we propose three contributions. First, we propose a better heuristic algorithm to solve the p-center problem in an area. Second, pwe ffiffiffi propose a better heuristic algorithm for the solution of the EL problem for a large enough value of r (for example r P 2 for a square area as suggested by Baron et al., 2007) so that the distance constraints are not active in the solution. Third, we propose and heuristically solve the MER problem of finding the smallest possible radius among all equitable loads solutions. The problems are defined for any demand area and any demand distribution. However, most of the analysis in this paper pertains to uniform demand and a square area. It can be generalized to any convex polygonal shape and non-uniform demand. Finding the Voronoi diagram inside a polygonal area is quite easy. If demand is not uniform, calculating the
θ
Fig. 1. EL and MER solutions.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
19
necessary derivatives may require numerical differentiation. Note that the distribution of demand is irrelevant for the p-center problem as long as demand’s density function is greater than zero at every point in the area. 2. Discussion A Voronoi diagram defined by a set of generator points (for example, facilities) is a partition of the area into Voronoi regions so that all the points in a region are closest to a generator point, see, Figs. 2–5. Each generator point defines its associated Voronoi area which is a polygon surrounding it. The sides of the polygon are perpendicular bisectors between generator points. If demand is uniform, then the load of each facility is the area of the Voronoi region associated with the facility. Baron et al. (2007) showed that, in a square area, for every p (if r is large enough so that every demand point is covered) there exists an optimal solution by creating p ‘‘thin” rectangles, whose length is ‘‘1”, with the same width of 1p side by side
MER Solution
p-Center Solution Fig. 2. p ¼ 3 solutions.
MER Solution
p-Center Solution Fig. 3. p ¼ 5 solutions.
20
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
MER Solution
p-Center Solution Fig. 4. p ¼ 7 solutions.
MER Solution
p-Center Solution Fig. 5. p ¼ 8 solutions.
and placing a facility at the center of each rectangle. This configuration yields p equal loads and thus must be optimal. Since we know that there exists an optimal solution for which all loads are equal, the MER problem must have a solution as well. We conjecture that for every p, when r is large enough, there is a continuum of optimal solutions. If our conjecture is true, then finding the smallest possible r among the continuum of optimal solutions (the MER problem) can be a challenging one. The following two theorems provide insight into the MER problem and its relationship to the p-center problem: Theorem 1. The minimum radius for the p-center solution is a lower bound for the MER solution. Proof. Every feasible solution for the MER problem is a feasible solution to the p-center problem. Therefore, the optimal solution to the MER problem cannot be better than the optimal solution to the p-center problem. h In the rest of the paper we assume that the MER problem is defined in a square of side 1.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Theorem 2. The optimal solution to the MER problem is lower than or equal to Proof. The radius enclosing the ‘‘thin” rectangles described above is
0:5 p
0:5 p
21
pffiffiffiffiffiffiffiffiffiffiffiffiffi p 2 þ 1.
pffiffiffiffiffiffiffiffiffiffiffiffiffi p 2 þ 1. h
Note that by Theorem 1 the solution to the p ¼ 4 problem is to split the square into four small squares because this is the optimal solution to the p-center problem (Suzuki and Drezner, 1996) and it is feasible for the MER problem. However, the optimal solution to the p-center problem for p ¼ 9 is not dividing the square into 9 small squares, which yields a radius of 0.23570. A solution with r ¼ 0:23064 is given by Suzuki and Drezner (1996). Therefore, the nine small squares solution may not be the optimal solution to the MER problem. The MER solution and the p-center solution solve the EL problem (Baron et al., 2007) for most values of r. If the optimal radius of the MER problem does not exceed the given r, it is a solution to the EL problem. If the optimal radius to the MER problem exceeds the given r, then there is no solution to the EL problem with equal loads. If r is less than the p-center solution, then there is no feasible solution to the EL problem. The EL problem remains of interest only for r between the Table 1 Comparing p-center and MER solutions p
MER
p-center
Diff.
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50714 0.35355 0.34874 0.30046 0.29243 0.26909 0.23570 0.23702 0.21917 0.20833 0.20480 0.19339 0.18623 0.17677 0.17780 0.16976 0.16417 0.15986 0.15653 0.15266 0.14853 0.14520 0.14142 0.14012 0.13764 0.13369 0.13190 0.12954 0.12736 0.12446 0.12311 0.12130 0.11953 0.11785 0.11800 0.11507 0.11392 0.11238 0.11079 0.10952 0.10894 0.10723 0.10515 0.10444 0.10317 0.10206 0.10102 0.10073
0.50389 0.35355 0.32616 0.29873 0.27429 0.26030 0.23064 0.21823 0.21252 0.20228 0.19431 0.18551 0.17966 0.16943 0.16568 0.16064 0.15818 0.15225 0.14895 0.14369 0.14124 0.13830 0.13446 0.13223 0.13087 0.12761 0.12616 0.12204 0.12208 0.11841 0.11789 0.11498 0.11336 0.11168 0.11110 0.10856 0.10770 0.10702 0.10508 0.10400 0.10234 0.10126 0.10028 0.09922 0.09793 0.09639 0.09558 0.09451
0.00325 0.00000 0.02258 0.00173 0.01814 0.00879 0.00506 0.01879 0.00666 0.00605 0.01049 0.00788 0.00657 0.00735 0.01212 0.00912 0.00600 0.00761 0.00758 0.00896 0.00729 0.00690 0.00696 0.00788 0.00677 0.00608 0.00574 0.00750 0.00528 0.00604 0.00522 0.00632 0.00617 0.00618 0.00690 0.00651 0.00622 0.00536 0.00571 0.00552 0.00661 0.00597 0.00487 0.00522 0.00524 0.00567 0.00545 0.00622
22
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
p-center solution and the MER solution for which there is no solution to the EL problem with equal areas. In Table 1, the solutions found in the computational experiments section below for the MER and p-center problems are summarized. For pffiffiffiffi example, assuming that the values in Table 1 are optimal, for p ¼ 3: if r < 0:50389 ð¼ 1665 which is optimal by Suzuki and Drezner, 1996), there is no feasible solution to the EL problem, and if r P 0:50714, the MER solution is optimal for the EL problem. In the range 0:50389 6 r < 0:50714, a range of only 0.00325, the solution to the EL problem is not with equal areas and the algorithms in Baron et al. (2007) need to be employed to find such a solution. These two solutions are depicted in Fig. 2. The p-center solution is the minimum possible radius, and the MER solution is the minimum possible radius among all partitions of the square into three equal areas. Note that the ranges for requiring such a solution are quite small for all values of p depicted in Table 1. 3. Solution algorithms A vector X of p generator points is defined as X ¼ fðxi ; y i Þg; i ¼ 1; . . . ; p. Each point defines the closest Voronoi region V i ðX Þ (see Figs. 2–5). Such a Voronoi diagram defines a set of KðX Þ vertices. Let the set P i ðX Þ be the set of all Voronoi vertices of the Voronoi region defined by generator point i. Let also S i ðX Þ be the area of the Voronoi region (polygon) associated with generator point i and Rik ðX Þ for i ¼ 1; . . . ; p and k 2 P i ðX Þ be the Euclidean distance between generator point i and vertex k 2 P i ðX Þ belong to its Voronoi region. ik ðX Þ dRik ðX Þ ; dy for i ¼ 1; . . . ; p; k 2 P i ðX Þ; j ¼ 1; . . . ; p can be The derivatives dSdxi ðXj Þ ; dSdyi ðX Þ for i ¼ 1; . . . ; p; j ¼ 1; . . . ; p and dRdx j j j explicitly found (see the Appendix). A good and useful approximation to these derivatives can be found numerically. In all the methods described below we start at a randomly generated (or generated by another method) starting solution X 0 . Let Dxi and Dy i be the movement of generator point i. We solve these problems iteratively by evaluating a better solution X 1 ¼ X 0 þ ðDx; DyÞ that becomes the ‘‘new” X 0 for the next iteration. The analyses below deal with finding ðDx; DyÞ for a given X 0 . 3.1. Solving the equitable load (EL) problem One way to solve the problem is to minimize the objective: 2 p X 1 S j ðX Þ ; p j¼1
ð1Þ
by using a gradient search. The step size of the gradient search is determined by a golden section search along the gradient. However, the following two approaches proved to be more efficient. 3.1.1. Using linear programming The problem can be solved iteratively by solving a linear program each iteration. We wish to move the generator points by the shortest possible distances so as to bring the areas to equal values. Assuming that the derivatives of the areas do not change much, this can be done by the following linear program: ( ) p X ðjDxi j þ jDy i jÞ ; ð2Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ i¼1
dxi
dS j ðX 0 Þ Dxi þ Dy i dy i
¼
1 p
j ¼ 1; . . . ; p:
The solution to the linear program is used to determine the next iterate. Note that if we assign Dxi ¼ Dxþ i Dxi to allow þ for negative displacements, then jDxi j ¼ Dxi þ Dxi , and similarly for Dy i , and (2) is indeed a linear program. The linear program may not have a feasible solution. In fact, we encountered many instances with no feasible solution. Therefore, a corrected linear program using a large parameter M to penalize the deviations from equality in the constraints (2), we used M ¼ 1000, is as follows: ( ) p p X X ðjDxi j þ jDy i jÞ þ M jsj j ð3Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ i¼1
where sj is a slack variable.
dxi
J ¼1
Dxi þ
dS j ðX 0 Þ Dy i dy i
þ sj ¼
1 p
j ¼ 1; . . . ; p;
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
23
3.1.2. The gradient search Since the solution to the linear program may be too large (for example, the generator points may move outside the area), we suggest a gradient search to find the next iterate. The next iterate is defined as X 0 þ hDX where DX is the solution to the linear program (3) and h is a positive scalar. A golden section search on the segment 0 6 h 6 1 is performed and the best h that minimizes the value of the objective function (1) is used for the next iterate. 3.1.3. Using quadratic programming Rather than formulating a linear program each iteration, we solve a quadratic program each iteration. This quadratic program is solved explicitly without a need for a simplex-like algorithm used to solve the linear program (3). The quadratic programming formulation: ( ) p X 2 2 ð½Dxi þ ½Dy i Þ ; ð4Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ
dxi
i¼1
Dxi þ
dS j ðX 0 Þ Dy i dy i
¼
1 p
j ¼ 1; . . . ; p:
Since there may be no feasible solution to (4) we suggest to solve for a large M: ( !) p p X X 2 2 ½Dxi þ ½Dy i þ M s2j ; min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ
dxi
i¼1
j¼1
dS j ðX 0 Þ Dxi þ Dy i dy i
þ sj ¼
1 p
ð5Þ
j ¼ 1; . . . ; p:
3.1.4. Solving the quadratic program For this section we define a matrix A of m rows and n columns, a vector b of m rows, a variable vector x of n columns and a scalar M P 0. Problem (5) can be written as T
min fxT x þ MðAx bÞ ðAx bÞg:
ð6Þ T
T
T
T
Equating the derivative of (6) by x to 0: 2x þ 2MfA Ax A bg ¼ 0, which yields fI þ MA Agx ¼ MA b. Thus, 1 1 x ¼ fI þ MAT Ag1 MAT b ¼ I þ AT A AT b: M
ð7Þ
Note that for any real matrix A the matrix AT A is positive semi-definite (all its eigenvalues are nonnegative) and thus, the matrix M1 I þ AT A is positive definite (all eigenvalues are greater than M1 ) and thus its inverse does exist. We suggest the same approach of using h as described in Section 3.1.2 for the iterations based on linear programming. We propose to perform a golden section search and identify the h that minimizes the value of the objective function (1). 3.2. Solving the p-center problem The objective is to reduce the value of maxfRik ðX Þg. Suppose that xj is changed to xj þ Dxj and y j is changed to y j þ Dy j for j ¼ 1; . . . ; p. By first order approximation ( ) p X dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j ð8Þ Rik ðX Þ Rik ðX 0 Þ þ i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ: dxj dy j j¼1 Assuming that the structure of the Voronoi diagram remains unchanged by the move and accepting approximation (8), the solution to the p-center problem can be found by the following linear program: min
fLg;
Subject to : Rik ðX 0 Þ þ
p X j¼1
(
dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j dxj dy j
ð9Þ
) 6L
i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ: Once the solution to the linear program (9) is obtained, the next iteration is ðxj þ Dxj ; y j þ Dy j Þ for j ¼ 1; . . . ; p.
24
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Notes: 1. According to Theorem 3 in the Appendix, the number of constraints needed in the formulation (9) is equal to the number of vertices in the Voronoi diagram because all radii centered Pp at the same vertex lead to identical constraints. 2. We suggest to find the ‘‘gradient” of L by adding a constraint j¼1 ½jDxj j þ jDy j j 6 d for a small d. Then, a golden section search along this gradient as described in Section 3.1.2 determines the next iterate. However, the value of h is not bounded by 1. To find the upper limit for h we used h ¼ 1 and if the value of the objective function is lower than its value at X 0 we keep doubling the value of h until the value of the objective function is greater than its value at X 0 . This h defines the upper bound for the golden section search. 3.3. Solving the minimum equitable radius (MER) problem To find the minimum possible radius among all EL solutions, we can combine the two linear programs (9) and (2) into one. fLg;
min
p X
(
)
dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j 6 L i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ; dx dy j j j¼1 p X dS j ðX 0 Þ dS j ðX 0 Þ 1 Dxi þ Dy i ¼ j ¼ 1; . . . ; p: S j ðX 0 Þ þ dxi dy i p i¼1
Subject to : Rik ðX 0 Þ þ
ð10Þ
ð11Þ
To avoid getting solutions with large values of the variables, it is recommended to first find a solution close to being equitable, by either the gradient approach for minimizing (1), or by solving (2) or by solving (7), and only then start iterations based on the linear program (10). In any case, to get the gradient of L we propose to add the constraint Pp j¼1 ½jDxj j þ jDy j j 6 d for a small d to the linear program and change constraint (11) to p X dS j ðX 0 Þ i¼1
dxi
Dxi þ
dS j ðX 0 Þ Dy i dy i
¼0
j ¼ 1; . . . ; p;
ð12Þ
and perform the gradient search that is described in Section 3.2. However, when performing a golden section search along the gradient, the areas deviate from being equal especially the move is relatively large. Therefore, we perform the Pwhen p 2 golden section search on the function maxfRik ðX Þg þ 100 i¼1 ðS i ðX Þ 1p Þ and apply the EL algorithm on the final point to assure that at the end of each iteration the areas are indeed equal to one another. 4. Computational experiments Programs were coded in FORTRAN using double precision arithmetic and compiled by Intel Fortran 9.0. Subroutines by Sugihara and Iri (1994) were used. The programs were run on a desktop Pentium IV 2.8 GHz computer with 256 MB RAM. Starting were randomly generated in a unit square according to the following procedure. For each value of pffiffiffiffiffiffiffiffiffiffiffiffiffiffisolutions ffi p; L ¼ b p 0:5c þ 1 was calculated. If p 6 LðL 1Þ, a grid of L 1 by L rectangles was constructed. Otherwise, a grid of L by L rectangles (squares) was constructed. A rectangle (that has no generator point) was randomly selected, and a generator point was randomly generated near its middle in a smaller rectangle centered at the middle of the rectangle and whose sides are 40% of the side of the rectangle. At the end of the process, some of the rectangles may have no generator point. One hundred starting solutions were randomly generated for each value of p. 4.1. Testing EL algorithms We compared the linear programming and the quadratic programming approaches to solve the EL problems. One hundred randomly generated starting solutions were generated for 3 6 p 6 50 and the problems solved until each area is within 0.001% of 1p. Both algorithms are very efficient. All 3 6 p 6 50 problems were solved within 0.4 seconds per problem using the linear programming approach and within 0.05 seconds using the quadratic programming approach. The quadratic programming approach was about 5–10 times faster than the linear programming approach. The quadratic programming approach is more efficient mainly because it is explicitly solved and there is no need for an iterative procedure (the simplex method) to solve the linear program.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
25
4.2. Solving p-center problems In Table 2 the results of solving the p-center problems are reported. For each value of p we report the best found solution, the number of times it was found (out of 100 runs), and the average time it took to solve one problem. The p-center problem has many local optima and therefore for a larger value of p the best known solution was found very few times out of 100 runs. It is likely that in many of the problems the optimum was missed. However, we believe that the best found solution is close to the optimal value. 4.3. Solving the MER problem The computational results of solving the MER problem are depicted in Table 3. The performance of the algorithm is similar to that of solving the p-center problem. The algorithm is very efficient for the problems of these sizes. The Table 2 Results for solving the p-center problem p
Best known
Times found
Time (second)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50389 0.35355 0.32616 0.29873 0.27429 0.26030 0.23064 0.21823 0.21252 0.20228 0.19431 0.18551 0.17966 0.16943 0.16568 0.16064 0.15818 0.15225 0.14895 0.14369 0.14124 0.13830 0.13446 0.13223 0.13087 0.12761 0.12616 0.12204 0.12208 0.11841 0.11789 0.11498 0.11336 0.11168 0.11110 0.10856 0.10770 0.10702 0.10508 0.10400 0.10234 0.10126 0.10028 0.09922 0.09793 0.09639 0.09558 0.09451
100 100 100 100 20 61 65 27 12 28 13 11 5 6 4 4 4 8 1 2 1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.00 0.01 0.01 0.02 0.05 0.04 0.05 0.08 0.07 0.09 0.16 0.18 0.19 0.21 0.32 0.38 0.42 0.44 0.68 0.68 0.75 0.79 0.85 1.21 1.30 1.36 1.49 1.53 2.24 2.43 2.49 2.48 2.75 2.73 3.80 3.81 3.90 4.50 4.66 5.05 6.63 7.01 7.03 7.41 7.80 7.94 7.95 11.48
26
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Table 3 Results for solving the MER problem p
Best known
Times found
Time (second)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50714 0.35355 0.34874 0.30046 0.29243 0.26909 0.23570 0.23702 0.21917 0.20833 0.20480 0.19339 0.18623 0.17677 0.17780 0.16976 0.16417 0.15986 0.15653 0.15266 0.14853 0.14520 0.14142 0.14012 0.13764 0.13369 0.13190 0.12954 0.12736 0.12446 0.12311 0.12130 0.11953 0.11785 0.11800 0.11507 0.11392 0.11238 0.11079 0.10952 0.10894 0.10723 0.10515 0.10444 0.10317 0.10206 0.10102 0.10073
65 99 97 100 39 90 100 23 35 100 17 23 62 100 1 2 19 49 1 2 3 50 96 1 2 1 10 5 1 1 2 5 12 59 1 1 1 1 1 5 2 1 1 2 1 2 14 1
0.02 0.02 0.03 0.04 0.05 0.07 0.09 0.11 0.13 0.17 0.20 0.25 0.28 0.35 0.45 0.50 0.60 0.64 0.84 0.90 1.00 1.15 1.22 1.51 1.71 1.87 2.00 1.98 2.31 2.39 2.57 2.83 2.87 2.82 3.68 4.03 4.24 4.34 4.23 4.28 5.63 5.89 6.43 6.60 6.88 7.08 7.15 8.05
MER problem has many local minima and thus many of the solutions, especially for larger values of p, are likely not optimal but we believe that they are close to the optimal solution. 5. Conclusions We propose efficient approaches to solve location problems with continuous demand. More efficient solution procedures are proposed for the solution of the p-center problem in an area (Suzuki and Drezner, 1996) and the equitable load problem (Baron et al., 2007). A new problem, the minimum equitable radius (MER) problem, is proposed. The objective is to find the minimum possible radius (the objective of the p-center problem) among all solutions with equal loads. Suzuki and Drezner (1996) for the p-center problem, and Baron et al. (2007) for the equitable load problem both proposed iterative procedures based on Voronoi diagrams that change the location of one facility at a time. The innovation in
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
27
our approach is that at every iteration the locations of all facilities are changed simultaneously based on a solution of either a linear program or a quadratic program. Suzuki and Drezner (1996) had to stop the iterations prematurely because of the slow convergence and added a ‘‘finishing up” separate procedure to complete the solution procedure. Baron et al. (2007) also terminated the iterations prematurely (1.5% above the optimum on an average) because of the slow convergence of their algorithm. Our approach converged very quickly to a local optimum for the p-center problem and to the global optimum for the equitable load problem. The equitable load problem converged to within 0.001% of the optimum in about five iterations in a very short time. It solved the p ¼ 50 problem in about 0.05 seconds. Convergence to the optimum is very fast near the optimum so that a couple of extra iterations would have solved the problem to the accuracy of the double precision capability of the computer. A similar algorithm was constructed for the solution of the MER problem with comparable performance. 5.1. Extensions and future research It is possible to construct similar algorithms for other continuous demand problems such as the p-median, or competitive objectives and/or problems with non-uniform demand. This requires the derivation of the derivatives of such objectives by a change in the location of a facility. Such derivatives can be found either analytically, if possible, or numerically. Numerical calculation of the derivatives means to recalculate the Voronoi diagram by moving a generator point by a small distance and recalculating the value of the objective function following the move. The derivatives are approximated by dividing the difference between the values by the distance of the generator point’s move. Double precision arithmetic is required for such a calculation. Variants of the MER problem are of interest, for example, finding the best p-median objective among all EL solutions. Appendix. In this Appendix, we explicitly calculate the derivatives of the Voronoi areas and the Voronoi radii by a movement of a generator point. Calculating the derivative of the voronoi area Calculating the area Suppose that a Voronoi region (polygon) is defined by an ordered sequence of k vertices U 1 ; . . . ; U k where U i ¼ ðui ; vi Þ. Also define U kþ1 ¼ U 1 and U 0 ¼ U k . The area S of the polygon is (Beyer, 1981): S¼
k 1X ui ðviþ1 vi1 Þ: 2 i¼1
ð13Þ
Calculating the derivative It is difficult to evaluate the change in a vertex U j when a generator point X i is changed infinitesimally. We therefore had to resort to a formula which is not based on (13) and we found the following alternative approach. When the location of generator point i is changed infinitesimally, the edges of its Voronoi region change. This affects the area of the Voronoi region defined by generator point i and all adjacent Voronoi regions. Once we evaluate the change in the area as a result of the movement of one edge, the change in the Voronoi areas is the sum of the relevant individual changes. Note that although the endpoints of the edge change by an infinitesimal change in the location of a generator point, their impact on the change in the area is of a second order and thus such changes do not affect the derivative. Consider two generator points X 1 ¼ ðx1 ; y 1 Þ and X 2 ¼ ðx2 ; y 2 Þ that have a common Voronoi edge E whose end points are U 1 and U 2 . We find the change in the area due to the change of this edge as x1 is changed by an infinitesimal dx1 . We define dE . this derivative due to the change in edge E as dx 1 The line defining the edge is 2
2
2
2
ðx x1 Þ þ ðy y 1 Þ ¼ ðx x2 Þ þ ðy y 2 Þ ;
ð14Þ
or y¼
x2 x1 x2 þ y 22 x21 y 21 : xþ 2 y2 y1 2ðy 2 y 1 Þ
ð15Þ
The derivative of y for a given x by x1 is: dy x x1 ¼ : dx1 y 2 y 1
ð16Þ
28
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
The change in the area is: Z u2 o dE dy ðu2 u21 Þ=2 x1 ðu2 u1 Þ u2 u1 nu1 þ u2 x1 : ¼ dx ¼ 2 ¼ y2 y1 2 dx1 y2 y1 u1 dx1
ð17Þ
Note that since the points U 1 and U 2 are on line (15), substituting them as ðx; yÞ in (15) leads to the relationship u2 u1 v2 v1 ¼ : ð18Þ y2 y1 x2 x1 Relationship (18) can be used in (17) in the case y 1 ¼ y 2 . Also, simple manipulations using (18) lead to: 2 2 u2 u1 v2 v1 ðu2 u1 Þ2 þ ðv2 v1 Þ2 ¼ ¼ 2 2 y2 y1 x2 x1 ðx2 x1 Þ þ ðy 2 y 1 Þ 1 which means that uy 2 u in (17) is the ratio of the edge’s length and the distance between the two generator points (it may be y 2
1
negative). To apply this relationship one only needs to determine the sign of the ratio. This approach has the advantage of not getting an expression of 0/0 unless the two generator points coincide.
Calculating the derivative of the radius Consider Figs. 2–5. There are three case: 1. The distance is measured to an internal vertex. 2. The distance is measured to a vertex on a side of the square, and 3. The distance is measured to the corner of the square.
Case 1: The distance is measured to an internal vertex Every internal vertex of the Voronoi diagram is at the center of the circle of radius R circumscribing a triangle with three generator points as its vertices. Let the three generator points be ðxi ; y i Þ; i ¼ 1; 2; 3 and we find the derivative of R by x1 . Let the three sides of the triangle be a; b; c, where a is the side opposite ðx1 ; y 1 Þ and a be the angle opposite side a. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ ðx2 x3 Þ2 þ ðy 2 y 3 Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 b ¼ ðx1 x3 Þ þ ðy 1 y 3 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 c ¼ ðx1 x2 Þ þ ðy 1 y 2 Þ : a The circumscribing radius is R ¼ 2 sin . Note that this formula is true also for a > p=2. However, in this case the center of a the circle is outside the triangle but the circle of radius R passes through the three generator points. The circumscribing circle in this case is actually at the center of the longest side of the triangle, and the third generator point is inside that circle. In this context we still call R as the radius of the circumscribing circle because it is equidistant from all three generator points. When the first generator point is moved, a may change, but a (the opposite side) remains unchanged. Therefore,
dR a cos a R ¼ ¼ : da tan a 2 sin2 a We now find
da . dx1
ð19Þ
By direct calculation:
db x1 x3 ¼ b dx1 dc x1 x2 : ¼ c dx1 By the cosine theorem: 2bc cos a ¼ b2 þ c2 a2 ;
ð20Þ
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
29
differentiating (20) by x1 yields: da dc db db dc ¼ 2b þ 2 cos a b þc þ 2c 2bc sin a dx1 dx1 dx1 dx1 dx1
da b c þ cos a ðx1 x2 Þ þ ðx1 x3 Þ ¼ 2ð2x1 x2 x3 Þ bc sin a dx1 c b h i x1 x3 2x1 x2 x3 x1 x2 da cos a c2 þ b2 bc : ¼ ¼ 2x1 x2 x3 sin a dx1
Combining (19) and (21) yields: dR dR da R ¼ ¼ dx1 da dx1 tan a
cos a
h
x1 x2 c2
i 3 þ x1bx 2x1 xbc2 x3 2 sin a
¼
ð21Þ
a cos a 2x1 x2 x3 x1 x2 x1 x3 ; cos a þ bc c2 b2 2 sin3 a
ð22Þ
and similarly
dR a cos a 2y 1 y 2 y 3 y1 y2 y1 y3 ; cos a ¼ þ bc c2 dy 1 2 sin3 a b2
ð23Þ
where cos a ¼
b2 þ c 2 a2 : 2bc
The following Theorem can save 23 of the effort required to calculate all the derivatives. Theorem 3. Let a vertex k be defined by three generator points i1 ; i2 ; i3 . Then, same is true for the derivatives by i2 and i3 .
dRi1 k dxi1
¼
dRi2 k dxi1
¼
dRi3 k dxi1
and
dRi1 k dy i1
¼
dRi2 k dy i1
¼
dRi3 k dy i1 .
The
Proof. Vertex k is the center of the circumscribing circle of the triangle whose vertices are i1 ; i2 , and i3 . Therefore, Ri1 k ¼ Ri2 k ¼ Ri3 k . Following a small change in the location of i1 , vertex k moves to the new center of the triangle and all three distances are still equal to one another. Therefore, the three derivatives are equal to one another. h Case 2: The distance is measured to a vertex on a side of the square Such a vertex is equally distant from two generator points ðx1 ; y 1 Þ and ðx2 ; y 2 Þ. As in Case 1 we find the derivative of the distance by x1 . Consider ðx3 ; y 3 Þ as the mirror image of the second generator point by the side of the square. This mirror image is outside the square. The side of the square is the perpendicular bisector between the second and third points. Therefore, the vertex is equi-distant from all three points. When x1 is moved and the other points remain unchanged, the vertex moves along the side of the square and Eqs. (22) and (23) hold. Case 3: The distance is measured to a corner of the square In this case the derivative is simply the derivative of the distance to a fixed point. Let the corner of the square be at qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dR 1 ðc1 ; c2 Þ. The distance is R ¼ ðx1 c1 Þ þ ðy 1 c2 Þ and the derivative is dx ¼ x1 c . R 1 References Baron, O., Berman, O., Krass, D., Wang, Q., 2007. The equitable location problem on the plane. European Journal of Operational Research 183, 578–590. Berman O., Drezner, Z., Tamir, A., Wesolowsky, G.O., in press. Optimal Location with Equitable Loads, Annals of Operations Research. Beyer, W.H., 1981. CRC Standard Mathematical Tables, 26th ed. CRC Press, Boca Raton, FL. Current, J., Daskin, M., Schilling, D., 2002. Discrete network location models, In: Drezner, Z., Hamacher, H.W. (Eds.), Facility Location: Applications and Theory, pp. 81–118. Daskin, M.S., 1995. Network and Discrete Location: Models, Algorithms, and Applications. John Wiley & Sons, New York. Drezner, Z., 1984. The p-center problem – heuristic and optimal algorithms. Journal of the Operational Research Society 35, 741–748. Drezner, Z., 1987. On the rectangular p-center problem. Naval Research Logistics Quarterly 34, 229–234. Drezner, T., Drezner, Z., 2006. Multiple facilities location in the plane using the gravity model. Geographical Analysis 38, 391–406. Handler, G.Y., 1990. p-Center problems. In: Mirchandani, P.B., Francis, R.L. (Eds.), Discrete Location Theory. Wiley Inter-Science, New York, pp. 305– 347. Ohya, T., Iri, M., Murota, K., 1984. Improvements of the incremental method of the Voronoi diagram with computational comparison of various algorithms. Journal of the Operations Research Society of Japan 27, 306–337.
30
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Okabe, A., Boots, B., Sugihara, K., Chiu, S.N., 2000. Spatial tessellations: Concepts and applications of Voronoi diagrams. Wiley Series in Probability and Statistics. Okabe, A., Suzuki, A., 1987. Stability of spatial competition for a large number of firms on a bounded two-dimensional space. Environment and Planning A 16, 107–114. Sugihara, K., Iri, M., 1994. A robust topology-oriented incremental algorithm for Voronoi diagram. International Journal of Computational Geometry and Applications 4, 179–228. Suzuki, A., Drezner, Z., 1996. The p-center location problem in an area. Location Science 4, 69–82. Suzuki, A., Okabe, A., 1995. Using Voronoi diagrams. In: Drezner, Z. (Ed.), Facility Location: A Survey of Applications and Methods. Springer, New York. Vijay, J., 1985. An algorithm for the p-center problem in the plane. Transportation Science 19, 235–245.
European Journal of Operational Research 195 (2009) 17–30 www.elsevier.com/locate/ejor
Continuous Optimization
The minimum equitable radius location problem with continuous demand Atsuo Suzuki a, Zvi Drezner b,* b
a Department of Mathematical Sciences, Nanzan University, 27 Seirei, Seto, Aichi 489-0863, Japan College of Business and Economics, California State University-Fullerton, Fullerton, CA 92834, United States
Received 6 September 2007; accepted 10 January 2008 Available online 26 January 2008
Abstract We analyze the location of p facilities satisfying continuous area demand. Three objectives are considered: (i) the p-center objective (to minimize the maximum distance between all points in the area and their closest facility), (ii) equalizing the load service by the facilities, and (iii) the minimum equitable radius – minimizing the maximum radius from each point to its closest facility subject to the constraint that each facility services the same load. The paper offers three contributions: (i) a new problem – the minimum equitable radius is presented and solved by an efficient algorithm, (ii) an improved and efficient algorithm is developed for the solution of the p-center problem, and (iii) an improved algorithm for the equitable load problem is developed. Extensive computational experiments demonstrated the superiority of the new solution algorithms. Ó 2008 Elsevier B.V. All rights reserved. Keywords: Location; p-Center; Equitable loads; Voronoi diagrams
1. Introduction Recently, problems of locating p facilities to satisfy continuous demand in an area were investigated. Ohya et al. (1984) formulated and solved the p-median problem in an area by applying an iterative procedure on the Voronoi diagram of the facilities (Okabe et al., 2000; Suzuki and Okabe, 1995). A similar approach was utilized by Okabe and Suzuki (1987) to formulate and solve the competitive location problem in an area and by Suzuki and Drezner (1996) to formulate and solve the p-center problem in an area. Baron et al. (2007) suggested an equitable load model (EL) for uniform continuous demand in an area. Their model seeks the location of p facilities in the area, each point in the area is serviced by the closest facility as long as it is within a given coverage distance r from the facility. The objective is equalizing the loads serviced by each facility which is formulated as minimizing the maximum load serviced by the facilities. If the value of r is too small, there may be no feasible solution. For a large enough value of r, the optimal solution is equal in areas for all facilities. They solved the EL problem by an iterative procedure based on Voronoi diagrams. The question of finding the minimum possible r among all optimal solutions to the equitable load problem, is a practical and interesting problem. We term this problem ‘‘the minimum equitable radius” (MER) problem. The EL solution is usually not unique. Therefore, the maximum distance to the nearest facility may vary from one optimal solution to another. The MER problem is to find the minimum value among those maximum distances. To illustrate the MER problem consider the location of four facilities in a square. If the four facilities are placed equally spaced on the periphery of a circle centered at the center of the square (the radius of the circle is irrelevant), they divide the *
Corresponding author. Tel.: +1 714 278 2712; fax: +1 714 278 5940. E-mail address: [email protected] (Z. Drezner).
0377-2217/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2008.01.022
18
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
square into four equal areas (see Fig. 1). This is because the perpendicular bisector between two facilities must pass through the center of the square, and the four areas are thus congruent. Each angle h and any radius lead to the locations of four facilities which are optimal solutions to the EL problem. Therefore, the number of EL solutions is infinite. It is easy to pffiffi show that the radius enclosing each shape is 2 cos 2ph . This radius is obtained when the facilities are located at the centers ð4 Þ of the lines connecting the intersection points of the perpendicular bisectors with the sides of the square. The optimal solution of the MER problem (the minimum possible radius), as proven in this paper, is obtained for h ¼ p4 when the facilities are located on the diagonals at the centers of the lines connecting the center of the square and its vertices. The square is divided into four equal small squares. The p-center problem is investigated in other environments as well. In the plane, it is investigated by Drezner (1984) and Vijay (1985) for Euclidean distances and by Drezner (1987) for rectilinear distances. In a network environment the p-center problem is extensively investigated (for example, Handler, 1990), see, Daskin (1995) and Current et al. (2002) for a review. The EL problem with no coverage radius constraints was investigated by Drezner and Drezner (2006) for planar problems and by Berman et al. (in press) in the network environment. We are not aware of any paper suggesting the MER problem. There are many situations in which equal loads of the facilities is desirable. When cell phone towers are established in an area, the busiest tower determines the efficiency of the system thus equalizing loads are desirable. This applies to any design of p identical M/M/1 facilities where minimizing the maximum waiting time at the facilities is equivalent to minimizing the maximum arrival rate at the facilities which is proportional to the number of customers getting service at the facility. When carving market territories to be assigned to marketing representatives, the objective is to assign to each territory an equitable market potential. Another example is the problem of designing machines with similar capacities in a production facility. The throughput of the system depends on the machine with the maximum load. Voter districting is another application. The objective is that each district will have about the same number of voters. The maximum distance to the closest facility is a measure of good service used in numerous p-center applications. The MER model combines the two objectives. We seek the minimum possible radius among all optimal solutions to the EL model thus achieving both the EL and the p-center objectives. In every application described above for the EL model, it is desirable to minimize the maximum distance to the closest facility in addition to equalizing the loads as much as possible. For example, it is desirable that voters in a district are close to the voting place. It is advantageous for customers to be close to the facilities. In this paper, we propose three contributions. First, we propose a better heuristic algorithm to solve the p-center problem in an area. Second, pwe ffiffiffi propose a better heuristic algorithm for the solution of the EL problem for a large enough value of r (for example r P 2 for a square area as suggested by Baron et al., 2007) so that the distance constraints are not active in the solution. Third, we propose and heuristically solve the MER problem of finding the smallest possible radius among all equitable loads solutions. The problems are defined for any demand area and any demand distribution. However, most of the analysis in this paper pertains to uniform demand and a square area. It can be generalized to any convex polygonal shape and non-uniform demand. Finding the Voronoi diagram inside a polygonal area is quite easy. If demand is not uniform, calculating the
θ
Fig. 1. EL and MER solutions.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
19
necessary derivatives may require numerical differentiation. Note that the distribution of demand is irrelevant for the p-center problem as long as demand’s density function is greater than zero at every point in the area. 2. Discussion A Voronoi diagram defined by a set of generator points (for example, facilities) is a partition of the area into Voronoi regions so that all the points in a region are closest to a generator point, see, Figs. 2–5. Each generator point defines its associated Voronoi area which is a polygon surrounding it. The sides of the polygon are perpendicular bisectors between generator points. If demand is uniform, then the load of each facility is the area of the Voronoi region associated with the facility. Baron et al. (2007) showed that, in a square area, for every p (if r is large enough so that every demand point is covered) there exists an optimal solution by creating p ‘‘thin” rectangles, whose length is ‘‘1”, with the same width of 1p side by side
MER Solution
p-Center Solution Fig. 2. p ¼ 3 solutions.
MER Solution
p-Center Solution Fig. 3. p ¼ 5 solutions.
20
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
MER Solution
p-Center Solution Fig. 4. p ¼ 7 solutions.
MER Solution
p-Center Solution Fig. 5. p ¼ 8 solutions.
and placing a facility at the center of each rectangle. This configuration yields p equal loads and thus must be optimal. Since we know that there exists an optimal solution for which all loads are equal, the MER problem must have a solution as well. We conjecture that for every p, when r is large enough, there is a continuum of optimal solutions. If our conjecture is true, then finding the smallest possible r among the continuum of optimal solutions (the MER problem) can be a challenging one. The following two theorems provide insight into the MER problem and its relationship to the p-center problem: Theorem 1. The minimum radius for the p-center solution is a lower bound for the MER solution. Proof. Every feasible solution for the MER problem is a feasible solution to the p-center problem. Therefore, the optimal solution to the MER problem cannot be better than the optimal solution to the p-center problem. h In the rest of the paper we assume that the MER problem is defined in a square of side 1.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Theorem 2. The optimal solution to the MER problem is lower than or equal to Proof. The radius enclosing the ‘‘thin” rectangles described above is
0:5 p
0:5 p
21
pffiffiffiffiffiffiffiffiffiffiffiffiffi p 2 þ 1.
pffiffiffiffiffiffiffiffiffiffiffiffiffi p 2 þ 1. h
Note that by Theorem 1 the solution to the p ¼ 4 problem is to split the square into four small squares because this is the optimal solution to the p-center problem (Suzuki and Drezner, 1996) and it is feasible for the MER problem. However, the optimal solution to the p-center problem for p ¼ 9 is not dividing the square into 9 small squares, which yields a radius of 0.23570. A solution with r ¼ 0:23064 is given by Suzuki and Drezner (1996). Therefore, the nine small squares solution may not be the optimal solution to the MER problem. The MER solution and the p-center solution solve the EL problem (Baron et al., 2007) for most values of r. If the optimal radius of the MER problem does not exceed the given r, it is a solution to the EL problem. If the optimal radius to the MER problem exceeds the given r, then there is no solution to the EL problem with equal loads. If r is less than the p-center solution, then there is no feasible solution to the EL problem. The EL problem remains of interest only for r between the Table 1 Comparing p-center and MER solutions p
MER
p-center
Diff.
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50714 0.35355 0.34874 0.30046 0.29243 0.26909 0.23570 0.23702 0.21917 0.20833 0.20480 0.19339 0.18623 0.17677 0.17780 0.16976 0.16417 0.15986 0.15653 0.15266 0.14853 0.14520 0.14142 0.14012 0.13764 0.13369 0.13190 0.12954 0.12736 0.12446 0.12311 0.12130 0.11953 0.11785 0.11800 0.11507 0.11392 0.11238 0.11079 0.10952 0.10894 0.10723 0.10515 0.10444 0.10317 0.10206 0.10102 0.10073
0.50389 0.35355 0.32616 0.29873 0.27429 0.26030 0.23064 0.21823 0.21252 0.20228 0.19431 0.18551 0.17966 0.16943 0.16568 0.16064 0.15818 0.15225 0.14895 0.14369 0.14124 0.13830 0.13446 0.13223 0.13087 0.12761 0.12616 0.12204 0.12208 0.11841 0.11789 0.11498 0.11336 0.11168 0.11110 0.10856 0.10770 0.10702 0.10508 0.10400 0.10234 0.10126 0.10028 0.09922 0.09793 0.09639 0.09558 0.09451
0.00325 0.00000 0.02258 0.00173 0.01814 0.00879 0.00506 0.01879 0.00666 0.00605 0.01049 0.00788 0.00657 0.00735 0.01212 0.00912 0.00600 0.00761 0.00758 0.00896 0.00729 0.00690 0.00696 0.00788 0.00677 0.00608 0.00574 0.00750 0.00528 0.00604 0.00522 0.00632 0.00617 0.00618 0.00690 0.00651 0.00622 0.00536 0.00571 0.00552 0.00661 0.00597 0.00487 0.00522 0.00524 0.00567 0.00545 0.00622
22
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
p-center solution and the MER solution for which there is no solution to the EL problem with equal areas. In Table 1, the solutions found in the computational experiments section below for the MER and p-center problems are summarized. For pffiffiffiffi example, assuming that the values in Table 1 are optimal, for p ¼ 3: if r < 0:50389 ð¼ 1665 which is optimal by Suzuki and Drezner, 1996), there is no feasible solution to the EL problem, and if r P 0:50714, the MER solution is optimal for the EL problem. In the range 0:50389 6 r < 0:50714, a range of only 0.00325, the solution to the EL problem is not with equal areas and the algorithms in Baron et al. (2007) need to be employed to find such a solution. These two solutions are depicted in Fig. 2. The p-center solution is the minimum possible radius, and the MER solution is the minimum possible radius among all partitions of the square into three equal areas. Note that the ranges for requiring such a solution are quite small for all values of p depicted in Table 1. 3. Solution algorithms A vector X of p generator points is defined as X ¼ fðxi ; y i Þg; i ¼ 1; . . . ; p. Each point defines the closest Voronoi region V i ðX Þ (see Figs. 2–5). Such a Voronoi diagram defines a set of KðX Þ vertices. Let the set P i ðX Þ be the set of all Voronoi vertices of the Voronoi region defined by generator point i. Let also S i ðX Þ be the area of the Voronoi region (polygon) associated with generator point i and Rik ðX Þ for i ¼ 1; . . . ; p and k 2 P i ðX Þ be the Euclidean distance between generator point i and vertex k 2 P i ðX Þ belong to its Voronoi region. ik ðX Þ dRik ðX Þ ; dy for i ¼ 1; . . . ; p; k 2 P i ðX Þ; j ¼ 1; . . . ; p can be The derivatives dSdxi ðXj Þ ; dSdyi ðX Þ for i ¼ 1; . . . ; p; j ¼ 1; . . . ; p and dRdx j j j explicitly found (see the Appendix). A good and useful approximation to these derivatives can be found numerically. In all the methods described below we start at a randomly generated (or generated by another method) starting solution X 0 . Let Dxi and Dy i be the movement of generator point i. We solve these problems iteratively by evaluating a better solution X 1 ¼ X 0 þ ðDx; DyÞ that becomes the ‘‘new” X 0 for the next iteration. The analyses below deal with finding ðDx; DyÞ for a given X 0 . 3.1. Solving the equitable load (EL) problem One way to solve the problem is to minimize the objective: 2 p X 1 S j ðX Þ ; p j¼1
ð1Þ
by using a gradient search. The step size of the gradient search is determined by a golden section search along the gradient. However, the following two approaches proved to be more efficient. 3.1.1. Using linear programming The problem can be solved iteratively by solving a linear program each iteration. We wish to move the generator points by the shortest possible distances so as to bring the areas to equal values. Assuming that the derivatives of the areas do not change much, this can be done by the following linear program: ( ) p X ðjDxi j þ jDy i jÞ ; ð2Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ i¼1
dxi
dS j ðX 0 Þ Dxi þ Dy i dy i
¼
1 p
j ¼ 1; . . . ; p:
The solution to the linear program is used to determine the next iterate. Note that if we assign Dxi ¼ Dxþ i Dxi to allow þ for negative displacements, then jDxi j ¼ Dxi þ Dxi , and similarly for Dy i , and (2) is indeed a linear program. The linear program may not have a feasible solution. In fact, we encountered many instances with no feasible solution. Therefore, a corrected linear program using a large parameter M to penalize the deviations from equality in the constraints (2), we used M ¼ 1000, is as follows: ( ) p p X X ðjDxi j þ jDy i jÞ þ M jsj j ð3Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ i¼1
where sj is a slack variable.
dxi
J ¼1
Dxi þ
dS j ðX 0 Þ Dy i dy i
þ sj ¼
1 p
j ¼ 1; . . . ; p;
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
23
3.1.2. The gradient search Since the solution to the linear program may be too large (for example, the generator points may move outside the area), we suggest a gradient search to find the next iterate. The next iterate is defined as X 0 þ hDX where DX is the solution to the linear program (3) and h is a positive scalar. A golden section search on the segment 0 6 h 6 1 is performed and the best h that minimizes the value of the objective function (1) is used for the next iterate. 3.1.3. Using quadratic programming Rather than formulating a linear program each iteration, we solve a quadratic program each iteration. This quadratic program is solved explicitly without a need for a simplex-like algorithm used to solve the linear program (3). The quadratic programming formulation: ( ) p X 2 2 ð½Dxi þ ½Dy i Þ ; ð4Þ min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ
dxi
i¼1
Dxi þ
dS j ðX 0 Þ Dy i dy i
¼
1 p
j ¼ 1; . . . ; p:
Since there may be no feasible solution to (4) we suggest to solve for a large M: ( !) p p X X 2 2 ½Dxi þ ½Dy i þ M s2j ; min i¼1
Subject to : S j ðX 0 Þ þ
p X dS j ðX 0 Þ
dxi
i¼1
j¼1
dS j ðX 0 Þ Dxi þ Dy i dy i
þ sj ¼
1 p
ð5Þ
j ¼ 1; . . . ; p:
3.1.4. Solving the quadratic program For this section we define a matrix A of m rows and n columns, a vector b of m rows, a variable vector x of n columns and a scalar M P 0. Problem (5) can be written as T
min fxT x þ MðAx bÞ ðAx bÞg:
ð6Þ T
T
T
T
Equating the derivative of (6) by x to 0: 2x þ 2MfA Ax A bg ¼ 0, which yields fI þ MA Agx ¼ MA b. Thus, 1 1 x ¼ fI þ MAT Ag1 MAT b ¼ I þ AT A AT b: M
ð7Þ
Note that for any real matrix A the matrix AT A is positive semi-definite (all its eigenvalues are nonnegative) and thus, the matrix M1 I þ AT A is positive definite (all eigenvalues are greater than M1 ) and thus its inverse does exist. We suggest the same approach of using h as described in Section 3.1.2 for the iterations based on linear programming. We propose to perform a golden section search and identify the h that minimizes the value of the objective function (1). 3.2. Solving the p-center problem The objective is to reduce the value of maxfRik ðX Þg. Suppose that xj is changed to xj þ Dxj and y j is changed to y j þ Dy j for j ¼ 1; . . . ; p. By first order approximation ( ) p X dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j ð8Þ Rik ðX Þ Rik ðX 0 Þ þ i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ: dxj dy j j¼1 Assuming that the structure of the Voronoi diagram remains unchanged by the move and accepting approximation (8), the solution to the p-center problem can be found by the following linear program: min
fLg;
Subject to : Rik ðX 0 Þ þ
p X j¼1
(
dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j dxj dy j
ð9Þ
) 6L
i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ: Once the solution to the linear program (9) is obtained, the next iteration is ðxj þ Dxj ; y j þ Dy j Þ for j ¼ 1; . . . ; p.
24
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Notes: 1. According to Theorem 3 in the Appendix, the number of constraints needed in the formulation (9) is equal to the number of vertices in the Voronoi diagram because all radii centered Pp at the same vertex lead to identical constraints. 2. We suggest to find the ‘‘gradient” of L by adding a constraint j¼1 ½jDxj j þ jDy j j 6 d for a small d. Then, a golden section search along this gradient as described in Section 3.1.2 determines the next iterate. However, the value of h is not bounded by 1. To find the upper limit for h we used h ¼ 1 and if the value of the objective function is lower than its value at X 0 we keep doubling the value of h until the value of the objective function is greater than its value at X 0 . This h defines the upper bound for the golden section search. 3.3. Solving the minimum equitable radius (MER) problem To find the minimum possible radius among all EL solutions, we can combine the two linear programs (9) and (2) into one. fLg;
min
p X
(
)
dRik ðX 0 Þ dRik ðX 0 Þ Dxj þ Dy j 6 L i ¼ 1; . . . ; p; k 2 P i ðX 0 Þ; dx dy j j j¼1 p X dS j ðX 0 Þ dS j ðX 0 Þ 1 Dxi þ Dy i ¼ j ¼ 1; . . . ; p: S j ðX 0 Þ þ dxi dy i p i¼1
Subject to : Rik ðX 0 Þ þ
ð10Þ
ð11Þ
To avoid getting solutions with large values of the variables, it is recommended to first find a solution close to being equitable, by either the gradient approach for minimizing (1), or by solving (2) or by solving (7), and only then start iterations based on the linear program (10). In any case, to get the gradient of L we propose to add the constraint Pp j¼1 ½jDxj j þ jDy j j 6 d for a small d to the linear program and change constraint (11) to p X dS j ðX 0 Þ i¼1
dxi
Dxi þ
dS j ðX 0 Þ Dy i dy i
¼0
j ¼ 1; . . . ; p;
ð12Þ
and perform the gradient search that is described in Section 3.2. However, when performing a golden section search along the gradient, the areas deviate from being equal especially the move is relatively large. Therefore, we perform the Pwhen p 2 golden section search on the function maxfRik ðX Þg þ 100 i¼1 ðS i ðX Þ 1p Þ and apply the EL algorithm on the final point to assure that at the end of each iteration the areas are indeed equal to one another. 4. Computational experiments Programs were coded in FORTRAN using double precision arithmetic and compiled by Intel Fortran 9.0. Subroutines by Sugihara and Iri (1994) were used. The programs were run on a desktop Pentium IV 2.8 GHz computer with 256 MB RAM. Starting were randomly generated in a unit square according to the following procedure. For each value of pffiffiffiffiffiffiffiffiffiffiffiffiffiffisolutions ffi p; L ¼ b p 0:5c þ 1 was calculated. If p 6 LðL 1Þ, a grid of L 1 by L rectangles was constructed. Otherwise, a grid of L by L rectangles (squares) was constructed. A rectangle (that has no generator point) was randomly selected, and a generator point was randomly generated near its middle in a smaller rectangle centered at the middle of the rectangle and whose sides are 40% of the side of the rectangle. At the end of the process, some of the rectangles may have no generator point. One hundred starting solutions were randomly generated for each value of p. 4.1. Testing EL algorithms We compared the linear programming and the quadratic programming approaches to solve the EL problems. One hundred randomly generated starting solutions were generated for 3 6 p 6 50 and the problems solved until each area is within 0.001% of 1p. Both algorithms are very efficient. All 3 6 p 6 50 problems were solved within 0.4 seconds per problem using the linear programming approach and within 0.05 seconds using the quadratic programming approach. The quadratic programming approach was about 5–10 times faster than the linear programming approach. The quadratic programming approach is more efficient mainly because it is explicitly solved and there is no need for an iterative procedure (the simplex method) to solve the linear program.
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
25
4.2. Solving p-center problems In Table 2 the results of solving the p-center problems are reported. For each value of p we report the best found solution, the number of times it was found (out of 100 runs), and the average time it took to solve one problem. The p-center problem has many local optima and therefore for a larger value of p the best known solution was found very few times out of 100 runs. It is likely that in many of the problems the optimum was missed. However, we believe that the best found solution is close to the optimal value. 4.3. Solving the MER problem The computational results of solving the MER problem are depicted in Table 3. The performance of the algorithm is similar to that of solving the p-center problem. The algorithm is very efficient for the problems of these sizes. The Table 2 Results for solving the p-center problem p
Best known
Times found
Time (second)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50389 0.35355 0.32616 0.29873 0.27429 0.26030 0.23064 0.21823 0.21252 0.20228 0.19431 0.18551 0.17966 0.16943 0.16568 0.16064 0.15818 0.15225 0.14895 0.14369 0.14124 0.13830 0.13446 0.13223 0.13087 0.12761 0.12616 0.12204 0.12208 0.11841 0.11789 0.11498 0.11336 0.11168 0.11110 0.10856 0.10770 0.10702 0.10508 0.10400 0.10234 0.10126 0.10028 0.09922 0.09793 0.09639 0.09558 0.09451
100 100 100 100 20 61 65 27 12 28 13 11 5 6 4 4 4 8 1 2 1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.00 0.01 0.01 0.02 0.05 0.04 0.05 0.08 0.07 0.09 0.16 0.18 0.19 0.21 0.32 0.38 0.42 0.44 0.68 0.68 0.75 0.79 0.85 1.21 1.30 1.36 1.49 1.53 2.24 2.43 2.49 2.48 2.75 2.73 3.80 3.81 3.90 4.50 4.66 5.05 6.63 7.01 7.03 7.41 7.80 7.94 7.95 11.48
26
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Table 3 Results for solving the MER problem p
Best known
Times found
Time (second)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
0.50714 0.35355 0.34874 0.30046 0.29243 0.26909 0.23570 0.23702 0.21917 0.20833 0.20480 0.19339 0.18623 0.17677 0.17780 0.16976 0.16417 0.15986 0.15653 0.15266 0.14853 0.14520 0.14142 0.14012 0.13764 0.13369 0.13190 0.12954 0.12736 0.12446 0.12311 0.12130 0.11953 0.11785 0.11800 0.11507 0.11392 0.11238 0.11079 0.10952 0.10894 0.10723 0.10515 0.10444 0.10317 0.10206 0.10102 0.10073
65 99 97 100 39 90 100 23 35 100 17 23 62 100 1 2 19 49 1 2 3 50 96 1 2 1 10 5 1 1 2 5 12 59 1 1 1 1 1 5 2 1 1 2 1 2 14 1
0.02 0.02 0.03 0.04 0.05 0.07 0.09 0.11 0.13 0.17 0.20 0.25 0.28 0.35 0.45 0.50 0.60 0.64 0.84 0.90 1.00 1.15 1.22 1.51 1.71 1.87 2.00 1.98 2.31 2.39 2.57 2.83 2.87 2.82 3.68 4.03 4.24 4.34 4.23 4.28 5.63 5.89 6.43 6.60 6.88 7.08 7.15 8.05
MER problem has many local minima and thus many of the solutions, especially for larger values of p, are likely not optimal but we believe that they are close to the optimal solution. 5. Conclusions We propose efficient approaches to solve location problems with continuous demand. More efficient solution procedures are proposed for the solution of the p-center problem in an area (Suzuki and Drezner, 1996) and the equitable load problem (Baron et al., 2007). A new problem, the minimum equitable radius (MER) problem, is proposed. The objective is to find the minimum possible radius (the objective of the p-center problem) among all solutions with equal loads. Suzuki and Drezner (1996) for the p-center problem, and Baron et al. (2007) for the equitable load problem both proposed iterative procedures based on Voronoi diagrams that change the location of one facility at a time. The innovation in
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
27
our approach is that at every iteration the locations of all facilities are changed simultaneously based on a solution of either a linear program or a quadratic program. Suzuki and Drezner (1996) had to stop the iterations prematurely because of the slow convergence and added a ‘‘finishing up” separate procedure to complete the solution procedure. Baron et al. (2007) also terminated the iterations prematurely (1.5% above the optimum on an average) because of the slow convergence of their algorithm. Our approach converged very quickly to a local optimum for the p-center problem and to the global optimum for the equitable load problem. The equitable load problem converged to within 0.001% of the optimum in about five iterations in a very short time. It solved the p ¼ 50 problem in about 0.05 seconds. Convergence to the optimum is very fast near the optimum so that a couple of extra iterations would have solved the problem to the accuracy of the double precision capability of the computer. A similar algorithm was constructed for the solution of the MER problem with comparable performance. 5.1. Extensions and future research It is possible to construct similar algorithms for other continuous demand problems such as the p-median, or competitive objectives and/or problems with non-uniform demand. This requires the derivation of the derivatives of such objectives by a change in the location of a facility. Such derivatives can be found either analytically, if possible, or numerically. Numerical calculation of the derivatives means to recalculate the Voronoi diagram by moving a generator point by a small distance and recalculating the value of the objective function following the move. The derivatives are approximated by dividing the difference between the values by the distance of the generator point’s move. Double precision arithmetic is required for such a calculation. Variants of the MER problem are of interest, for example, finding the best p-median objective among all EL solutions. Appendix. In this Appendix, we explicitly calculate the derivatives of the Voronoi areas and the Voronoi radii by a movement of a generator point. Calculating the derivative of the voronoi area Calculating the area Suppose that a Voronoi region (polygon) is defined by an ordered sequence of k vertices U 1 ; . . . ; U k where U i ¼ ðui ; vi Þ. Also define U kþ1 ¼ U 1 and U 0 ¼ U k . The area S of the polygon is (Beyer, 1981): S¼
k 1X ui ðviþ1 vi1 Þ: 2 i¼1
ð13Þ
Calculating the derivative It is difficult to evaluate the change in a vertex U j when a generator point X i is changed infinitesimally. We therefore had to resort to a formula which is not based on (13) and we found the following alternative approach. When the location of generator point i is changed infinitesimally, the edges of its Voronoi region change. This affects the area of the Voronoi region defined by generator point i and all adjacent Voronoi regions. Once we evaluate the change in the area as a result of the movement of one edge, the change in the Voronoi areas is the sum of the relevant individual changes. Note that although the endpoints of the edge change by an infinitesimal change in the location of a generator point, their impact on the change in the area is of a second order and thus such changes do not affect the derivative. Consider two generator points X 1 ¼ ðx1 ; y 1 Þ and X 2 ¼ ðx2 ; y 2 Þ that have a common Voronoi edge E whose end points are U 1 and U 2 . We find the change in the area due to the change of this edge as x1 is changed by an infinitesimal dx1 . We define dE . this derivative due to the change in edge E as dx 1 The line defining the edge is 2
2
2
2
ðx x1 Þ þ ðy y 1 Þ ¼ ðx x2 Þ þ ðy y 2 Þ ;
ð14Þ
or y¼
x2 x1 x2 þ y 22 x21 y 21 : xþ 2 y2 y1 2ðy 2 y 1 Þ
ð15Þ
The derivative of y for a given x by x1 is: dy x x1 ¼ : dx1 y 2 y 1
ð16Þ
28
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
The change in the area is: Z u2 o dE dy ðu2 u21 Þ=2 x1 ðu2 u1 Þ u2 u1 nu1 þ u2 x1 : ¼ dx ¼ 2 ¼ y2 y1 2 dx1 y2 y1 u1 dx1
ð17Þ
Note that since the points U 1 and U 2 are on line (15), substituting them as ðx; yÞ in (15) leads to the relationship u2 u1 v2 v1 ¼ : ð18Þ y2 y1 x2 x1 Relationship (18) can be used in (17) in the case y 1 ¼ y 2 . Also, simple manipulations using (18) lead to: 2 2 u2 u1 v2 v1 ðu2 u1 Þ2 þ ðv2 v1 Þ2 ¼ ¼ 2 2 y2 y1 x2 x1 ðx2 x1 Þ þ ðy 2 y 1 Þ 1 which means that uy 2 u in (17) is the ratio of the edge’s length and the distance between the two generator points (it may be y 2
1
negative). To apply this relationship one only needs to determine the sign of the ratio. This approach has the advantage of not getting an expression of 0/0 unless the two generator points coincide.
Calculating the derivative of the radius Consider Figs. 2–5. There are three case: 1. The distance is measured to an internal vertex. 2. The distance is measured to a vertex on a side of the square, and 3. The distance is measured to the corner of the square.
Case 1: The distance is measured to an internal vertex Every internal vertex of the Voronoi diagram is at the center of the circle of radius R circumscribing a triangle with three generator points as its vertices. Let the three generator points be ðxi ; y i Þ; i ¼ 1; 2; 3 and we find the derivative of R by x1 . Let the three sides of the triangle be a; b; c, where a is the side opposite ðx1 ; y 1 Þ and a be the angle opposite side a. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ ðx2 x3 Þ2 þ ðy 2 y 3 Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 b ¼ ðx1 x3 Þ þ ðy 1 y 3 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 c ¼ ðx1 x2 Þ þ ðy 1 y 2 Þ : a The circumscribing radius is R ¼ 2 sin . Note that this formula is true also for a > p=2. However, in this case the center of a the circle is outside the triangle but the circle of radius R passes through the three generator points. The circumscribing circle in this case is actually at the center of the longest side of the triangle, and the third generator point is inside that circle. In this context we still call R as the radius of the circumscribing circle because it is equidistant from all three generator points. When the first generator point is moved, a may change, but a (the opposite side) remains unchanged. Therefore,
dR a cos a R ¼ ¼ : da tan a 2 sin2 a We now find
da . dx1
ð19Þ
By direct calculation:
db x1 x3 ¼ b dx1 dc x1 x2 : ¼ c dx1 By the cosine theorem: 2bc cos a ¼ b2 þ c2 a2 ;
ð20Þ
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
29
differentiating (20) by x1 yields: da dc db db dc ¼ 2b þ 2 cos a b þc þ 2c 2bc sin a dx1 dx1 dx1 dx1 dx1
da b c þ cos a ðx1 x2 Þ þ ðx1 x3 Þ ¼ 2ð2x1 x2 x3 Þ bc sin a dx1 c b h i x1 x3 2x1 x2 x3 x1 x2 da cos a c2 þ b2 bc : ¼ ¼ 2x1 x2 x3 sin a dx1
Combining (19) and (21) yields: dR dR da R ¼ ¼ dx1 da dx1 tan a
cos a
h
x1 x2 c2
i 3 þ x1bx 2x1 xbc2 x3 2 sin a
¼
ð21Þ
a cos a 2x1 x2 x3 x1 x2 x1 x3 ; cos a þ bc c2 b2 2 sin3 a
ð22Þ
and similarly
dR a cos a 2y 1 y 2 y 3 y1 y2 y1 y3 ; cos a ¼ þ bc c2 dy 1 2 sin3 a b2
ð23Þ
where cos a ¼
b2 þ c 2 a2 : 2bc
The following Theorem can save 23 of the effort required to calculate all the derivatives. Theorem 3. Let a vertex k be defined by three generator points i1 ; i2 ; i3 . Then, same is true for the derivatives by i2 and i3 .
dRi1 k dxi1
¼
dRi2 k dxi1
¼
dRi3 k dxi1
and
dRi1 k dy i1
¼
dRi2 k dy i1
¼
dRi3 k dy i1 .
The
Proof. Vertex k is the center of the circumscribing circle of the triangle whose vertices are i1 ; i2 , and i3 . Therefore, Ri1 k ¼ Ri2 k ¼ Ri3 k . Following a small change in the location of i1 , vertex k moves to the new center of the triangle and all three distances are still equal to one another. Therefore, the three derivatives are equal to one another. h Case 2: The distance is measured to a vertex on a side of the square Such a vertex is equally distant from two generator points ðx1 ; y 1 Þ and ðx2 ; y 2 Þ. As in Case 1 we find the derivative of the distance by x1 . Consider ðx3 ; y 3 Þ as the mirror image of the second generator point by the side of the square. This mirror image is outside the square. The side of the square is the perpendicular bisector between the second and third points. Therefore, the vertex is equi-distant from all three points. When x1 is moved and the other points remain unchanged, the vertex moves along the side of the square and Eqs. (22) and (23) hold. Case 3: The distance is measured to a corner of the square In this case the derivative is simply the derivative of the distance to a fixed point. Let the corner of the square be at qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dR 1 ðc1 ; c2 Þ. The distance is R ¼ ðx1 c1 Þ þ ðy 1 c2 Þ and the derivative is dx ¼ x1 c . R 1 References Baron, O., Berman, O., Krass, D., Wang, Q., 2007. The equitable location problem on the plane. European Journal of Operational Research 183, 578–590. Berman O., Drezner, Z., Tamir, A., Wesolowsky, G.O., in press. Optimal Location with Equitable Loads, Annals of Operations Research. Beyer, W.H., 1981. CRC Standard Mathematical Tables, 26th ed. CRC Press, Boca Raton, FL. Current, J., Daskin, M., Schilling, D., 2002. Discrete network location models, In: Drezner, Z., Hamacher, H.W. (Eds.), Facility Location: Applications and Theory, pp. 81–118. Daskin, M.S., 1995. Network and Discrete Location: Models, Algorithms, and Applications. John Wiley & Sons, New York. Drezner, Z., 1984. The p-center problem – heuristic and optimal algorithms. Journal of the Operational Research Society 35, 741–748. Drezner, Z., 1987. On the rectangular p-center problem. Naval Research Logistics Quarterly 34, 229–234. Drezner, T., Drezner, Z., 2006. Multiple facilities location in the plane using the gravity model. Geographical Analysis 38, 391–406. Handler, G.Y., 1990. p-Center problems. In: Mirchandani, P.B., Francis, R.L. (Eds.), Discrete Location Theory. Wiley Inter-Science, New York, pp. 305– 347. Ohya, T., Iri, M., Murota, K., 1984. Improvements of the incremental method of the Voronoi diagram with computational comparison of various algorithms. Journal of the Operations Research Society of Japan 27, 306–337.
30
A. Suzuki, Z. Drezner / European Journal of Operational Research 195 (2009) 17–30
Okabe, A., Boots, B., Sugihara, K., Chiu, S.N., 2000. Spatial tessellations: Concepts and applications of Voronoi diagrams. Wiley Series in Probability and Statistics. Okabe, A., Suzuki, A., 1987. Stability of spatial competition for a large number of firms on a bounded two-dimensional space. Environment and Planning A 16, 107–114. Sugihara, K., Iri, M., 1994. A robust topology-oriented incremental algorithm for Voronoi diagram. International Journal of Computational Geometry and Applications 4, 179–228. Suzuki, A., Drezner, Z., 1996. The p-center location problem in an area. Location Science 4, 69–82. Suzuki, A., Okabe, A., 1995. Using Voronoi diagrams. In: Drezner, Z. (Ed.), Facility Location: A Survey of Applications and Methods. Springer, New York. Vijay, J., 1985. An algorithm for the p-center problem in the plane. Transportation Science 19, 235–245.