Vertex Weighting-Based Tabu Search for p-Center Problem

1 Introduction

As a classical combinatorial optimization problem introduced by Hakimi [1964], the p-center problem consists of choosing p centers from a set of candidate centers to serve a set of clients, where each client is served by one of its closest centers. There is a serving arc if a center serves a client, and the length of the longest serving arc is the covering radius. The objective of this problem is to minimize the covering radius.

The p-center problem has a wide range of real-world applications. For instance, one can formulate several important problems in city planning as the p-center problem, such as determining the locations of emergency centers [Toregas et al., 1971] and fire stations [Tansel et al., 1983; Badri et al., 1998]. Regarding the supply-chain management, production plant location and warehouse distribution can be modeled as the p-center problem [Amiri, 2006]. In telecommunication, the p-center problem can formulate the reliable hub network design problem [Tran et al., 2016].

∗Corresponding author

As a challenging NP-hard problem [Kariv and Hakimi, 1979], the p-center problem has attracted much attention from academic society during last decades. The solution methods of the p-center problem can be mainly categorized into exact algorithms and metaheuristic algorithms. Regarding the exact algorithms, Minieka [1970] converted the p-center problem into a series of set covering problems. Based on Minieka [1970], Daskin [2002] solved the p-center problem using Lagrangian relaxation. Ilhan et al. [2002] proposed an algorithm which finds the lower bound by optimizing the linear relaxation of the original problem, then it improves this lower bound by proving the infeasibility of a series of decision problems via integer programming. Elloumi et al. [2004] independently designed a similar approach as Ilhan et al. [2002] and obtained better results. Calik and Tansel [2013] proposed an integer programming formulation and a decompositionbased exact algorithm for tackling the p-center problem.

Apart from the exact methods, various metaheuristic algorithms have been proposed for solving the p-center problem. Mladenovi´c et al. [2003] presented a tabu search and a variable neighborhood search. Caruso et al. [2003] proposed a hybrid algorithm which efficiently solves some smallscale p-center instances. Pullan [2008] presented a memetic algorithm which combines a local search procedure with a population-based metaheuristic algorithm. Davidovi´c et al. [2011] solved the p-center problem by bee colony optimization. Ferone et al. [2017] proposed a GRASP-based algorithm and Yin et al. [2017] combined path-relinking with GRASP to tackle the p-center problem.

This paper presents a vertex weighting-based tabu search (VWTS) algorithm which incorporates a vertex weighting technique and a tabu search procedure for solving the p-center problem. Distinguished from the previous metaheuristics which directly tackle the original problem, VWTS converts the p-center problem into a series of set covering problems. Specifically, for each possible covering radius r, it checks if each client can be served by any center within the covering radius r. Our contributions can be summarized as follows:

1. We transform the p-center problem to a series of decision subproblems, each of which is further converted to a new optimization problem to reduce the complexity.

1. By combining a vertex weighting technique with a tabu search procedure, and using incremental neighborhood evaluation, VWTS reaches better balance between exploitation and exploration despite its simplicity.
2. Tested on two sets of totally 138 most commonly used instances of the p-center problem, the proposed VWTS algorithm improves the best known results on 14 instances and matches the records on the remaining ones. Moreover, our computational time is much shorter than that of the state-of-the-art algorithms in the literature.
3. We conduct extensive experiments to justify the merits of the key components of our algorithm, which are the vertex weighting technique, the tabu search strategy, and the incremental neighborhood evaluation, respectively.

2 ProblemDescription and Transformation

The p-center problem is defined on an undirected complete graph G = (V,C,E), where V ∪ C is the vertex set and E is the edge set. Each vertex vi ∈ V corresponds to one of n clients to be served by one of the p centers chosen from the set of m candidate centers C. For each edge connecting vertices i and j, dij gives its length. The solution vector can be defined as (x, y, r), where variable xj = 1 iff a candidate center j ∈ C is opened, variable yij = 1 iff client i ∈ V is served by candidate center j ∈ C, and variable r ∈ R+ is the upper bound of the covering radius. Then, the classical mixed-integer programming (MIP) model for the p-center problem can be formulated to the following model (PC).

