A Decomposition Heuristic for the Maximal Covering Location Problem

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.

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. Section 2 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:(MCLP)𝑣(MCLP)=Max_𝑖∈𝑁𝑤_𝑖𝑥_𝑖subjectto_{𝑗∈𝑆𝑖}𝑦_𝑗≥𝑥_𝑖,∀𝑖∈𝑁,_𝑗∈𝑀𝑦_𝑗𝑥=𝑝,_𝑖𝑦∈{0,1},∀𝑖∈𝑁,_𝑗∈{0,1},∀𝑗∈𝑀, (2.1)(2.2)(2.3)(2.4)(2.5) where(i)𝑀={1,2,…,𝑚} is the set of potential facility locations;(ii)𝑁={1,2,…,𝑛} is the set of clients to be covered;(iii)𝐷=[𝑑_𝑖𝑗] is the Euclidean distance matrix between each pair of nodes 𝑖∈𝑁 and 𝑗∈𝑀;(iv)𝑈 is the covering distance;(v)𝑆_𝑖={𝑗∈𝑀∣𝑑_𝑖𝑗≤𝑈} is the set of facilities that can attend each client 𝑖∈𝑁;(vi)𝑤_𝑖 is the demand (a positive integer value) for each client 𝑖∈𝑁;(vii)𝑝 is the number of facilities to be located;(viii)𝑥_𝑖 is a decision variable, with 𝑥_𝑖=1 if the demand of the client 𝑖 is covered, and 𝑥_𝑖=0, otherwise;(ix)𝑦_𝑗 is a decision variable, with 𝑦_𝑗=1 if a facility was installed at the location 𝑗, and 𝑦_𝑗=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 𝜇_𝑖≥0, 𝑖∈𝑁, obtaining(LR_𝜇)𝑣(LR_𝜇)=Max_𝑖∈𝑁(𝑤_𝑖−𝜇_𝑖)𝑥_𝑖+_{𝑖∈𝑁𝑗∈𝑆𝑖}𝜇_𝑖𝑦_𝑗 (2.6)subject to (2.3), (2.4), and (2.5).

