# A Decomposition Heuristic for the Maximal Covering Location Problem

Authors: Edson Luiz França Senne; Marcos Antonio Pereira; Luiz Antonio Nogueira Lorena
Journal: Advances in Operations Research
 (2010)

---

## Abstract

This paper proposes a cluster partitioning technique to calculate improved upper bounds to the optimal solution of maximal covering location problems. Given a covering distance, a graph is built considering as vertices the potential facility locations, and with an edge connecting each pair of facilities that attend a same client. Coupling constraints, corresponding to some edges of this graph, are identified and relaxed in the Lagrangean way, resulting in disconnected subgraphs representing smaller subproblems that are computationally easier to solve by exact methods. The proposed technique is compared to the classical approach, using real data and instances from the available literature.


---

## Body

## 1. Introduction

The covering class of facility location problems deals with the maximum distance between any client and the facility designed to attend an associated demand. These problems are known as covering problems and the maximum service distance is known as covering distance. The Set Covering Problem [1] determines the minimal number of facilities that are necessary to attend all clients, for a given covering distance. Due to formulation restrictions, this model does not consider the individual demand of each client. In addition, the number of needed facilities can be large, incurring high fixed installation costs. An alternative formulation considers the installation of a limited number of facilities, even if this amount is unable to attend the total demand. In this formulation, the condition that all clients must be served is relaxed and the objective is changed to locate p facilities such that the most part of the existing demand can be attended, for a given covering distance. This model corresponds to the Maximal Covering Location Problem (MCLP).Covering models are often found in problems of public organizations for the location of emergency services. Early techniques for solving the MCLP tried to obtain integer solutions from the linear relaxation equivalent of the model proposed by Church and ReVelle [2]. This pioneer work formalizes the MCLP and presents a greedy heuristic based on vertices exchange.MCLP applications range from emergency services [3, 4], hierarchical health services [5], air pollution control [6], to congested systems [7–9]. Solution methods for the MCLP include the linear programming relaxation [2], greedy heuristics [10], and Lagrangean relaxation [11]. Lorena and Pereira [12] report results obtained with a Lagrangean/surrogate heuristic using a subgradient optimization method, in complement to the dissociated Lagrangean and surrogate heuristics presented in [13]. Arakaki and Lorena [14] present a constructive genetic algorithm to solve real case instances with up to 500 vertices. Surveys can be found in [15–18].In this paper is presented a cluster relaxation technique to solve large-scale maximal covering location problems. The proposed approach requires the identification of a graph related to a set of constraints. If some of these constraints are relaxed, this graph can be partitioned into subgraphs (clusters), corresponding to smaller problems that can be solved independently.This paper is organized as follows. Section2 presents a decomposition approach to obtain improved upper bounds to the optimal solution of maximal covering location problems. Section 3 describes computational results. Some conclusions are given in Section 4.


## 2. The Proposed Approach for the MCLP