In model (PC), objective (1) aims to minimize the covering radius. Constraints (2)-(5) requires that there is at most p opened centers, each client must be served by exactly one center, only the opened centers can serve the clients, and the covering radius r is not shorter than any serving arc from each client to its assigned center, respectively. It is obvious that the optimal covering radius must be the same as the length of a certain edge dij . Hence, let Γ = {r1, r2, ..., rk} be an ordered list of distinct edge lengths where r1 < r2 < ... < rk, the pcenter problem can be considered as seeking for the smallest rank q such that model (PC) is still feasible when r ≥ rq, while it becomes infeasible if constraint r ≤ rq-1 is added.

When the optimal edge length rank q of an instance is given, the p-center problem is equivalent to the set covering problem [Chvatal, 1979]. Specifically, we use a vertex set V q j = {i ∈ V |dij ≤ rq} to denote the set of clients that candidate center j ∈ C can serve within covering radius rq. Thus, we can obtain a set covering instance V q = {V q 1 , V q 2 , ..., V qm}. If there exists p sets out of V q whose union contains all vertices, we can claim that the corresponding p centers are able to serve all clients within the covering radius rq. Therefore, the set covering problem defined by Eqs. (7)-(10) in the following model (SCq) is equivalent to model (PC) when the optimal edge length rank q is specified, where ui is a binary variable which is 1 if client i is not covered by any center, and xj is the same decision variable in model (PC).

Model (SCq) is slightly different from the standard set covering model or the MIP model presented by Elloumi et al. [2004]. Instead of minimizing the number of chosen sets, it minimizes the number of uncovered clients by exactly opening p centers as shown in Eq. (7). Constraints (8) ensure that for each client there is at least one center within the covering radius r. Constraint (9) limits the number of opened centers.

Unfortunately, we are usually lack of the priori knowledge about the optimal covering radius, so we need to iterate over the distinct edge length list Γ to examine each possible radius. In fact, by utilizing the bounds of the original problem, the number of subproblems can be greatly reduced.

3 VertexWeighting-Based Tabu Search

In order to tackle the p-center problem, the proposed VWTS algorithm combines the tabu search strategy and the vertex weighting technique to solve a series of subproblems. Starting from an upper bound rq0 obtained by executing solvers for model (PC) such as the PBS algorithm [Pullan, 2008] under limited time, VWTS solves models (SCq0-1), (SCq0-2), ..., (SC1) in turn until it fails to find any feasible solution for model (SCq∗-1) within the given time limit. Then, rq∗ is the best covering radius found. Therefore, our focus will be on the solution method for model (SCq) with a given radius rq.

The main framework of the proposed VWTS algorithm is presented in Algorithm 1. It generates an initial solution X by a greedy algorithm (line 1), and iteratively improves the incumbent solution by a tabu search procedure (lines 4-14). At each iteration of the VWTS algorithm, it first evaluates the neighborhood of the current solution and records the best neighborhood move while respecting their tabu states (line 5). Then, it makes the best move and replaces the current solution with the resulting neighboring solution (line 6). If the current solution X improves the best solution found so far, then X∗ is updated with X (lines 7-8). Otherwise, we need to check whether the current solution falls into a local optimum, that is, the best move returned by function FindPair() cannot reduce the number of uncovered clients (line 9). If stagnation occurs, the weight of each uncovered client is adjusted (line

Algorithm 1 The main framework of the VWTS algorithm

Input: A graph G, a center number p, a covering radius rq Output: The best solution found so far X∗

• 1: A set of p centers X ← Init(G, p, rq) /* (Section 3.1) */
• 2: X∗ ←X, X′ ←X, tabu list TL←∅, iter ← 1
• 3: Vertex weights wi ← 1, ∀i ∈ V
• /* (Section 3.2) */
• 4: while termination condition is not met do
• 5: (i, j) ← FindPair(X′ , TL, iter) /* (Algorithm 2) */
• 6: MakeMove(i, j) /* (Algorithm 4) */
• 7: if |U(X)| < |U(X∗ )| then /* U(X) is the set of */
• 8: X∗ ← X /* clients uncovered by X */ 9: else if |U(X)| ≥ |U(X′ )| then
• 10: wv ← wv + 1, ∀v ∈ U(X) /* (Section 3.2) */
• 11: end if /* more uncovered clients than last solution */
• 12: TL← {i, j} /* update tabu list (Section 3.4) */
• 13: X′ ←X, iter ← iter +1
• 14: end while

10). Finally, when the specified termination condition is met, the algorithm terminates and returns the best solution X∗ .

