- Email: [email protected]
Siting a facility in continuous space to maximize coverage of a region
Socio-Economic Planning Sciences 43 (2009) 131–139
Contents lists available at ScienceDirect
Socio-Economic Planning Sciences journal homepage: www.elsevier.com/locate/seps
Siting a facility in continuous space to maximize coverage of a region Timothy C. Matisziw a, b, *, Alan T. Murray c a b c
Department of Geography, University of Missouri-Columbia, Columbia, MO 65211-6170, United States Department of Civil and Environmental Engineering, University of Missouri-Columbia, Columbia, MO 65211-6170, United States School of Geographical Sciences, Arizona State University, Tempe, AZ 85287-0104, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Available online 12 June 2008
Siting facilities in continuous space such that continuously distributed demand within a region is optimally served is a challenging location problem. This problem is further complicated by the non-convexity of regions typically encountered in practice. In this paper a model for maximizing the service coverage of continuously distributed demand through the location of a single service facility in continuous space is proposed. To address this problem, theoretical conditions are established and associated methods are proposed for optimally siting a service facility in a region (convex or non-convex) with uniformly distributed demand. Through the use of geographic information systems (GIS), the developed approach is applied to identify facility sites that maximize regional coverage provided limitations on facility service ability. 2008 Elsevier Ltd. All rights reserved.
Keywords: Medial axis Maximally enclosed disks Non-convex region Spatial optimization GIS
1. Introduction Many planning problems involve siting facilities in order to serve (or cover) regional demand. Often complete coverage of all demand in a region is not possible due to budgetary limitations on the number of facilities that can be sited. Thus, limited resources must be efficiently managed and regional demand should be served/covered to the greatest extent possible (see Ref. [7]). Siting studies where this planning dilemma occurs include placement of emergency warning sirens [9], bus stops [14], cellular towers [17], emergency response stations [1,7,36], automatic meter reading stations [13], radio receiving stations and other types of signal transmitting equipment [3,23]. Of interest in this paper, is locating a facility in continuous space to serve continuously distributed demand. In other words, demand is considered present everywhere within a region. As an example, when locating emergency
* Corresponding author. Department of Geography, University of Missouri-Columbia, Columbia, MO 65211-6170, United States Tel.: þ1 573 882 8370; fax: þ1 573 884 4239. E-mail addresses: [email protected] (T.C. Matisziw), atmurray@ asu.edu (A.T. Murray). 0038-0121/$ – see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.seps.2008.02.009
warning sirens, it is important for sirens to be audible in residential and commercial areas, but schools, places for outdoor recreation, transportation corridors, and anticipated development areas may be equally important to serve [9,26]. Further, it is assumed that the service facility can be sited anywhere within a region as opposed to traditional approaches where candidate sites are identified a priori. This is a realistic assumption given the relatively small geographic footprint of many types of service facilities. For instance, siting emergency warning sirens and cellular equipment is very flexible since these types of facilities can be mounted on posts or existing structures. Provided these assumptions, the planning problem of interest here is to locate a single facility in continuous space to maximize coverage of continuously distributed demand. While there is considerable literature focusing on special cases of this problem, as will be discussed next, it has not been explicitly approached in prior research. In this paper, a model is developed for addressing the location problem of siting a facility in a region such that coverage of continuous demand is optimized. First, the ways in which this problem has been approached in the past and the limitations that have been encountered are discussed. In the next section, the location problem is
132
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
formalized. Following this, several geometric properties of regions are explored, contributing to a solution approach for the problem. To illustrate the solution approach, the developed methodology is applied and the results are discussed for locating a single facility in a non-convex region. Finally, conclusions are given and future research directions are discussed. 2. Background There has been considerable work oriented toward siting facilities that have specific geographic ranges within which their service is most effective. Additionally, prior research has recognized that complete coverage of regional demand is not always possible given budgetary limitations. To address these concerns when facility and demand locations are discrete, Church and ReVelle [7] proposed the maximal covering location problem (MCLP), a linearinteger model that seeks to maximize demand coverage subject to limited resources. The objective of the MCLP is to maximize suitable coverage of demand. Constraints are employed to limit the number of facilities sited and to ensure that a demand location cannot be considered served unless a facility capable of serving it is sited. Suitable service is typically established using a time or distance standard S from candidate facility sites, within which effective service can be provided to demand locations (i.e., dij S, where dij is the shortest distance or travel time between demand i and candidate site j). Application of the MCLP typically assumes that candidate facility sites and demand locations are represented as sets of discrete points. An extension for siting facilities in continuous space in order to maximally cover discrete demand locations has been approached by Mehrez and Stulman [22,23], Mehrez [21], and Church [8]. In particular, Church [8] refers to this problem as the planar maximal covering problem (PMCE) when the Euclidean metric is used to demarcate service provision. In order to solve the PMCE, properties of coverage are exploited to discretize this continuous space problem. To identify candidate facility locations in the continuous domain, the geometry of the location problem is used to identify a set of potential locations that will contain optimal facility sites. This set of candidate locations is formed by finding the intersection of all circles of radius S centered at each discrete demand location. These points are referred to as the circle intersection point set (CIPS) [8]. Other problem variants, such as minimal covering objectives and the use of varying service radii, are discussed in Plastria [30]. Though the CIPS approach addresses continuous facility location, demand is still assumed to be discrete. Representing demand in a discrete fashion is known to introduce a variety of inaccuracies into the location modeling process [12,20]. Models, such as the continuous p-center problem, have been introduced for locating facilities in the plane (where facilities may be sited anywhere) to address continuously distributed demand [25,34,35]. The objective of the p-center problem is to minimize the maximum distance from any demand location i to its nearest sited facility j. However, when sited facilities have specific service
standards, complete coverage cannot always be ensured for a given value of p. Complete coverage can only be guaranteed if min fdij g is less than or equal to S for each demand j area i. In the case where demand is not completely covered for a specified value of p, no attempt is made to maximize the demand that is covered. Though the PMCE and p-center problems do incorporate aspects of continuous facility location, they do not simultaneously address service facility siting anywhere in a region and coverage maximization of continuously distributed demand.
3. Maximizing coverage of continuous demand Representing demand as continuously distributed, in contrast to discretely defined in space, is a major complicating feature in modeling this particular location problem. The goal here is to identify the location(s) where a facility with a specific service area would provide the most coverage of a region (G). Fig. 1 depicts a non-convex region G, not unlike that encountered in planning applications, where maximal coverage of regional demand is desired. It is assumed that demand is distributed within region G according to some continuous function d (G). All locations within G are considered potential facility sites. Thus, continuously distributed demand is in need of coverage by a facility located in the plane. The problem may be stated as follows: Continuous maximal covering problem: 1-facility (CMCP-1)
Max ðx;yÞ
XZZ k
dðGÞ dG;
(1)
Uk
where (x,y) ¼ decision variables for geographic coordinates of facility, V ¼ service area of facility sited at (x,y), U ¼ VXG, k ¼ index for continuous components of U (entire set K), Uk ¼ kth continuous component of UðW Uk ¼ UÞ. k The objective of the CMCP-1, Eq. (1), is to maximize the amount of continuous demand d(G) in G that is suitably served (or covered). The decision to be made in the CMCP-1 is where to site a single facility (determine the facility’s x, y coordinates). Demand is considered covered if it is within the service area (V) of the sited facility. U represents the portion(s) of the region that are provided coverage by the facility. This non-linear, non-convex model involves integrating continuous demand over the coverage area U. Since U need not be a single, continuous area, Uk references individual components of U should discontinuities exist. The primary difference between the CMCP-1 and its discrete counterparts is that it allows demand and candidate facility locations to exist everywhere within a region. Thus, in the CMCP-1, the cardinality of the sets of demand (i ˛ I) and facility locations (j ˛ J) is infinite (jIj / N and jJj / N), in contrast to the discrete and finite number of demand and facility locations in the MCLP and discrete demand in the PMCE. In order to solve the maximization problem in Eq. (1), one faces the impossible task of needing to evaluate every candidate facility location (x,y) within the plane for service
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
133
Fig. 1. A non-convex region.
areas V and determining the demand covered. After all locations have been evaluated, the location(s) providing coverage to the greatest amount of demand would represent the optimal site(s). Clearly, the CMCP-1 cannot be approached using standard optimization techniques due to the infinite number of facility and demand locations involved, and the fact that the model can be non-convex, non-linear, and discontinuous! Drezner and Drezner [12] do show how some location problems considering continuous demand can be solved using double integration when the region is convex (square in their case). However, the problem of interest in this paper deals with regions (either convex or non-convex) and their intersection with a facility’s service area, whose integral calculations would be unwieldy, even if approached via quadrature. As a result, it is evident why researchers have focused on special cases of the CMCP-1, like the MCLP and PMCE, that can be discretized, enabling a linear-integer solution approach. Fortunately, it may be possible to exploit geometric properties of a region to be covered in order to establish optimality properties/conditions for the CMCP-1. It is assumed for the remainder of this paper that demand is uniformly distributed in region G (i.e., s ¼ d(G)). As discussed previously, uniformly distributed demand often reflects the realities of many planning problems, like siting warning sirens or radar stations in an urban area, as all areas are equally important to serve. Another assumption is that the service area of the facility (V) is defined as the area within S (as measured in Euclidean distance) of a point centered at (x,y). There are two distinct cases of the CMCP-1. First, the facility’s service area can be completely contained within the region G. Second, the facility’s service area cannot be completely contained within the region. Solution approaches based on these two cases are detailed in Section 5. 4. Geometric representation of regions: the medial axis Before getting into the specifics of the two cases of the CMCP-1, it is necessary to review a geometrical
representation of a region known as the medial axis. Brandt [4] details some important properties of polygons in the plane, and they will be used to define the medial axis. First, let L be defined as the set of all pairs (q,r), where q is a point and r is the maximum radius of a disk centered on the point that is completely contained in region G, or equivalently, the minimum distance from the point to the complement of the region G c (d(q,G c)). In other words, L ¼ fðq; rÞj Cðq; rÞ4Gg, meaning that every location in G is considered to be a center of the largest disk possible that is still completely contained within G. For example, Fig. 2b shows a subset of those points defined in L whose disks are fully contained within the regional boundary given in Fig. 2a. From this definition, it can be shown that the union of the disks centered on these points (and completely contained in G) will exactly regenerate G [4]. In other words,
G [ W Cðq; rÞ: ðq;rÞ˛L
(2)
As can be seen in Fig. 2b, every point in the region will be the center of a disk having at least one point in common with the region’s boundary. Although only a sampling of points are depicted in this figure, it is evident that the edges of their associated disks are already beginning to approximate the general shape of G. As defined above, every point q˛L has an associated radius and it follows that every resulting disk will also be contained within or dominated by a disk of greater or equivalent radius, ðq; rÞ ðq 0 ; r 0 Þ5Cðq; rÞ4Cðq 0 ; r 0 Þ. Given this relationship among points centered in G, the set of points having disks that are maximal is identified as [4]
L [ fðq; rÞ˛L : ðq; rÞ £ ðq 0 ; r 0 Þ˛L0ðq; rÞ [ ðq 0 ; r 0 Þg: Fig. 2c shows a subset of maximally enclosed disks within this region. L* maintains the additional property that each of the associated disks have at least two points in common with the region’s boundary, meaning that their center is equidistant to their two nearest boundary locations. Since all of the disks in L are inferior or equivalent
134
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
and cartographic label placement in map polygons [11]. In the following sections the properties of the medial axis are relied upon in the solution of the CMCP-1.
a
5. Coverage properties for single facility siting The definition of a region’s medial axis is important for the CMCP-1 because every point on the medial axis is the center of a disk of maximal size completely contained within the region. This definition ensures that, under certain conditions, a disk centered on the medial axis will definitely cover as much or more area (uniformly distributed continuous demand in this case) than any other disk within the region. This property of the medial axis is consistent with the objective of the CMCP-1 and can be exploited to solve the CMCP-1. At this point several important observations are necessary relative to the special cases of the CMCP-1.
b
c
d
Fig. 2. Development of the medial axis. (a) A simple, non-convex region, (b) enclosed disks centered at points within the region, (c) maximally enclosed disks centered at points within the region and (d) the union of the centers of all maximally enclosed disks (the medial axis).
to those in L*, the reduced set L* can also be used to exactly regenerate region G as follows:
G[
W Cðq; rÞ: ðq;rÞ˛L
(3)
This relationship establishes that a region can be represented as a locus of points, while still retaining all of the information necessary to exactly regenerate the region. Fig. 2d shows the centers of the maximally enclosed disks for the region in Fig. 2a. This unique set of disk centers, L*, is commonly referred to as the skeleton or medial axis in computational geometry. Brandt [4] formally defines the medial axis as MðGÞ ¼ fq : dr; ðq; rÞ˛L)g. Various techniques exist to generate the medial axis for a region, such as wavefront propagation [2,24] and Voronoi diagram based approaches [6,18,19,28,31]. Given that the medial axis of a polygon is an alternative representation of space, it has found many useful applications in geographical analysis. Some of these uses include point-in-polygon tests [29], statistical analysis of points within a region [33], routing objects through an enclosed space [37], topology generation for scanned maps [15,16],
Observation 1: The smallest service standard needed to ensure complete coverage of a region can be obtained by solving a 1-center problem to identify the center and radius of the smallest disk that can completely enclose the region (see Refs. [10,35]). That is, for service standards greater than or equal to the radius of the smallest enclosing disk, complete regional coverage is always ensured when siting a facility at the 1-center. Observation 2: If S is less than or equal to the maximum radius of the maximally enclosed disks centered on the medial axis (rmax), at least one optimal facility site exists in G. To see this, consider any facility j sited in G. If the service standard, S, for the facility is less than the minimum distance between j and G c, the service area of the facility will be completely contained in G and will provide optimal coverage. Move the facility anywhere in G such that S dðj; G c Þ rmax and the solution is equivalent. Since a disk along the medial axis is maximal with respect to all other locally enclosed disks in G, the site(s) j along the medial axis with the maximum d(j,G c) will correspond to the largest service area that can be completely contained within the region. Observation 3: If the service standard is greater than the radius of the maximum maximally enclosed disk(s) and less than the radius of the smallest enclosing disk associated with the 1-center solution for the region, sites providing optimal coverage will be influenced by Gc. Consider any site j on the medial axis. If the service standard associated with this point is d(j,Gc), any increase to its radius will result in the service area no longer being enclosed in G. A service area, V, not enclosed in G may not necessarily be convex or continuous with respect to VXG. Given this, assessment of whether j maximizes coverage of G depends not only on its ability to provide maximal local coverage, but also on its ability to cover Gc. In this regard, the local properties of the medial axis do not necessarily coincide with optimal siting locations. 6. Solution procedure This section proposes techniques to address the two cases outlined in Observations 2 and 3 for the CMCP-1. An
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
optimal solution to the CMCP-1 can be found for cases where S rmax (Observation 2). This task involves finding the enclosed disk(s) and associated facility location(s) that provides maximal demand coverage of the region. If a facility’s service standard is less than or equal to rmax, the medial axis will always contain an optimal solution. Additionally, when S ¼ rmax where d(j,G c) is locally maximal, the optimal solution will exist only on the medial axis. A solution approach for identifying an optimal location along the medial axis is as follows: CMCP-1 Algorithm I. (Case 1: S rmax) 1 (i) Generate a line-Voronoi diagram for G to obtain medial axis J. (ii) Search J assessing d(j,G c) to find optimal point(s) or line segment(s) where d(j,G c) S. (iii) Apply line-search (or golden section search) technique to identify rmax. In order to find solutions to the CMCP-1 when S > rmax, but less than that the radius of the smallest enclosing disk (Observation 3), a different technique is needed. In this case, any optimal service area having a radius smaller than that of the smallest enclosing disk but greater than that of the maximum maximally enclosed disk must intersect the region’s boundary at multiple points. A reasonable approach for identifying sites that maximize coverage is to search the boundary for these points since an intersection of circles of radius S centered on them will return the center of the optimal service area. The following procedure is one way for addressing this case: CMCP-1 Algorithm II. (Case 2: S > rmax) 2 (i) Given the region’s boundary H, approximate H by the set of points Ha, where a ¼ sampling density. (ii) For each h ˛ Ha, generate a disk of radius S (C(h, S)). (iii) Identify the intersection points of all resulting disk boundaries (denoted vCðh; SÞ) IðHa Þ ¼ X vCðh; SÞ. h˛Ha
(iv) For each intersection point c˛IðHa Þ, generate a disk of radius S ðCðc; SÞÞ. (v) For each resulting disk, compute the area of intersection between the disk and the region G ðCðc; SÞXGÞ, tracking the c generating the disk that covers the most regional demand (cmax). (vi) As a/M (where M is a very large number) the center of the disk Cðcmax ; SÞ is a location maximizing demand coverage.
135
demand (boundary shown in Fig. 1). The objective, Eq. (1), is to cover as much of this region’s demand in locating a single service facility. Implementing the procedures discussed in the previous section, optimal sites for cases where S rmax were computed using a commercial GIS.3 Fig. 3a–d depict sites in G that are optimal for S ¼ 25, 75, 100, and 150 m. For instance, optimal locations for a service standard of 25 m are those sites at least 25 m from any location in the region’s complement (Gc), as established in Observation 2 (Fig. 3a). Given that multiple optima are often present, a GIS is ideal for visualizing such solutions. It is apparent that as S increases, the number of optimal facility sites decreases until only locations along some portions of the medial axis are optimal. In other words, the medial axis will always contain optimal sites for service standards that are less than or equal to rmax. The panels in Fig. 3 also depict optimal portions of the region’s medial axis. Since the region is non-convex, the medial axis comprised line segments and parabolas (caused by concave angles on the region’s boundary). For convex regions, the medial axis consists of only line segments since all angles on the boundary are convex [5]. Searching the medial axis for the sites with the maximum d(j,Gc) results in the identification of a line segment comprising the central trunk of the medial axis (Fig. 4). These points represent the sites at which the largest, fully enclosed service areas may be centered in G. That is, no other sites in G have a larger service area contained in the region than those along the bold line segment shown in Fig. 4. The radii of the maximally enclosed disks, rmax, centered on these points are 163 m. Accordingly, if the service standard for a facility to be sited was 163 m, these locations would represent the set of optimal candidate facility sites, anyone of which optimizes Eq. (1). The union of all of the maximal service areas is also shown in Fig. 4. To establish an upper bound on the largest service standard for which complete regional coverage cannot be achieved, a 1-center problem was solved for region G using the MINIDISC algorithm described in de Berg et al. [10].4 In this case, the maximum distance between the 1-center and the region’s boundary is 306.56 m, representing the radius of the smallest disk capable of containing every point in G (Fig. 5). For service standards less than the upper bound of 306.56 but greater than 163 m, the CMCP-1 Algorithm II (Case 2) outlined in the previous section is necessary to search for a facility site maximizing regional coverage. Fig. 6 illustrates the best sites found for service standards in the range
7. Assessing regional coverage The developed approaches for solving the CMCP-1 are applied to a single non-convex region with uniform
1 Given an n-edge polygon, computational effort associated with obtaining the medial axis is on the order of O(n log n) [19]. Line searches converge linearly with respect to computational effort [32]. 2 Given m ¼ jHaj, then disk intersections may be constructed in O (m2) time. For two polygons involving n vertices and k edge intersections, overlay can be computed in O(n log n þ k) time [27].
3 A true line-Voronoi diagram generally cannot be computed in a GIS (since arcs are generalized as line segments), so an approximation must be used. Here, a suitably dense point-based approximation of the region’s boundary was used and a Delaunay triangulation was computed from these points. The connected centers of the triangles falling completely within the region’s boundary constitute the elements of the medial axis (see Ref. [28]). Approximations of the medial axis, such as that employed here, are shown to converge to the exact medial axis as the number of sampling points along the boundary approaches infinity [5]. 4 1-Center problem programmed in Avenue ArcView GIS version 3.2.
136
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
a
b
c
d
Fig. 3. Optimal candidate facility sites for selected values of S. (a) Optimal sites for S ¼ 25 m, (b) optimal sites for S ¼ 75 m, (c) optimal sites for S ¼ 100 m and (d) optimal sites for S ¼ 150 m.
163 < S < 306.56. In this case, the results show that as the service standard increases above the maximum radius of the maximally enclosed disks (163 m), the center of the selected service area moves in a non-linear trajectory toward the 1-center site. This type of non-linear movement clearly illustrates the influence of the region’s form on the location of the best site.
8. Conclusion This paper developed a methodology for modeling continuous facility location and continuously distributed demand to maximize service coverage in a region when locating a single service facility. The 1-facility continuous maximal covering problem was proposed as
Legend Region Boundary Medial Axis Centers of Maximal Disks Maximal Disks (r = 163 meters) 0
70
140
210
Meters
Fig. 4. The union of the maximum maximally enclosed disks identified on the medial axis.
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
137
Legend Region Boundary 1-Center Location 1-Center Service Area 0
70
140
210
Meters
Fig. 5. 1-center solution and associated service area.
a generalization of its discrete counterpart, the maximal covering location problem. Unlike the MCLP, the CMCP-1 is a non-linear and non-convex model. Given that complete enumeration, linear, or non-linear programming are not feasible options for identifying exact solutions to this type of problem, a means for deriving the optimal solution to the CMCP-1 was proposed. The solution approach is based on the idea that it is not feasible to search the entire region for the facility sites that provide maximal coverage. Instead, an equivalent line representation of a region, known as the medial axis or skeleton, is exploited. The medial axis is an exact representation of a polygon and it is shown that under certain conditions optimal solutions to the CMCP-1 exist on the medial axis. Furthermore, this technique allows for the identification of multiple optima (if present), which often are
desirable for planning purposes. However, in the presence of multiple optima, siting on the medial axis could present additional coverage opportunities given future upgrading of a service facility’s service standard is possible. Consideration of a region’s medial axis definitely provides needed spatial structure to the search for optimal solutions. Dealing with the case where S is greater than the maximum radius of the maximally enclosed disks presents yet another challenge for assessing regional coverage of sited facilities. Here, this issue was addressed through the generation of intersection points derived from circles centered on the region’s boundary and enumerating the regional coverage of each resulting service area. As was described, this procedure need only be performed for service standards greater than the radius of the largest
Fig. 6. Identified facility sites for service standards greater than the maximum radius of the maximally enclosed disks.
138
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
maximally enclosed disk and less than the radius of the smallest enclosing disk. All of the procedures discussed in this paper are based upon the geometrical properties of a region and are easily assessed using GIS as implemented in the example provided. Hence, an additional benefit of this spatial optimization approach is that it is readily accessible to anyone with a GIS, further highlighting important contributions from GIScience for addressing spatial planning problems. This work represents an initial step in addressing the solution properties/techniques for continuous maximal covering problems, and further work is possible. Future work could involve extending the developed methodology for siting multiple facilities in continuous space to serve continuous demand. Furthermore, this paper only considers facility service areas that are circular in shape. Therefore, exploring facility siting given a broader interpretation of service area potential is of interest. Finally, relaxing the assumption of uniform demand would be a major contribution to this area of research.
Acknowledgements This material is based upon work supported by the National Science Foundation under Grant No. 0518967. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. References [1] Adenso-Diaz B, Rodriguez F. A simple search heuristic for the MCLP: application to the location of ambulance bases in a rural region. Omega 1997;25(2):181–7. [2] Blum H. A transformation for extracting new descriptors of shape. In: Walthen-Dunn W, editor. Models for the perception of speech and visual form. Cambridge: MIT Press; 1967. p. 362– 80. [3] Brady SD, Rosenthal RE. Interactive computer graphical solutions of constrained minimax location problems. AIIE Transactions 1980 September:241–8. [4] Brandt JW. Convergence and continuity criteria for discrete approximations of the continuous planar skeleton. CVGIP: Image Understanding 1994;59(1):116–24. [5] Brandt JW, Algazi VR. Continuous skeleton computation by Voronoi diagram. CVGIP: Image Understanding 1992;55(3):329–38. [6] Chou JJ. Voronoi diagrams for planar shapes. IEEE Computer Graphics and Applications 1995 March:52–9. [7] Church RL, ReVelle CS. The maximal covering location problem. Papers – Regional Science Association 1974;32:101–18. [8] Church RL. The planar maximal covering location problem. Journal of Regional Science 1984;24(2):185–201. [9] Current J, O’Kelly M. Locating emergency warning sirens. Decision Sciences 1992;23:221–34. [10] de Berg M, van Kreveld M, Overmars M, Schwarzkopf O. Computational geometry: algorithms and applications. 2nd ed. New York: Springer-Verlag; 2000. [11] Do¨rschlag D, Petzold I, Plu¨mer L. Placing objects automatically in areas of maps. In: Proceedings of the Durban ICC 2003, South Africa. [12] Drezner T, Drezner Z. Replacing continuous demand with discrete demand in a competitive location model. Naval Research Logistics 1997;44:81–95. [13] Gavirneni S, Clark L, Pataki G. Schlumberger optimizes receiver location for automated meter reading. Interfaces 2004;34(3): 208–14. [14] Gleason JM. A set covering approach to bus stop location. Omega 1975;3:605–8.
[15] Gold C, Snoeyink J. A one-step crust and skeleton extraction algorithm. Algorithimica 2001;30:144–63. [16] Gold CM, Nantel J, Yang W. Outside-in: an alternative approach to forest map digitizing. International Journal of Geographical Information Systems 1996;10(3):291–310. [17] Kalvenes J, Kennington J, Olinick E. Hierarchical cellular network design with channel allocation. European Journal of Operational Research 2004;160(1):3–18. [18] Kirkpatrick DG. Efficient computation of continuous skeletons. In: Proceedings of the 20th annual symposium on foundations of computer science; 1979. p. 18–27. [19] Lee DT. Medial axis transformation of a planar shape. IEEE Transactions on Pattern Analysis and Machine Intelligence 1982;PAMI-4(4): 363–9. [20] Love RF. A computational procedure for optimally locating a facility with respect to several rectangular regions. Journal of Regional Science 1972;12(2):233–42. [21] Mehrez A. A note on the linear integer formulation of the maximal covering location problem with facility placement on the entire plane. Journal of Regional Science 1983;23(4):553–5. [22] Mehrez A, Stulman A. The maximal covering location problem with facility placement on the entire plane. Journal of Regional Science 1982;22(3):361–5. [23] Mehrez A, Stulman A. An extended continuous maximal covering location problem with facility placement. Computers & Operations Research 1984;11(1):19–23. [24] Montanari U. Continuous skeletons from digitized images. Journal of the Association for Computing Machinery 1969;16(4): 534–49. [25] Murray AT, O’Kelly ME, Church RL. Regional service coverage modeling. Computers & Operations Research 2008;35(2):339– 55. [26] O’Kelly ME, Murray AT. A lattice covering model for evaluating existing service facilities. Papers in Regional Science 2004;83: 565–80. [27] O’Rourke J. Computational geometry in C. 2nd ed. Cambridge: Cambridge University Press; 1998. [28] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts & applications of Voronoi diagrams. 2nd ed. New York: John Wiley & Sons, Ltd.; 2000. [29] Pfaltz JL, Rosenfeld A. Computer representation of planar regions by their skeletons. Communications of the ACM 1967;10(2): 119–25. [30] Plastria F. Continuous covering location problems. In: Drezner Z, Hamacher H, editors. Facility location: applications and theory. New York: Springer; 2002. [31] Preparata FP. The medial axis of a simple polygon. In: Goos G, Hartmanis J, editors. Mathematical foundations of computer science, lecture notes in computer science, vol. 53. New York: SpringerVerlag; 1977. p. 443–50. [32] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C: the art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [33] Sadahiro Y, Takami K. A statistical method for analyzing the relative location of points in a bounded region. Geographical Analysis 2001; 33(1):40–57. [34] Suzuki A, Okabe A. Using Voronoi diagrams. In: Drezner Z, editor. Facility location: a survey of applications and methods. New York: Springer-Verlag; 1995. p. 103–18. [35] Suzuki A, Drezner Z. The p-center location problem in an area. Location Science 1996;4(1/2):69–82. [36] Toregas C, Swain R, ReVelle C, Bergman L. The location of emergency service facilities. Operations Research 1971;19: 1363–73. [37] Xu Y, Mattikalli R, Khosla P. Determining two circumscribing discs for a polygon. Computer-Aided Design 1995;27(3):223–34.
Tim Matisziw received his B.A. and M.A. in Geography from the University of Missouri – Columbia and a Ph.D. in Geography from The Ohio State University. Dr. Matisziw is currently Assistant Professor, Department of Geography and Department of Civil & Environmental Engineering at the University of Missouri – Columbia. His current research addresses modeling issues in transportation, location science, and spatial optimization and has appeared in journals such as the European Journal of Operational Research, Computers & Operations Research, Decision Support Systems, International Regional Science Review, Computers, Environment and Urban Systems, Journal of Geographical Systems, Transportation Research C, and Regional Studies.
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139 Alan Murray obtained a B.S. in mathematical sciences, an M.A. in statistics and applied probability and a Ph.D. in geography, all from the University of California at Santa Barbara. Dr. Murray is Professor, School of Geographical Sciences, at Arizona State University. His research and teaching interests are in spatial optimization modeling, geographic information science, urban/regional planning and development, and transportation. He has
139
published on a range of technical and application oriented topics in journals such as the Operations Research, Discrete Applied Mathematics, International Journal of Geographic Information Sciences, Geographical Analysis, Journal of Geographical Systems, Journal of Urban Planning and Development, European Journal of Operational Research, Computers and Operations Research, Transportation Research, and Urban Studies, among others.
Contents lists available at ScienceDirect
Socio-Economic Planning Sciences journal homepage: www.elsevier.com/locate/seps
Siting a facility in continuous space to maximize coverage of a region Timothy C. Matisziw a, b, *, Alan T. Murray c a b c
Department of Geography, University of Missouri-Columbia, Columbia, MO 65211-6170, United States Department of Civil and Environmental Engineering, University of Missouri-Columbia, Columbia, MO 65211-6170, United States School of Geographical Sciences, Arizona State University, Tempe, AZ 85287-0104, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Available online 12 June 2008
Siting facilities in continuous space such that continuously distributed demand within a region is optimally served is a challenging location problem. This problem is further complicated by the non-convexity of regions typically encountered in practice. In this paper a model for maximizing the service coverage of continuously distributed demand through the location of a single service facility in continuous space is proposed. To address this problem, theoretical conditions are established and associated methods are proposed for optimally siting a service facility in a region (convex or non-convex) with uniformly distributed demand. Through the use of geographic information systems (GIS), the developed approach is applied to identify facility sites that maximize regional coverage provided limitations on facility service ability. 2008 Elsevier Ltd. All rights reserved.
Keywords: Medial axis Maximally enclosed disks Non-convex region Spatial optimization GIS
1. Introduction Many planning problems involve siting facilities in order to serve (or cover) regional demand. Often complete coverage of all demand in a region is not possible due to budgetary limitations on the number of facilities that can be sited. Thus, limited resources must be efficiently managed and regional demand should be served/covered to the greatest extent possible (see Ref. [7]). Siting studies where this planning dilemma occurs include placement of emergency warning sirens [9], bus stops [14], cellular towers [17], emergency response stations [1,7,36], automatic meter reading stations [13], radio receiving stations and other types of signal transmitting equipment [3,23]. Of interest in this paper, is locating a facility in continuous space to serve continuously distributed demand. In other words, demand is considered present everywhere within a region. As an example, when locating emergency
* Corresponding author. Department of Geography, University of Missouri-Columbia, Columbia, MO 65211-6170, United States Tel.: þ1 573 882 8370; fax: þ1 573 884 4239. E-mail addresses: [email protected] (T.C. Matisziw), atmurray@ asu.edu (A.T. Murray). 0038-0121/$ – see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.seps.2008.02.009
warning sirens, it is important for sirens to be audible in residential and commercial areas, but schools, places for outdoor recreation, transportation corridors, and anticipated development areas may be equally important to serve [9,26]. Further, it is assumed that the service facility can be sited anywhere within a region as opposed to traditional approaches where candidate sites are identified a priori. This is a realistic assumption given the relatively small geographic footprint of many types of service facilities. For instance, siting emergency warning sirens and cellular equipment is very flexible since these types of facilities can be mounted on posts or existing structures. Provided these assumptions, the planning problem of interest here is to locate a single facility in continuous space to maximize coverage of continuously distributed demand. While there is considerable literature focusing on special cases of this problem, as will be discussed next, it has not been explicitly approached in prior research. In this paper, a model is developed for addressing the location problem of siting a facility in a region such that coverage of continuous demand is optimized. First, the ways in which this problem has been approached in the past and the limitations that have been encountered are discussed. In the next section, the location problem is
132
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
formalized. Following this, several geometric properties of regions are explored, contributing to a solution approach for the problem. To illustrate the solution approach, the developed methodology is applied and the results are discussed for locating a single facility in a non-convex region. Finally, conclusions are given and future research directions are discussed. 2. Background There has been considerable work oriented toward siting facilities that have specific geographic ranges within which their service is most effective. Additionally, prior research has recognized that complete coverage of regional demand is not always possible given budgetary limitations. To address these concerns when facility and demand locations are discrete, Church and ReVelle [7] proposed the maximal covering location problem (MCLP), a linearinteger model that seeks to maximize demand coverage subject to limited resources. The objective of the MCLP is to maximize suitable coverage of demand. Constraints are employed to limit the number of facilities sited and to ensure that a demand location cannot be considered served unless a facility capable of serving it is sited. Suitable service is typically established using a time or distance standard S from candidate facility sites, within which effective service can be provided to demand locations (i.e., dij S, where dij is the shortest distance or travel time between demand i and candidate site j). Application of the MCLP typically assumes that candidate facility sites and demand locations are represented as sets of discrete points. An extension for siting facilities in continuous space in order to maximally cover discrete demand locations has been approached by Mehrez and Stulman [22,23], Mehrez [21], and Church [8]. In particular, Church [8] refers to this problem as the planar maximal covering problem (PMCE) when the Euclidean metric is used to demarcate service provision. In order to solve the PMCE, properties of coverage are exploited to discretize this continuous space problem. To identify candidate facility locations in the continuous domain, the geometry of the location problem is used to identify a set of potential locations that will contain optimal facility sites. This set of candidate locations is formed by finding the intersection of all circles of radius S centered at each discrete demand location. These points are referred to as the circle intersection point set (CIPS) [8]. Other problem variants, such as minimal covering objectives and the use of varying service radii, are discussed in Plastria [30]. Though the CIPS approach addresses continuous facility location, demand is still assumed to be discrete. Representing demand in a discrete fashion is known to introduce a variety of inaccuracies into the location modeling process [12,20]. Models, such as the continuous p-center problem, have been introduced for locating facilities in the plane (where facilities may be sited anywhere) to address continuously distributed demand [25,34,35]. The objective of the p-center problem is to minimize the maximum distance from any demand location i to its nearest sited facility j. However, when sited facilities have specific service
standards, complete coverage cannot always be ensured for a given value of p. Complete coverage can only be guaranteed if min fdij g is less than or equal to S for each demand j area i. In the case where demand is not completely covered for a specified value of p, no attempt is made to maximize the demand that is covered. Though the PMCE and p-center problems do incorporate aspects of continuous facility location, they do not simultaneously address service facility siting anywhere in a region and coverage maximization of continuously distributed demand.
3. Maximizing coverage of continuous demand Representing demand as continuously distributed, in contrast to discretely defined in space, is a major complicating feature in modeling this particular location problem. The goal here is to identify the location(s) where a facility with a specific service area would provide the most coverage of a region (G). Fig. 1 depicts a non-convex region G, not unlike that encountered in planning applications, where maximal coverage of regional demand is desired. It is assumed that demand is distributed within region G according to some continuous function d (G). All locations within G are considered potential facility sites. Thus, continuously distributed demand is in need of coverage by a facility located in the plane. The problem may be stated as follows: Continuous maximal covering problem: 1-facility (CMCP-1)
Max ðx;yÞ
XZZ k
dðGÞ dG;
(1)
Uk
where (x,y) ¼ decision variables for geographic coordinates of facility, V ¼ service area of facility sited at (x,y), U ¼ VXG, k ¼ index for continuous components of U (entire set K), Uk ¼ kth continuous component of UðW Uk ¼ UÞ. k The objective of the CMCP-1, Eq. (1), is to maximize the amount of continuous demand d(G) in G that is suitably served (or covered). The decision to be made in the CMCP-1 is where to site a single facility (determine the facility’s x, y coordinates). Demand is considered covered if it is within the service area (V) of the sited facility. U represents the portion(s) of the region that are provided coverage by the facility. This non-linear, non-convex model involves integrating continuous demand over the coverage area U. Since U need not be a single, continuous area, Uk references individual components of U should discontinuities exist. The primary difference between the CMCP-1 and its discrete counterparts is that it allows demand and candidate facility locations to exist everywhere within a region. Thus, in the CMCP-1, the cardinality of the sets of demand (i ˛ I) and facility locations (j ˛ J) is infinite (jIj / N and jJj / N), in contrast to the discrete and finite number of demand and facility locations in the MCLP and discrete demand in the PMCE. In order to solve the maximization problem in Eq. (1), one faces the impossible task of needing to evaluate every candidate facility location (x,y) within the plane for service
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
133
Fig. 1. A non-convex region.
areas V and determining the demand covered. After all locations have been evaluated, the location(s) providing coverage to the greatest amount of demand would represent the optimal site(s). Clearly, the CMCP-1 cannot be approached using standard optimization techniques due to the infinite number of facility and demand locations involved, and the fact that the model can be non-convex, non-linear, and discontinuous! Drezner and Drezner [12] do show how some location problems considering continuous demand can be solved using double integration when the region is convex (square in their case). However, the problem of interest in this paper deals with regions (either convex or non-convex) and their intersection with a facility’s service area, whose integral calculations would be unwieldy, even if approached via quadrature. As a result, it is evident why researchers have focused on special cases of the CMCP-1, like the MCLP and PMCE, that can be discretized, enabling a linear-integer solution approach. Fortunately, it may be possible to exploit geometric properties of a region to be covered in order to establish optimality properties/conditions for the CMCP-1. It is assumed for the remainder of this paper that demand is uniformly distributed in region G (i.e., s ¼ d(G)). As discussed previously, uniformly distributed demand often reflects the realities of many planning problems, like siting warning sirens or radar stations in an urban area, as all areas are equally important to serve. Another assumption is that the service area of the facility (V) is defined as the area within S (as measured in Euclidean distance) of a point centered at (x,y). There are two distinct cases of the CMCP-1. First, the facility’s service area can be completely contained within the region G. Second, the facility’s service area cannot be completely contained within the region. Solution approaches based on these two cases are detailed in Section 5. 4. Geometric representation of regions: the medial axis Before getting into the specifics of the two cases of the CMCP-1, it is necessary to review a geometrical
representation of a region known as the medial axis. Brandt [4] details some important properties of polygons in the plane, and they will be used to define the medial axis. First, let L be defined as the set of all pairs (q,r), where q is a point and r is the maximum radius of a disk centered on the point that is completely contained in region G, or equivalently, the minimum distance from the point to the complement of the region G c (d(q,G c)). In other words, L ¼ fðq; rÞj Cðq; rÞ4Gg, meaning that every location in G is considered to be a center of the largest disk possible that is still completely contained within G. For example, Fig. 2b shows a subset of those points defined in L whose disks are fully contained within the regional boundary given in Fig. 2a. From this definition, it can be shown that the union of the disks centered on these points (and completely contained in G) will exactly regenerate G [4]. In other words,
G [ W Cðq; rÞ: ðq;rÞ˛L
(2)
As can be seen in Fig. 2b, every point in the region will be the center of a disk having at least one point in common with the region’s boundary. Although only a sampling of points are depicted in this figure, it is evident that the edges of their associated disks are already beginning to approximate the general shape of G. As defined above, every point q˛L has an associated radius and it follows that every resulting disk will also be contained within or dominated by a disk of greater or equivalent radius, ðq; rÞ ðq 0 ; r 0 Þ5Cðq; rÞ4Cðq 0 ; r 0 Þ. Given this relationship among points centered in G, the set of points having disks that are maximal is identified as [4]
L [ fðq; rÞ˛L : ðq; rÞ £ ðq 0 ; r 0 Þ˛L0ðq; rÞ [ ðq 0 ; r 0 Þg: Fig. 2c shows a subset of maximally enclosed disks within this region. L* maintains the additional property that each of the associated disks have at least two points in common with the region’s boundary, meaning that their center is equidistant to their two nearest boundary locations. Since all of the disks in L are inferior or equivalent
134
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
and cartographic label placement in map polygons [11]. In the following sections the properties of the medial axis are relied upon in the solution of the CMCP-1.
a
5. Coverage properties for single facility siting The definition of a region’s medial axis is important for the CMCP-1 because every point on the medial axis is the center of a disk of maximal size completely contained within the region. This definition ensures that, under certain conditions, a disk centered on the medial axis will definitely cover as much or more area (uniformly distributed continuous demand in this case) than any other disk within the region. This property of the medial axis is consistent with the objective of the CMCP-1 and can be exploited to solve the CMCP-1. At this point several important observations are necessary relative to the special cases of the CMCP-1.
b
c
d
Fig. 2. Development of the medial axis. (a) A simple, non-convex region, (b) enclosed disks centered at points within the region, (c) maximally enclosed disks centered at points within the region and (d) the union of the centers of all maximally enclosed disks (the medial axis).
to those in L*, the reduced set L* can also be used to exactly regenerate region G as follows:
G[
W Cðq; rÞ: ðq;rÞ˛L
(3)
This relationship establishes that a region can be represented as a locus of points, while still retaining all of the information necessary to exactly regenerate the region. Fig. 2d shows the centers of the maximally enclosed disks for the region in Fig. 2a. This unique set of disk centers, L*, is commonly referred to as the skeleton or medial axis in computational geometry. Brandt [4] formally defines the medial axis as MðGÞ ¼ fq : dr; ðq; rÞ˛L)g. Various techniques exist to generate the medial axis for a region, such as wavefront propagation [2,24] and Voronoi diagram based approaches [6,18,19,28,31]. Given that the medial axis of a polygon is an alternative representation of space, it has found many useful applications in geographical analysis. Some of these uses include point-in-polygon tests [29], statistical analysis of points within a region [33], routing objects through an enclosed space [37], topology generation for scanned maps [15,16],
Observation 1: The smallest service standard needed to ensure complete coverage of a region can be obtained by solving a 1-center problem to identify the center and radius of the smallest disk that can completely enclose the region (see Refs. [10,35]). That is, for service standards greater than or equal to the radius of the smallest enclosing disk, complete regional coverage is always ensured when siting a facility at the 1-center. Observation 2: If S is less than or equal to the maximum radius of the maximally enclosed disks centered on the medial axis (rmax), at least one optimal facility site exists in G. To see this, consider any facility j sited in G. If the service standard, S, for the facility is less than the minimum distance between j and G c, the service area of the facility will be completely contained in G and will provide optimal coverage. Move the facility anywhere in G such that S dðj; G c Þ rmax and the solution is equivalent. Since a disk along the medial axis is maximal with respect to all other locally enclosed disks in G, the site(s) j along the medial axis with the maximum d(j,G c) will correspond to the largest service area that can be completely contained within the region. Observation 3: If the service standard is greater than the radius of the maximum maximally enclosed disk(s) and less than the radius of the smallest enclosing disk associated with the 1-center solution for the region, sites providing optimal coverage will be influenced by Gc. Consider any site j on the medial axis. If the service standard associated with this point is d(j,Gc), any increase to its radius will result in the service area no longer being enclosed in G. A service area, V, not enclosed in G may not necessarily be convex or continuous with respect to VXG. Given this, assessment of whether j maximizes coverage of G depends not only on its ability to provide maximal local coverage, but also on its ability to cover Gc. In this regard, the local properties of the medial axis do not necessarily coincide with optimal siting locations. 6. Solution procedure This section proposes techniques to address the two cases outlined in Observations 2 and 3 for the CMCP-1. An
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
optimal solution to the CMCP-1 can be found for cases where S rmax (Observation 2). This task involves finding the enclosed disk(s) and associated facility location(s) that provides maximal demand coverage of the region. If a facility’s service standard is less than or equal to rmax, the medial axis will always contain an optimal solution. Additionally, when S ¼ rmax where d(j,G c) is locally maximal, the optimal solution will exist only on the medial axis. A solution approach for identifying an optimal location along the medial axis is as follows: CMCP-1 Algorithm I. (Case 1: S rmax) 1 (i) Generate a line-Voronoi diagram for G to obtain medial axis J. (ii) Search J assessing d(j,G c) to find optimal point(s) or line segment(s) where d(j,G c) S. (iii) Apply line-search (or golden section search) technique to identify rmax. In order to find solutions to the CMCP-1 when S > rmax, but less than that the radius of the smallest enclosing disk (Observation 3), a different technique is needed. In this case, any optimal service area having a radius smaller than that of the smallest enclosing disk but greater than that of the maximum maximally enclosed disk must intersect the region’s boundary at multiple points. A reasonable approach for identifying sites that maximize coverage is to search the boundary for these points since an intersection of circles of radius S centered on them will return the center of the optimal service area. The following procedure is one way for addressing this case: CMCP-1 Algorithm II. (Case 2: S > rmax) 2 (i) Given the region’s boundary H, approximate H by the set of points Ha, where a ¼ sampling density. (ii) For each h ˛ Ha, generate a disk of radius S (C(h, S)). (iii) Identify the intersection points of all resulting disk boundaries (denoted vCðh; SÞ) IðHa Þ ¼ X vCðh; SÞ. h˛Ha
(iv) For each intersection point c˛IðHa Þ, generate a disk of radius S ðCðc; SÞÞ. (v) For each resulting disk, compute the area of intersection between the disk and the region G ðCðc; SÞXGÞ, tracking the c generating the disk that covers the most regional demand (cmax). (vi) As a/M (where M is a very large number) the center of the disk Cðcmax ; SÞ is a location maximizing demand coverage.
135
demand (boundary shown in Fig. 1). The objective, Eq. (1), is to cover as much of this region’s demand in locating a single service facility. Implementing the procedures discussed in the previous section, optimal sites for cases where S rmax were computed using a commercial GIS.3 Fig. 3a–d depict sites in G that are optimal for S ¼ 25, 75, 100, and 150 m. For instance, optimal locations for a service standard of 25 m are those sites at least 25 m from any location in the region’s complement (Gc), as established in Observation 2 (Fig. 3a). Given that multiple optima are often present, a GIS is ideal for visualizing such solutions. It is apparent that as S increases, the number of optimal facility sites decreases until only locations along some portions of the medial axis are optimal. In other words, the medial axis will always contain optimal sites for service standards that are less than or equal to rmax. The panels in Fig. 3 also depict optimal portions of the region’s medial axis. Since the region is non-convex, the medial axis comprised line segments and parabolas (caused by concave angles on the region’s boundary). For convex regions, the medial axis consists of only line segments since all angles on the boundary are convex [5]. Searching the medial axis for the sites with the maximum d(j,Gc) results in the identification of a line segment comprising the central trunk of the medial axis (Fig. 4). These points represent the sites at which the largest, fully enclosed service areas may be centered in G. That is, no other sites in G have a larger service area contained in the region than those along the bold line segment shown in Fig. 4. The radii of the maximally enclosed disks, rmax, centered on these points are 163 m. Accordingly, if the service standard for a facility to be sited was 163 m, these locations would represent the set of optimal candidate facility sites, anyone of which optimizes Eq. (1). The union of all of the maximal service areas is also shown in Fig. 4. To establish an upper bound on the largest service standard for which complete regional coverage cannot be achieved, a 1-center problem was solved for region G using the MINIDISC algorithm described in de Berg et al. [10].4 In this case, the maximum distance between the 1-center and the region’s boundary is 306.56 m, representing the radius of the smallest disk capable of containing every point in G (Fig. 5). For service standards less than the upper bound of 306.56 but greater than 163 m, the CMCP-1 Algorithm II (Case 2) outlined in the previous section is necessary to search for a facility site maximizing regional coverage. Fig. 6 illustrates the best sites found for service standards in the range
7. Assessing regional coverage The developed approaches for solving the CMCP-1 are applied to a single non-convex region with uniform
1 Given an n-edge polygon, computational effort associated with obtaining the medial axis is on the order of O(n log n) [19]. Line searches converge linearly with respect to computational effort [32]. 2 Given m ¼ jHaj, then disk intersections may be constructed in O (m2) time. For two polygons involving n vertices and k edge intersections, overlay can be computed in O(n log n þ k) time [27].
3 A true line-Voronoi diagram generally cannot be computed in a GIS (since arcs are generalized as line segments), so an approximation must be used. Here, a suitably dense point-based approximation of the region’s boundary was used and a Delaunay triangulation was computed from these points. The connected centers of the triangles falling completely within the region’s boundary constitute the elements of the medial axis (see Ref. [28]). Approximations of the medial axis, such as that employed here, are shown to converge to the exact medial axis as the number of sampling points along the boundary approaches infinity [5]. 4 1-Center problem programmed in Avenue ArcView GIS version 3.2.
136
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
a
b
c
d
Fig. 3. Optimal candidate facility sites for selected values of S. (a) Optimal sites for S ¼ 25 m, (b) optimal sites for S ¼ 75 m, (c) optimal sites for S ¼ 100 m and (d) optimal sites for S ¼ 150 m.
163 < S < 306.56. In this case, the results show that as the service standard increases above the maximum radius of the maximally enclosed disks (163 m), the center of the selected service area moves in a non-linear trajectory toward the 1-center site. This type of non-linear movement clearly illustrates the influence of the region’s form on the location of the best site.
8. Conclusion This paper developed a methodology for modeling continuous facility location and continuously distributed demand to maximize service coverage in a region when locating a single service facility. The 1-facility continuous maximal covering problem was proposed as
Legend Region Boundary Medial Axis Centers of Maximal Disks Maximal Disks (r = 163 meters) 0
70
140
210
Meters
Fig. 4. The union of the maximum maximally enclosed disks identified on the medial axis.
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
137
Legend Region Boundary 1-Center Location 1-Center Service Area 0
70
140
210
Meters
Fig. 5. 1-center solution and associated service area.
a generalization of its discrete counterpart, the maximal covering location problem. Unlike the MCLP, the CMCP-1 is a non-linear and non-convex model. Given that complete enumeration, linear, or non-linear programming are not feasible options for identifying exact solutions to this type of problem, a means for deriving the optimal solution to the CMCP-1 was proposed. The solution approach is based on the idea that it is not feasible to search the entire region for the facility sites that provide maximal coverage. Instead, an equivalent line representation of a region, known as the medial axis or skeleton, is exploited. The medial axis is an exact representation of a polygon and it is shown that under certain conditions optimal solutions to the CMCP-1 exist on the medial axis. Furthermore, this technique allows for the identification of multiple optima (if present), which often are
desirable for planning purposes. However, in the presence of multiple optima, siting on the medial axis could present additional coverage opportunities given future upgrading of a service facility’s service standard is possible. Consideration of a region’s medial axis definitely provides needed spatial structure to the search for optimal solutions. Dealing with the case where S is greater than the maximum radius of the maximally enclosed disks presents yet another challenge for assessing regional coverage of sited facilities. Here, this issue was addressed through the generation of intersection points derived from circles centered on the region’s boundary and enumerating the regional coverage of each resulting service area. As was described, this procedure need only be performed for service standards greater than the radius of the largest
Fig. 6. Identified facility sites for service standards greater than the maximum radius of the maximally enclosed disks.
138
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139
maximally enclosed disk and less than the radius of the smallest enclosing disk. All of the procedures discussed in this paper are based upon the geometrical properties of a region and are easily assessed using GIS as implemented in the example provided. Hence, an additional benefit of this spatial optimization approach is that it is readily accessible to anyone with a GIS, further highlighting important contributions from GIScience for addressing spatial planning problems. This work represents an initial step in addressing the solution properties/techniques for continuous maximal covering problems, and further work is possible. Future work could involve extending the developed methodology for siting multiple facilities in continuous space to serve continuous demand. Furthermore, this paper only considers facility service areas that are circular in shape. Therefore, exploring facility siting given a broader interpretation of service area potential is of interest. Finally, relaxing the assumption of uniform demand would be a major contribution to this area of research.
Acknowledgements This material is based upon work supported by the National Science Foundation under Grant No. 0518967. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. References [1] Adenso-Diaz B, Rodriguez F. A simple search heuristic for the MCLP: application to the location of ambulance bases in a rural region. Omega 1997;25(2):181–7. [2] Blum H. A transformation for extracting new descriptors of shape. In: Walthen-Dunn W, editor. Models for the perception of speech and visual form. Cambridge: MIT Press; 1967. p. 362– 80. [3] Brady SD, Rosenthal RE. Interactive computer graphical solutions of constrained minimax location problems. AIIE Transactions 1980 September:241–8. [4] Brandt JW. Convergence and continuity criteria for discrete approximations of the continuous planar skeleton. CVGIP: Image Understanding 1994;59(1):116–24. [5] Brandt JW, Algazi VR. Continuous skeleton computation by Voronoi diagram. CVGIP: Image Understanding 1992;55(3):329–38. [6] Chou JJ. Voronoi diagrams for planar shapes. IEEE Computer Graphics and Applications 1995 March:52–9. [7] Church RL, ReVelle CS. The maximal covering location problem. Papers – Regional Science Association 1974;32:101–18. [8] Church RL. The planar maximal covering location problem. Journal of Regional Science 1984;24(2):185–201. [9] Current J, O’Kelly M. Locating emergency warning sirens. Decision Sciences 1992;23:221–34. [10] de Berg M, van Kreveld M, Overmars M, Schwarzkopf O. Computational geometry: algorithms and applications. 2nd ed. New York: Springer-Verlag; 2000. [11] Do¨rschlag D, Petzold I, Plu¨mer L. Placing objects automatically in areas of maps. In: Proceedings of the Durban ICC 2003, South Africa. [12] Drezner T, Drezner Z. Replacing continuous demand with discrete demand in a competitive location model. Naval Research Logistics 1997;44:81–95. [13] Gavirneni S, Clark L, Pataki G. Schlumberger optimizes receiver location for automated meter reading. Interfaces 2004;34(3): 208–14. [14] Gleason JM. A set covering approach to bus stop location. Omega 1975;3:605–8.
[15] Gold C, Snoeyink J. A one-step crust and skeleton extraction algorithm. Algorithimica 2001;30:144–63. [16] Gold CM, Nantel J, Yang W. Outside-in: an alternative approach to forest map digitizing. International Journal of Geographical Information Systems 1996;10(3):291–310. [17] Kalvenes J, Kennington J, Olinick E. Hierarchical cellular network design with channel allocation. European Journal of Operational Research 2004;160(1):3–18. [18] Kirkpatrick DG. Efficient computation of continuous skeletons. In: Proceedings of the 20th annual symposium on foundations of computer science; 1979. p. 18–27. [19] Lee DT. Medial axis transformation of a planar shape. IEEE Transactions on Pattern Analysis and Machine Intelligence 1982;PAMI-4(4): 363–9. [20] Love RF. A computational procedure for optimally locating a facility with respect to several rectangular regions. Journal of Regional Science 1972;12(2):233–42. [21] Mehrez A. A note on the linear integer formulation of the maximal covering location problem with facility placement on the entire plane. Journal of Regional Science 1983;23(4):553–5. [22] Mehrez A, Stulman A. The maximal covering location problem with facility placement on the entire plane. Journal of Regional Science 1982;22(3):361–5. [23] Mehrez A, Stulman A. An extended continuous maximal covering location problem with facility placement. Computers & Operations Research 1984;11(1):19–23. [24] Montanari U. Continuous skeletons from digitized images. Journal of the Association for Computing Machinery 1969;16(4): 534–49. [25] Murray AT, O’Kelly ME, Church RL. Regional service coverage modeling. Computers & Operations Research 2008;35(2):339– 55. [26] O’Kelly ME, Murray AT. A lattice covering model for evaluating existing service facilities. Papers in Regional Science 2004;83: 565–80. [27] O’Rourke J. Computational geometry in C. 2nd ed. Cambridge: Cambridge University Press; 1998. [28] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts & applications of Voronoi diagrams. 2nd ed. New York: John Wiley & Sons, Ltd.; 2000. [29] Pfaltz JL, Rosenfeld A. Computer representation of planar regions by their skeletons. Communications of the ACM 1967;10(2): 119–25. [30] Plastria F. Continuous covering location problems. In: Drezner Z, Hamacher H, editors. Facility location: applications and theory. New York: Springer; 2002. [31] Preparata FP. The medial axis of a simple polygon. In: Goos G, Hartmanis J, editors. Mathematical foundations of computer science, lecture notes in computer science, vol. 53. New York: SpringerVerlag; 1977. p. 443–50. [32] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C: the art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [33] Sadahiro Y, Takami K. A statistical method for analyzing the relative location of points in a bounded region. Geographical Analysis 2001; 33(1):40–57. [34] Suzuki A, Okabe A. Using Voronoi diagrams. In: Drezner Z, editor. Facility location: a survey of applications and methods. New York: Springer-Verlag; 1995. p. 103–18. [35] Suzuki A, Drezner Z. The p-center location problem in an area. Location Science 1996;4(1/2):69–82. [36] Toregas C, Swain R, ReVelle C, Bergman L. The location of emergency service facilities. Operations Research 1971;19: 1363–73. [37] Xu Y, Mattikalli R, Khosla P. Determining two circumscribing discs for a polygon. Computer-Aided Design 1995;27(3):223–34.
Tim Matisziw received his B.A. and M.A. in Geography from the University of Missouri – Columbia and a Ph.D. in Geography from The Ohio State University. Dr. Matisziw is currently Assistant Professor, Department of Geography and Department of Civil & Environmental Engineering at the University of Missouri – Columbia. His current research addresses modeling issues in transportation, location science, and spatial optimization and has appeared in journals such as the European Journal of Operational Research, Computers & Operations Research, Decision Support Systems, International Regional Science Review, Computers, Environment and Urban Systems, Journal of Geographical Systems, Transportation Research C, and Regional Studies.
T.C. Matisziw, A.T. Murray / Socio-Economic Planning Sciences 43 (2009) 131–139 Alan Murray obtained a B.S. in mathematical sciences, an M.A. in statistics and applied probability and a Ph.D. in geography, all from the University of California at Santa Barbara. Dr. Murray is Professor, School of Geographical Sciences, at Arizona State University. His research and teaching interests are in spatial optimization modeling, geographic information science, urban/regional planning and development, and transportation. He has
139
published on a range of technical and application oriented topics in journals such as the Operations Research, Discrete Applied Mathematics, International Journal of Geographic Information Sciences, Geographical Analysis, Journal of Geographical Systems, Journal of Urban Planning and Development, European Journal of Operational Research, Computers and Operations Research, Transportation Research, and Urban Studies, among others.