It is easy to see by the integrality property that 𝑣(LR_𝜇)≥𝑣(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 Figure 1, where the dots correspond to the clients to be covered (𝑁={1,…,12}) and the small squares correspond to the potential facility locations (𝑀={1,…,7}). In this figure, the values in parenthesis are the demand values 𝑤_𝑖, for all 𝑖∈𝑁. For the chosen covering distance, the sets 𝑆_𝑖, for all 𝑖∈𝑁, are defined as follows:𝑆₁={1},𝑆₄𝑆={4},₂={1},𝑆₅𝑆={3,4},₃={1,2},𝑆₆𝑆={1,2,3,5},⁷={2,6},𝑆₁₀𝑆={4},₈={2,5,6},𝑆₁₁𝑆={5,7},₉={3,4,7},𝑆₁₂={7}. (2.7)

Figure 1 

A MCLP instance.

Assuming the number of facilities to be installed as p = 3, this instance can be formulated as𝑣(MCLP)=Max3𝑥₁+𝑥₂+2𝑥₃+4𝑥₄+5𝑥₅+2𝑥₆+3𝑥₇+7𝑥₈+𝑥₉+𝑥₁₀+2𝑥₁₁+2𝑥₁₂subjectto𝑦₁≥𝑥₁𝑦₁≥𝑥₂𝑦₁,+𝑦₂≥𝑥₃𝑦₄,≥𝑥₄𝑦₃,+𝑦₄≥𝑥₅𝑦₁,+𝑦₂+𝑦₃+𝑦₅≥𝑥₆𝑦₂,+𝑦₆≥𝑥₇𝑦₂,+𝑦₅+𝑦₆≥𝑥₈𝑦₃,+𝑦₄+𝑦₇≥𝑥₉𝑦₄,≥𝑥₁₀𝑦₅,+𝑦₇≥𝑥₁₁𝑦₇,≥𝑥₁₂𝑦₁,+𝑦₂+𝑦₃+𝑦₄+𝑦₅+𝑦₆+𝑦₇𝑥=3,_𝑖𝑦∈{0,1},∀𝑖∈𝑁,_𝑗∈{0,1},∀𝑗∈𝑀. (2.8)(2.9)(2.10)(2.11)(2.12)(2.13)(2.14)(2.15)(2.16)(2.17)(2.18)(2.19)(2.20)(2.21)(2.22)(2.23) The proposed approach considers the MCLP as a covering graph, which is defined as follows.

Let 𝐺(𝑀,𝐴) be a graph where M is the set of vertices corresponding to the potential facility locations and 𝐴={(𝑝,𝑞)∶𝑝and 𝑞∈𝑆_𝑖, 𝑖=1,…,𝑁} is the set of edges. So, in a covering graph there exists an edge (𝑝,𝑞) connecting two potential facility locations, if 𝑝 and 𝑞 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. Figure 3(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).

(a) A possible partition

(b) The corresponding subgraphs

Figure 3 

Partitioning a covering graph.

This partition corresponds to relax in the Lagrangean way the constraints (2.14) and (2.19) of the MCLP formulation, using Lagrangean multipliers 𝜆₆ and 𝜆₁₁, respectively. If constraint (2.3) is also relaxed with Lagrangean multiplier 𝜇, the relaxed problem will be𝑣(MCLP_𝑅)=Max3𝑥₁+𝑥₂+2𝑥₃+4𝑥₄+5𝑥₅+2𝑥₆+3𝑥₇+7𝑥₈+𝑥₉+𝑥₁₀+2𝑥₁₁+2𝑥₁₂+𝜆₆(𝑦₁+𝑦₂+𝑦₃+𝑦₅−𝑥₆)+𝜆₁₁(𝑦₅+𝑦₇−𝑥₁₁)+𝜇(𝑦₁+𝑦₂+𝑦₃+𝑦₄+𝑦₅+𝑦₆+𝑦₇−3) (2.24) subject to (2.9)–(2.13), (2.15)–(2.18), (2.20), (2.22), (2.23), and𝜆₆,𝜆₁₁≥0,𝜇∈𝑅. (2.25) The objective function can now be rewritten as𝑣(MCLP_𝑅)=Max3𝑥₁+𝑥₂+2𝑥₃+4𝑥₄+5𝑥₅+(2−𝜆₆)𝑥₆+3𝑥₇+7𝑥₈+𝑥₉+𝑥₁₀+(2−𝜆₁₁)𝑥₁₁+2𝑥₁₂+(𝜆₆+𝜇)𝑦₁+(𝜆₆+𝜇)𝑦₂+(𝜆₆+𝜇)𝑦₃+𝜇𝑦₄+(𝜆₆+𝜆₁₁+𝜇)𝑦₅+𝜇𝑦₆+(𝜆₁₁+𝜇)𝑦₇−3𝜇. (2.26) Then, the problem can be decomposed in two subproblems:𝑣(MCLP_𝑅)=𝑣(SP₁)+𝑣(SP₂)+(2−𝜆₆)𝑥₆+(2−𝜆₁₁)𝑥₁₁−3𝜇, (2.27) where𝑣(SP₁)=Max3𝑥₁+𝑥₂+2𝑥₃+3𝑥₇+7𝑥₈+(𝜆₆+𝜇)𝑦₁+(𝜆₆+𝜇)𝑦₂+(𝜆₆+𝜆₁₁+𝜇)𝑦₅+𝜇𝑦₆ (2.28) subject to (2.9)–(2.11), (2.15), (2.16), (2.22), (2.23), and𝜆₆,𝜆₁₁≥0,𝜇∈𝑅, (2.29) and 𝑣(SP₂)=Max4𝑥₄+5𝑥₅+𝑥₉+𝑥₁₀+2𝑥₁₂+(𝜆₆+𝜇)𝑦₃+𝜇𝑦₄+(𝜆₁₁+𝜇)𝑦₇ (2.30) subject to (2.12), (2.13), (2.17), (2.18), (2.20), (2.22), (2.23), and𝜆₆,𝜆₁₁≥0,𝜇∈𝑅. (2.31)

Note that these subproblems correspond to the following clusters (which are associated with the subgraphs of the covering graph):(i)Cluster 1: 𝑀₁={1,2,3,7,8} and 𝑁₁={1,2,5,6};(ii)Cluster 2: 𝑀₂={4,5,6,9,10,11,12} and 𝑁₂={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 𝜆₆, 𝜆₁₁, 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 𝑥₆ and 𝑥₁₁ 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:𝑥_𝑖=1,if𝑤_𝑖−𝜆_𝑖>0,0,otherwise (2.32) for all 𝑖∈𝐼={𝑖∈𝑁∣constraint (2.2) containing variable 𝑥_𝑖 is relaxed}.

Therefore, the proposed decomposition approach can be established in the following steps.(a)Build a covering graph 𝐺(𝑀,𝐴) corresponding to the MCLP.(b)Apply a graph partitioning heuristic to divide the covering graph 𝐺 into 𝑘 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 set 𝐼) and also relax constraint (2.3).(d)The resulting Lagrangean relaxation is decomposed into 𝑘 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, 𝜆_𝑖=𝑤_𝑖,if𝑖∈𝐼,0,otherwise,𝜇=0,𝜃=2,LB=−∞,UB=+∞, (2.33) While (the stop conditions are not satisfied) do the following Solve subproblems SP_𝑘 for𝑥 and𝑦, using the current values of𝜆 and𝜇. Calculate 𝑣(MCLP_𝑅)=_𝑘𝑣(SP_𝑘)+_𝑖∈𝐼max{0,(𝑤_𝑖−𝜆_𝑖)}−𝑝𝜇. (2.34) Update UB=min{UB,𝑣(MCLP_𝑅)}. If ∑_𝑗∈𝑁𝑦_𝑗=𝑝 (a feasible solution for MCLP is found), then do the following.   Calculate 𝑣(MCLP) using the obtained values for 𝑥 and 𝑦.   Update LB=max{LB,𝑣(MCLP)}. Calculate 𝑔^𝜆_𝑖=𝑥_𝑖−_{𝑗∈𝑆𝑖}𝑦_𝑗𝑔,𝑖∈𝐼,_𝜇=𝑝−_𝑗∈𝑁𝑦_𝑗. (2.35) Update the step size 𝜃. Update 𝜆_𝑖=min{0,𝜆_𝑖+𝜃𝑔^𝜆_𝑖},𝑖∈𝐼,𝜇=𝜇+𝜃𝑔^𝜇. (2.36)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 vector 𝑔=[𝑔^𝜆_𝑖,𝑔^𝜇]=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 G and a predefined number k of clusters, METIS divides the graph in k clusters minimizing the number of edges with terminations in different clusters.

The results obtained are shown in Tables 1 to 3. These tables use the following legend in the columns(i)𝑘: number of clusters;(ii)𝑛: number of potential facilities locations and clients to be covered;(iii)𝑝: 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, 𝑈=150.

The values marked with an asterisk in Table 3 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 for 𝑘=10 (SJC instances) and 𝑘=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 Tables 1, 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, 𝑈=200.

Table 3 

Computation times for TSPLIB PCB3038 instance, 𝑈=400.

Comparing the values of CPUSeq and CPUPar for each value of 𝑘, 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 fixed a 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.

Acknowledgments

The authors acknowledge CNPq—Conselho Nacional de Desenvolvimento Científico e Tecnológico for partial financial research support. They also would like to thank the valuable suggestions of two anonymous referees.