3.1 Initialization

The VWTS algorithm employs a constructive heuristic for generating a set of centers X as the initial solution. Let Vj denote the set of clients that candidate center j can serve within the current covering radius. Ci gives the set of candidate centers which are able to serve client i, that is, Ci = {j ∈ C|i ∈ Vj}. U(X) = V \ ⋃j∈X Vj represents the set of clients that are not served in solution X. The constructive heuristic opens centers one by one under a maximal coverage principle. Precisely, it iteratively selects a candidate center j = argmaxj∈C\X |Vj ∩ U(X)| which covers most uncovered clients, and inserts it into the current solutionX. If there are multiple candidate centers covering the same number of uncovered clients, ties are broken randomly. The time complexity of the construction is O(np).

3.2 VertexWeighting Technique

The vertex weighting technique helps the search to escape from the local optima by altering the objective function. It can be regarded as a variant of guided local search [Voudouris and Tsang, 2003] and has been successfully applied to many problems, such as unicost set covering problem [Gao et al., 2015], minimum vertex cover problem [Cai et al., 2011] and satisfiability problem [Luo et al., 2012]. In the p-center problem, we adapt objective (7) as Eq. (11) presents.

Therefore, the VWTS algorithm actually works on model (SC w q ). Note that it is a dynamic model where wi varies as the search proceeds. If the VWTS algorithm keeps failing to cover a client, it implies that this vertex is hard to cover and we should treat it with higher priority. Specifically, when the tabu search is trapped in local optimal solution X, the VWTS algorithm increases the weight wi of each uncovered client i ∈ U(X) by one unit (Algorithm 1, lines 9-11). Nevertheless, the objective value of an optimal solution is always zero regardless of the configuration of the weights. This process changes the landscape of the solution space such that X is no longer a local optimum, and the search will be able to continue to explore other search areas. The more frequent a client appears in U(X) when encountering stagnations, the greater its weight will be. On the one hand, the vertex weighting technique is able to prevent vertices from being repeatedly uncovered and diversify the search in an adaptive manner. On the other hand, it modifies the solution space in a smooth way, which can guide the search to promising search regions.

3.3 Neighborhood Structure and Evaluation

In order to improve the initial solution under model (SC w q ), VWTS employs a swap-based neighborhood structure inspired by most local search-based metaheuristics for the pcenter problem [Pullan, 2008], while the neighborhood evaluation is significantly different due to the reformulation presented in Section 2 and the following acceleration strategies. Denoted by Swap(i, j), a swap move produces a neighboring solution X ⊕ Swap(i, j) = X ∪ {i} \ {j}, by opening a candidate i ∈ C \ X, and closing an opened center j ∈ X.

The neighborhood evaluation is the most time consuming procedure in a trajectory-based metaheuristic algorithm. For a typical tabu search algorithm which adopts the best improvement policy, it evaluates all feasible moves at each iteration, and performs one of the best neighborhood moves which improves the objective value as much as possible. Since there are O(p(n - p)) swap moves, the size of the neighborhood could be huge on some large instances. Therefore, we use a neighborhood sampling strategy and an incremental evaluation technique to accelerate the evaluation.

On the one hand, the objective value can only be improved by covering some uncovered clients, so the VWTS algorithm will only evaluate a swap move Swap(i, j) if i covers some uncovered vertices in U(X). Since each client must be eventually covered, the VWTS algorithm randomly picks a vertex k ∈ U(X), and only evaluates neighborhood moves Swap(i, j) where i ∈ Ck and j ∈ X. Apart from the reduced time consumption for the neighborhood evaluation, the diversification of the search is improved as a side effect, reaching a better balance between intensification and diversification.