The MCLP was formulated in [2] as the following 0-1 linear programming problem:(2.1)(MCLP)v(MCLP)=Max∑i∈Nwixi(2.2)subjectto∑j∈Siyj≥xi,∀i∈N,(2.3)∑j∈Myj=p,(2.4)xi∈{0,1},∀i∈N,(2.5)yj∈{0,1},∀j∈M,
where(i)
M={1,2,…,m} is the set of potential facility locations;(ii)
N={1,2,…,n} is the set of clients to be covered;(iii)
D=[dij] is the Euclidean distance matrix between each pair of nodes i∈N and j∈M;(iv)
U is the covering distance;(v)
Si={j∈M∣dij≤U} is the set of facilities that can attend each client i∈N;(vi)
wi is the demand (a positive integer value) for each client i∈N;(vii)
p is the number of facilities to be located;(viii)
xi is a decision variable, with xi=1 if the demand of the client i is covered, and xi=0, otherwise;(ix)
yj is a decision variable, with yj=1 if a facility was installed at the location j, and yj=0, otherwise.The objective function maximizes the covered demand. Constraints (2.2) state that a client will be covered if there is at least one facility located within the covering distance. Constraint (2.3) limits to exactly p the number of located facilities and (2.4) and (2.5) express the binary conditions.The traditional Lagrangean relaxation approach [11] relaxes the set of constraints (2.2) with a vector μ of multipliers μi≥0, i∈N, obtaining(2.6)(LRμ)v(LRμ)=Max∑i∈N(wi-μi)xi+∑i∈N∑j∈Siμiyjsubject to (2.3), (2.4), and (2.5).It is easy to see by the integrality property thatv(LRμ)≥v(MCLP) and the Lagrangean bound cannot be better than the linear relaxation of (MCLP).In this paper, a decomposition approach based on the Lagrangean relaxation with clusters (LagClus) of Ribeiro and Lorena [19–22] is presented. LagClus is a stronger relaxation that can be useful for several theoretical and practical large-scale problems. The first application of the LagClus was performed on point-feature instances. Later, Ribeiro and Lorena applied this relaxation on pallet loading instances obtaining good results. Besides, the authors proposed a column generation for that problem using this cluster relaxation idea. Another interesting application was performed on wood pulp stowage context. This problem consists of arranging items into holds of dedicated maritime international ships. Recently, [23] applied the LagClus to uncapacitated facility location instances providing better bounds than the ones presented in the literature for a set of difficult instances.Consider the MCLP instance represented in Figure1, where the dots correspond to the clients to be covered (N={1,…,12}) and the small squares correspond to the potential facility locations (M={1,…,7}). In this figure, the values in parenthesis are the demand values wi, for all i∈N. For the chosen covering distance, the sets Si, for all i∈N, are defined as follows:(2.7)S1={1},S4={4},S2={1},S5={3,4},S3={1,2},S6={1,2,3,5},S7={2,6},S10={4},S8={2,5,6},S11={5,7},S9={3,4,7},S12={7}.Figure 1
A MCLP instance.Assuming the number of facilities to be installed asp = 3, this instance can be formulated as(2.8)v(MCLP)=Max3x1+x2+2x3+4x4+5x5+2x6+3x7+7x8+x9+x10+2x11+2x12(2.9)subjecttoy1≥x1(2.10)y1≥x2,(2.11)y1+y2≥x3,(2.12)y4≥x4,(2.13)y3+y4≥x5,(2.14)y1+y2+y3+y5≥x6,(2.15)y2+y6≥x7,(2.16)y2+y5+y6≥x8,(2.17)y3+y4+y7≥x9,(2.18)y4≥x10,(2.19)y5+y7≥x11,(2.20)y7≥x12,(2.21)y1+y2+y3+y4+y5+y6+y7=3,(2.22)xi∈{0,1},∀i∈N,(2.23)yj∈{0,1},∀j∈M.
The proposed approach considers the MCLP as a covering graph, which is defined as follows.LetG(M,A) be a graph where M is the set of vertices corresponding to the potential facility locations and A={(p,q):pand q∈Si, i=1,…,N} is the set of edges. So, in a covering graph there exists an edge (p,q) connecting two potential facility locations, if p and q share, at least, one covered client. Figure 2 shows the covering graph associated with the above MCLP instance.Figure 2
A covering graph.It is easy to note that the edges in a covering graph are related to the set of constraints (2.2) in the MCLP formulation. For example, the edges shown in Figure 2 correspond to the constraints (2.11), (2.13)–(2.17), and (2.19).Now, consider that a covering graph is partitioned in some way. Figure3(a) shows a possible partition. If the edges (1,3), (2,3), (3,5), and (5,7) are removed from the graph, two subgraphs are obtained, as shown in Figure 3(b).Partitioning a covering graph.
(a)
A possible partition(b)
The corresponding subgraphsThis partition corresponds to relax in the Lagrangean way the constraints (2.14) and (2.19) of the MCLP formulation, using Lagrangean multipliers λ6 and λ11, respectively. If constraint (2.3) is also relaxed with Lagrangean multiplier μ, the relaxed problem will be(2.24)v(MCLPR)=Max3x1+x2+2x3+4x4+5x5+2x6+3x7+7x8+x9+x10+2x11+2x12+λ6(y1+y2+y3+y5-x6)+λ11(y5+y7-x11)+μ(y1+y2+y3+y4+y5+y6+y7-3)
subject to (2.9)–(2.13), (2.15)–(2.18), (2.20), (2.22), (2.23), and(2.25)λ6,λ11≥0,μ∈R.
The objective function can now be rewritten as(2.26)v(MCLPR)=Max3x1+x2+2x3+4x4+5x5+(2-λ6)x6+3x7+7x8+x9+x10+(2-λ11)x11+2x12+(λ6+μ)y1+(λ6+μ)y2+(λ6+μ)y3+μy4+(λ6+λ11+μ)y5+μy6+(λ11+μ)y7-3μ.
Then, the problem can be decomposed in two subproblems:(2.27)v(MCLPR)=v(SP1)+v(SP2)+(2-λ6)x6+(2-λ11)x11-3μ,
where(2.28)v(SP1)=Max3x1+x2+2x3+3x7+7x8+(λ6+μ)y1+(λ6+μ)y2+(λ6+λ11+μ)y5+μy6
subject to (2.9)–(2.11), (2.15), (2.16), (2.22), (2.23), and(2.29)λ6,λ11≥0,μ∈R,
and
(2.30)v(SP2)=Max4x4+5x5+x9+x10+2x12+(λ6+μ)y3+μy4+(λ11+μ)y7
subject to (2.12), (2.13), (2.17), (2.18), (2.20), (2.22), (2.23), and(2.31)λ6,λ11≥0,μ∈R.Note that these subproblems correspond to the following clusters (which are associated with the subgraphs of the covering graph):(i)
Cluster 1:M1={1,2,3,7,8} and N1={1,2,5,6};(ii)
Cluster 2:M2={4,5,6,9,10,11,12} and N2={3,4,7}.The resulting Lagrangean relaxation does not have the integrality property, being stronger than (LRμ) of Galvão and ReVelle [11]. As the clusters are smaller than the original covering graph, exact methods can be employed to solve each corresponding subproblem, obtaining better quality bounds in shorter computational times.For the above example, one can apply a subgradient optimization method in order to determine the values of the dual variablesλ6, λ11, and μ. At each iteration of this method, the values of the dual variables are used to update the coefficients of the objective function of each subproblem, which can be solved by exact methods.It is interesting to observe that, due to the relaxation of constraints (2.14) and (2.19), the variables x6 and x11 disappeared from the formulations of the subproblems. According to the objective function (2.27) their values will be set as 0 or 1, depending on their respective coefficient at each iteration:(2.32)xi={1,ifwi-λi>0,0,otherwise
for all i∈I={i∈N∣constraint (2.2) containing variable xi is relaxed}.Therefore, the proposed decomposition approach can be established in the following steps.(a)
Build a covering graphG(M,A) corresponding to the MCLP.(b)
Apply a graph partitioning heuristic to divide the covering graphG into k clusters. The MCLP can be written through the objective function defined in (2.1) subject to (2.2)–(2.5) where the constraints (2.2) are now divided into two groups: one with constraints corresponding to intracluster edges and other formed by constraints that correspond to edges connecting the clusters.(c)
Using distinct nonnegative multipliers, relax in the Lagrangean way the constraints corresponding to the edges connecting the clusters (defining the setI) and also relax constraint (2.3).(d)
The resulting Lagrangean relaxation is decomposed intok subproblems.(e)
Apply the standard subgradient method in order to optimize the dual variablesλ and μ.The subgradient method used in the step (e) can be written as follows.Set, initially,(2.33)λi={wi,ifi∈I,0,otherwise,μ=0,θ=2,LB=-∞,UB=+∞,While (the stop conditions are not satisfied) do the followingSolve subproblems SPk forx andy, using the current values ofλ andμ.Calculate(2.34)v(MCLPR)=∑kv(SPk)+∑i∈Imax{0,(wi-λi)}-pμ.UpdateUB=min{UB,v(MCLPR)}.If∑j∈Nyj=p (a feasible solution for MCLP is found), then do the following.Calculatev(MCLP) using the obtained values for x and y.UpdateLB=max{LB,v(MCLP)}.Calculate(2.35)giλ=xi-∑j∈Siyj,i∈I,gμ=p-∑j∈Nyj.Update the step sizeθ.Update(2.36)λi=min{0,λi+θgiλ},i∈I,μ=μ+θgμ.End-While.The step sizeθ used in this algorithm is the one proposed in [24], beginning with θ=2 and halving it whenever the upper bound does not decrease for a certain number of successive iterations. The stopping tests used are the following: (a)
θ≤0.005, or(b)
UB-LB<1, or(c)
the subgradient vectorg=[giλ,gμ]=0.


## 3. Computational Results

The LagClus algorithm was coded in C and the tests were conducted on a notebook with Intel Core 2 Duo 2.0 GHz processor and 2.0 GB RAM, running Windows XP (Service Pack 3), and ILOG CPLEX 10.1.1 [25]. The data were obtained from TSPLIB PCB3038 [26] and real case instances for facility location problems in São José dos Campos city, Brazil (available for download at http://www.lac.inpe.br/~lorena/instancias.html).For the graph partitioning task was used the well-known METIS heuristic for graph partitioning problems [27], with default values. Given a covering graph Gand a predefined number kof clusters, METIS divides the graph in kclusters minimizing the number of edges with terminations in different clusters.The results obtained are shown in Tables1 to 3. These tables use the following legend in the columns(i)
k: number of clusters;(ii)
n: number of potential facilities locations and clients to be covered;(iii)
p: number of facilities to be installed;(iv)
Optimal: optimal solution of the corresponding MCLP obtained by CPLEX;(v)
GapLP: linear relaxation gap provided by CPLEX (in percentage);(vi)
CPU: computing time of the corresponding optimal solution obtained by CPLEX (in seconds);(vii)
Cuts: number of relaxed constraints;(viii)
Gap = 100%  × (UB − Optimal)/Optimal;(ix)
CPUSeq: sum of the computing times for every cluster, at each iteration (in seconds);(x)
CPUPar: the largest computing time for a cluster, at each iteration (in seconds).Table 1
Computation times for SJC instances,U=150.
k = 10k = 50npOptimalGapLPCPUCutsGapCPUSeqCPUParCutsGapCPUSeqCPUPar3242073020.0000.0152960.0002.5432.53714490.0009.8759.7393091270.0270.0470.02724.65024.6350.02744.59641.10140104430.1560.1880.10825.98525.9690.15745.12242.35350113970.1800.3910.13824.45224.4220.19543.74843.45960119910.1840.2350.02444.51444.4210.22148.04947.90980121520.0000.0310.0008.8768.8160.003109.129108.115108121520.0000.0160.0001.5951.5952.97726.53326.48750040133400.0000.0471080.0003.4533.3618040.00017.18611.49950147730.0140.0470.0004.9384.5140.01459.01936.29360159190.0480.0630.0008.2337.2430.04857.15734.73770169080.0000.0310.0003.7233.3700.00022.89114.42180177490.0000.0150.0005.4064.7660.00026.68616.697100189120.0980.1090.00010.2767.1710.05662.07137.748130196640.0410.2970.01530.82724.5880.04169.93443.373167197060.0050.0470.00314.60014.2350.00546.07835.41481880233250.0550.1401660.00345.56421.88016490.06185.92245.81990244550.1230.2660.04156.38824.7470.14387.79747.001100254350.1270.3440.01287.27934.4810.14096.12452.060120269820.0840.2970.01569.65831.3680.062105.54754.658140280020.1400.3590.09552.96626.2710.128121.12763.713160286990.1280.3910.10758.45324.9040.12696.01750.828200291530.0180.2340.01161.53128.3010.039253.908135.048273291680.0000.0310.0003.3432.5450.55446.76637.178The values marked with an asterisk in Table3 denote instances where CPLEX failed to produce an optimal solution within the time limit of 20000 seconds. The presented figures are suboptimal values.From these results one can observe that the smaller the number of clusters is, the better are the upper bounds obtained (smaller are the gaps). On the other hand, as the number of clusters increases, the computational effort for solving the subproblems is reduced.In the results shown that fork=10 (SJC instances) and k=5 (TSPLIB PCB3038 instance), the gaps obtained for 92.3% of the instances (48 out of 52) are equal or smaller (values presented in bold face) than the gaps obtained by the linear relaxation, which demonstrates the effectiveness of the decomposition heuristic. However, improved gaps could be obtained by reducing the number of clusters, at the cost of larger computational times.Therefore, as shown in Tables1, 2, and 3, the bounds obtained by the proposed approach are better than the ones produced by the linear relaxation and, consequently, by any Lagrangean relaxation as the ones presented in [11–13]. Nevertheless, comparing the values of the CPU times presented by this method and those presented by CPLEX, it is clear that the proposed decomposition approach is appropriate only for large-scale problems.Table 2
Computation times for SJC instances,U=200.
k=10k=50npOptimalGapLPCPUCutsGapCPUSeqCPUParCutsGapCPUSeqCPUPar3242096700.3340.1729800.24319.29319.29330080.34732.94332.91130117370.0870.4840.06028.94328.9430.09469.95969.95940121510.0050.0940.00831.06631.0660.13877.87277.61250121520.0000.0150.0009.9269.9260.00061.45661.39560121520.0000.0470.0004.5754.5750.00014.37614.37680121520.0000.0160.0003.6703.6700.00033.93633.936108121520.0000.0310.24811.34311.3430.001897.830893.40550040170770.4530.2036570.38724.66820.66926250.46955.00151.04850183610.0140.1090.00339.10932.2480.02567.62662.59660191530.0350.0630.00552.63935.3630.11285.37476.57870195510.1101.0780.06943.94629.8170.17076.67168.72180197030.0130.1560.00835.49527.9270.150102.05695.253100197070.0000.0780.00016.62416.5010.00189.85886.698130197070.0000.0470.0001.9861.8640.00126.54625.809167197070.0000.0160.01622.37922.2250.85922.31422.13381880279450.0700.2038400.06957.83527.42349100.121147.155115.60590285190.1381.1410.071114.14545.5360.177128.57499.585100289100.1031.3910.03688.88533.1530.175101.87580.758120291650.0021.2340.00255.71031.1800.117141.246115.434140291680.0000.1250.00011.6438.9400.021171.961143.737160291680.0000.0620.0009.7387.5980.87839.34336.244200291680.0000.0320.0005.7624.6100.84737.20534.871273291680.0000.0310.20724.69818.5052.90425.28223.772Table 3
Computation times for TSPLIB PCB3038 instance,U=400.
k=5k=10npOptimalGapLPCPUCutsGapCPUSeqCPUParCutsGapCPUSeqCPUPar3038171253200.368802.3901655790.205843.838235.5412913630.470582.528223.245181300040.51710265.0160.372817.076283.4000.712634.402243.74719134262*0.60520000.0490.3821483.237598.6530.793576.821222.08720139028*0.69820000.1560.5001712.078798.9110.973628.288236.76721141279*0.85320000.0940.6543117.1741448.7301.128646.765243.30222143809*1.19620000.1230.9926656.2673094.4101.598615.783231.525Comparing the values ofCPUSeq and CPUPar for each value of k, it is also possible to note that CPUPar values are becoming relatively smaller than CPUSeq values as the size of the instances increases, indicating that the proposed decomposition approach can substantially reduce the time for solving large instances of MCLP in parallel (or multicore CPU) computers.


## 4. Conclusions

This paper presents a decomposition approach based on cluster partitioning to calculate improved upper bounds to the optimal solution of maximal covering location problems. The partitioning is based on the covering graph of potential facility locations that attend a same client. The corresponding coupling constraints are identified and some of them are relaxed in the Lagrangean way, resulting in subproblems that can be solved independently. Each subproblem represents a cluster smaller than the original problem and can be solved by exact methods in smaller computational times. Computational tests using real data and instances from the available literature were conducted and confirmed the effectiveness of the proposed approach.An important characteristic of large-scale problems addressed by the proposed approach is the tradeoff between gap values and CPU times. Depending on the application, one may choose to sacrifice the quality of the bounds (increasing the number of clusters) in order to obtain shorter processing times. On the other hand, if quality is the issue, the processing times needed to solve instances with only a few clusters may be longer.In this study, the number of clusters was fixeda priori and then the number of intercluster edges was minimized by the partitioning algorithm. In this way, for a chosen k, the number of relaxed constraints is minimized and, consequently, better bounds are obtained. The partitioning is clearly crucial and other choices could be considered, for example, minimize the maximal size of the clusters, so as to obtain the smallest possible subproblems. This has yet to be investigated.The heuristic presented in this article can be used in a branch-and-bound exact method. As, in general, the upper bounds obtained with this heuristic are better than those obtained by the linear relaxation, one would expect many more nodes be pruned, with a possibly significant size reduction of the search tree.Advances in applied mathematics and computer science have resulted in high-performance tools for mathematical programming, allowing tough optimization problems to be solved. However, as the problem size increases, the computational time may grow excessively, making the problem intractable even for the most efficient tools. In such a case an approach in which a large-scale problem is divided into a number of smaller-scale subproblems can be a nice solution possibility.



---

*Source: 120756-2010-05-03.xml*