On the other hand, we try to accelerate the neighborhood evaluation by reusing some intermediate results for calculating the objective values. Instead of naively summing the weights of the covered clients according to objective (11), VWTS incrementally evaluates all neighborhood moves by storing and maintaining δj, the objective value increase (decrease) by closing (opening) center j. For each center j ∈ X, δj = ∑i∈(Vj∩U(X\{j})) wi is the sum of the weights of the vertices which can only be served by center j. For each candidate center j /∈ X, δj = ∑i∈(Vj∩U(X)) wi is the sum of the weights of all uncovered clients in Vj. Then, the evaluation for each neighborhood move can be implemented in O(1) time complexity, at the cost of updating the affected δ values every time after opening or closing a center. Specifically, it calculates the objective value for opening center i by f(X ∪ {i}) = f(X) - δi. Then, we need to update the

Algorithm 2 Find the best swap pair

1: function FINDPAIR(X, TL, iter) 2: The set of best swap moves M ←∅ 3: The best objective value obj ←+∞ 4: v ← a randomly picked uncovered vertex in U(X) 5: δ′j ← δj, ∀j ∈ C /* backup before trial moves */ 6: for all i ∈ Cv do /* Cv: candidates covering v */ 7: TryToOpenCenter(i) /* (Algorithm 3) */ 8: for all j ∈ X do /* evaluate closing center j */ 9: if {i, j} ∩ TL = ∅ then /* not tabu move */ 10: if f(X ⊕ Swap(i, j)) < obj then 11: obj ← f (X ⊕ Swap(i, j)) 12: M ← {Swap(i, j)} 13: else if f(X ⊕ Swap(i, j)) = obj then 14: M ← M ∪ {Swap(i, j)} 15: end if 16: end if 17: end for 18: δj ← δ′j , ∀j ∈ C /* restore after trial moves */ 19: end for /* v ∈ U(X) ⇔ Cv ∩ X = ∅ */ 20: return a randomly picked move in M 21: end function

Algorithm 3 Open a center virtually

1: function TRYTOOPENCENTER(i) 2: for all v ∈ Vi do /* |X ∩ Cv|: number of centers */ 3: if |X ∩ Cv| = 1 then /* covering v in X */ 4: /* cancel penalty for making v uncovered */ 5: δl ← δl - wv, for l ∈ X ∩ Cv /* O(1) */ 6: end if /* l was the only center covering v but */ 7: end for /* it will not be the only one if i opens so */ 8: end function /* closing l does not make v uncovered */

δ values affected by opening center i which eliminates the penalties for closing some centers due to the overlapped coverage. After that, we can evaluate the neighboring solutions by f(X ⊕ Swap(i, j)) = f(X ∪ {i}) + δj . As tabu search procedure evaluates a large number of moves but only makes a single move at each iteration, it is worthy of maintaining and querying the cache rather than computing the objective function from scratch for every swap move.

Algorithm 2 describes the neighborhood evaluation procedure. It randomly selects a client k ∈ U(X) (line 4) and tries to open each candidate center which covers vertex k (lines 67). The sub-routine TryToOpenCenter(i) keeps each δj upto-date to accelerate the calculation of the objective function f(X ⊕ Swap(i, j)). For each non-tabu neighboring solution, if its objective value is smaller than obj, the corresponding neighborhood move is saved as the best move (lines 9-16). When all trial moves regarding opening candidate center i are evaluated, we need to restore the δ values (line 18).

Apparently, covering a client v which has already been covered (|X ∩ Cv| ≥ 1) will not improve the objective value, and closing the center which serves a client v will not worsen the objective value if there are multiple centers which can cover vertex v (|X ∩ Cv| ≥ 2). Therefore, as presented in Al-

Algorithm 4 Make a swap move

1: function MAKEMOVE(i, j) 2: for all v ∈ Vi do /* consequences of opening i */ 3: if |X ∩ Cv| = 1 then /* (Algorithm 3) */ 4: δl ← δl - wv, for l ∈ X ∩ Cv 5: else if |X ∩ Cv| = 0 then 6: δl ← δl - wv, ∀l ∈ Cv \ {i} 7: end if /* cancel reward for covering v */ 8: end for 9: X ←X ∪ {i} \ {j} 10: for all v ∈ Vj do /* consequences of closing j */ 11: if |X ∩ Cv| = 0 then /* add reward for */ 12: δl ← δl +wv, ∀l ∈ Cv \ {j} /* covering v */ 13: else if |X ∩ Cv| = 1 then 14: δl ← δl + wv, for l ∈ X ∩ Cv 15: end if /* add penalty for uncovering v */ 16: end for 17: end function

gorithm 3, if there is exactly one center l which covers client v (|X ∩ Cv| = 1) before opening center i, then the δ value for closing center l which used to make client v uncovered will decrease by wv (lines 4-7), because closing center l will not make client v uncovered once center i is opened. We are only interested in δj(∀j ∈ X) for closing a center, since lines 8-17 in Algorithm 2 only concern closing another center. As a result, the worst-case time complexity of Algorithm 2 and Algorithm 3 are O(n2 ) and O(n), respectively.

When the best swap move is performed (Algorithm 1, line 6), we need to update the affected data structures by calling Algorithm 4. In addition to updating the δ values as Algorithm 3 (lines 3-4), if there is no center covering client v (|X ∩ Cv| = 0) before opening center i, then the δ value for opening each candidate center l to cover client v is decreased by wv (lines 5-7), because opening center l to cover the already covered client v in the future will no longer improve the objective value by wv. Then, it updates the set of opened centers (line 9). The influences on δ values for closing center j are taken into account in a similar manner (lines 10-16).

3.4 Tabu Search

Tabu search usually incorporates a recency-based tabu list to prohibit revisiting recently visited solutions. The tabu strategy prevents closing newly opened centers or reopening newly closed ones immediately. We fix the parameter of the tabu tenure in the tabu strategy to one iteration, so the proposed algorithm keeps simple and parameter-free. Specifically, if we open (close) a center at the current iteration iter, it is forbidden to close (open) it again at the the next iteration. Thus, Swap(i, j) at iteration iter introduces two vertices {i, j} in tabu list TL (Algorithm 1, line 12) and neither i nor j can be involved at iteration iter + 1 (Algorithm 2, line 9).

4 Experimental Results and Comparisons

In order to evaluate the effectiveness of the proposed VWTS algorithm, we conduct extensive experiments on the wellknown datasets, and compare VWTS with the state-of-theart algorithms including two exact algorithms (ELP [Elloumi et al., 2004] and DBR2 [Calik and Tansel, 2013]) and three metaheuristic algorithms (PBS [Pullan, 2008], GRASP+PS [Ferone et al., 2017], and GRASP/PR [Yin et al., 2017]).

4.1 Experimental Protocol

Our proposed VWTS algorithm is programmed in C++ and compiled with Visual Studio 2017. All experiments are carried out on Windows Server 2012 x64 with Intel Xeon E52609v2 2.50 GHz CPU. We run a reimplemented PBS algorithm [Pullan, 2008] under a 5-millisecond time limit to obtain a good initial upper bound of the covering radius rq0 for each instance. We performed 20 independent runs on each instance, and the cutoff time for each subproblem induced by each covering radius is 6 minutes. Since the programs of the reference algorithms are not available, we report their computational time in the corresponding papers and normalize the run time according to CPU frequency for a fair comparison. The run time for ELP is divided by 6.25 for it runs on a 400 MHz Pentium II CPU. The remaining algorithms are tested on faster (PBS, GRASP+PS, GRASP/PR) or unknown (DBR2) CPUs, so we report their computational time as is.

There are mainly two sets of benchmark instances for the p-center problem. The first set consists of 40 small instances (pmed) from OR-Library [Beasley, 1990]. The second set contains 98 instances from TSPLIB which are based on planar graphs and derived from real-world applications [Reinelt, 1991]. The TSPLIB dataset can be further divided into 44 small instances (sTSP) with less than 1000 vertices, and 54 large ones (u1060, rl1323, u1817, pcb3038). A client is also a candidate center (V = C) in these instances. We compute the all-pair shortest paths by Floyd algorithm [Floyd, 1962] for the pmed instances, and round the Euclidean distances to the nearest hundredths for the TSPLIB instances, as Pullan [2008] did. Note that algorithm ELP and DBR2 adopt integer distances which are too rough to figure out the optimal q∗ .

4.2 Computational Results

In this section, we provide comprehensive investigations on the performance of our VWTS algorithm and report the detailed results in Tables 1 and 2. They follow the same convention described below. Columns Dataset and Instance give the names of datasets and instances, respectively. Columns n and p indicate the numbers of vertices and centers, respectively. Column Count presents the number of instances in the dataset. Column CPU gives the normalized average CPU time in seconds. For VWTS, it includes the total time for computing rq0 and sequentially solving all subproblems from (SCq0) to (SCq∗). Column CPU-q∗ reports the run time in seconds for VWTS to solve the subproblem induced by the best known covering radius rq∗, which shows the limit of VWTS when parallel computing is available. Column fbest presents the best objective values obtained by the corresponding algorithms. Column #best shows the number of instances where the corresponding algorithms match the best known results. The #better, #equal, and #worse are the number of instances where VWTS obtains better, equal, and worse results comparing to the corresponding algorithms, respectively.

Table 1 compares the overall results on each dataset obtained by ELP, DBR2, PBS, GRASP/PR, and our VWTS. We can observe that VWTS obtains the best known results for all the instances, while no reference algorithm can obtain this overall outcome. Moreover, the average computational time of VWTS on each dataset is over 4 times shorter than the best algorithms in the literature. If parallel computing is available, the performance gain can be over 30 times. The time difference is so big that the comparison is not biased by the normalization. In addition, VWTS is very stable since it always obtains the reported results in all 20 independent runs.

Table 2 compares the experimental results produced by ELP, DBR2, PBS, GRASP+PS, GRASP/PR, and our VWTS on 30 largest and most challenging instances. Note that the real upper bounds for dataset pcb3038 have not been updated for decades (although DBR2 improved some integer bounds). From Table 2 we can observe that, the proposed VWTS algorithm improves the previous best known results obtained by GRASP/PR and PBS on rl1323 with p = 100, u1817 with p = 80, 130, 150, and pcb3038 with p ≥ 50. Apart from the improved solution quality, VWTS outperforms all reference algorithms on u1817 and pcb3038 datasets in terms of the computational time. Regarding the run time for the subproblem corresponding to the optimal or best known covering radius, VWTS is able to obtain the new best known results on 23 out of 30 instances within 3 seconds. Furthermore, if VWTS stops at the previous best covering radius, it only takes 3 seconds on 28 instances. In sum, these statistics show that the our VWTS is highly effective, efficient, robust, and concurrency-friendly for solving the p-center problem.

4.3 Importance of VWTS Ingredients

In order to evaluate the merits of the vertex weighting technique, the tabu search strategy, and the incremental neighborhood evaluation, we conducted experiments on four largest instances (pcb3038, p = 350, 400, 450, 500) to compare VWTS with its simplified versions obtained by disabling the tabu strategy (VW), deactivating the vertex weighting technique (TS), and using the naive neighborhood evaluation instead of the incremental one (VWTS′ ), respectively.

Figure 1 depicts the evolution of the uncovered vertex numbers as the search proceeds by VWTS, VW, TS, and VWTS′ . Each point (x, y) means that there are 10y - 1 clients which are not served by any center at 10x microsecond. We can observe that, although TS can obtain better solutions than VWTS at first milliseconds, VWTS and VW overtake after 0.1 seconds. The reason for this phenomena might be that, the tabu search pays too much attention on prohibiting a candidate center from being centers, instead of directly preventing a client from being uncovered. In extreme conditions, some hard-to-serve vertices will never get a chance to be served at all. However, the tabu strategy accelerates the convergence of VWTS. According to Figure 1, VWTS is over 100.2 = 1.5 times faster than VW for finding a feasible covering. Moreover, when the incremental evaluation is disabled, VWTS′ is 100 times slower than the complete algorithm. These observations confirm that the vertex weighting technique, the tabu strategy, and the incremental evaluation are all essential for VWTS in terms of effectiveness and efficiency.

5 Conclusion

This paper presented a vertex weighting-based tabu search algorithm to solve the p-center problem. We decompose the pcenter problem into a series of set covering subproblems, and solve them by combining the vertex weighting technique and the tabu search strategy. Tested on 138 widely used datasets, VWTS improves the best known results on 14 instances and matches the records in the literature for all remaining ones within less run time. Thus, it would be appealing to investigate the combination of the proposed strategies for solving other optimization problems in the future.

Acknowledgments

Research partially supported by Alkaid lab of Huawei Cloud.