A Weighting-Based Tabu Search Algorithm for the p-Next Center Problem

1 Introduction

The p-center problem is a classical discrete optimization problem and has been proven to be NP-hard [Kariv and Hakimi, 1979]. It consists of selecting p centers from a set of candidate centers and serving a set of clients in order to minimize the maximum distance between each client and its closest center. The p-center problem has a wide range of applications, such as the problems of determining the locations of emergency centers [Toregas et al., 1971], hospitals [Hakimi, 1964], and fire stations [Drezner, 1987] to serve the communities. In the last 50 years, the p-center problem has been widely concerned in academic society. There are mainly three kinds of methodologies, including exact algorithms (Minieka

∗Co-corresponding authors.

[1970], Elloumi et al. [2004], Calik and Tansel [2013], Liu et al. [2020]), approximation algorithms (Hochbaum and Shmoys [1985], Martinich [1988]), and metaheuristic algorithms (Mladenovi´c et al. [2003], Pullan [2008], Irawan et al. [2016], Yin et al. [2017], Zhang et al. [2020]), respectively.

In real-world applications, a center may be forced to be closed due to an unpredictable incident. In that case, the clients served by the center have to turn to another center from the current center as quickly as possible. Therefore, Albareda-Sambola et al. [2015] introduced the p-next center problem (pNCP), which occurs in many industrial applications and has been proven to be NP-hard. It requires locating p centers from a set of candidate centers and assigning a reference center and a backup center to each client. The objective of pNCP is to minimize the maximum cost of all clients, which is calculated by the distance of a client to its closest center (reference center), plus the distance from the reference center to its closest center (backup center). Figure 1 shows an illustrative diagram of pNCP.

Albareda-Sambola et al. [2015] also proposed a series of different formulations for pNCP and compared their performances. L´opez-S´anchez et al. [2019] proposed a hybrid metaheuristic algorithm, which incorporates Greedy Randomized Adaptive Search Procedure (GRASP) and Variable Neighborhood Search algorithm (VNS) to solve the pNCP. Londe et al. [2021] introduced a Biased Random-Key Genetic Algorithm (BRKGA), which combines a multi-parent (MP) strategy and a path-relinking (PR) procedure. They conducted experiments on 413 instances derived from p-median [Mladenovi´c et al., 2007] in the OR-Library [Beasley, 1990].

In this paper, we present a weighting-based tabu search (WTS) algorithm for solving the pNCP. WTS uses the weighting technique and a solution-based tabu search framework to strengthen the search capability. Instead of directly optimizing the original problem like previous metaheuristics, WTS converts the original pNCP into a series of decision subproblems and solves them with a fast tabu search procedure which adopts a neighborhood structure based on the closingopening swap and an incremental neighborhood evaluation technique based on a hierarchical penalty strategy.

Our main contributions can be summarized as follows:

1. We propose a new solution approach which transforms the pNCP from an optimization problem to a series of decision subproblems and solves each of them in turn to obtain the solution of the original optimization problem. It improves the smoothness of the solution space, which makes it easier for the local search-based algorithms to solve the pNCP.
2. We present a weighting technique for solving the pNCP, which diversifies the search trajectory and helps the search to escape from the local optima trap.
3. WTS algorithm employs a new neighborhood structure based on a dedicated swap move which is made up of a closing operation and an opening operation. In addition, we use a new hierarchical penalty strategy for neighborhood evaluation and an incremental evaluation strategy to improve the efficiency of search.
4. Tested on 413 classical instances of the pNCP, the proposed WTS algorithm improves the best known results on 12 instances and obtains the optimal solutions on 10 instances for the first time, and matches the records on the remaining ones. In addition, the computational time of WTS is much shorter than that of the state-of-the-art algorithms in the literature.
5. We conduct additional experiments to analyze the importance of the key components of our algorithm, such as the solution-based tabu strategy, the weighting technique, and the hierarchical penalty strategy.

2 ProblemDescription and Transformation

Given an undirected complete graph G = (V,E), where V = {1, 2, ..., n} is the vertex set and E = {(i, j)|i, j ∈ V } is the edge set. The distance between each pair of vertices i and j is dij. Each vertex i ∈ V corresponds to a client to be served by a reference and a backup center, where the two centers are both one of the p centers, which should be selected from a set of candidate centers C ⊆ V (usually C = V ). The goal of pNCP is to choose p best centers from C such that the maximum length of all serving paths can be minimized. We define dijk as the total distance from a vertex i to a center j and then to another center k, where dijk = dij+djk, i ∈ V, j, k ∈ C. The solution can be defined as (x, y, z). In detail, x = {xj|j ∈ C} and xj is a decision variable which equals to 1 if and only if a candidate center j ∈ C is opened as a center. y = {yijk|i ∈ V, j, k ∈ C} and yijk is a binary variable, where yijk = 1 iff the reference center and backup center of client vertex i ∈ V are j, k ∈ C, respectively. z ∈ R+ is the upper bound of the assignment cost, i.e., the length of the longest path. Based on the above notations, the classical mixed-integer programming (MIP) model [Albareda-Sambola et al., 2015] for the pNCP can be formulated as follows.

(5)

Objective (1) aims to minimize the assignment cost. Constraint (2) limits the number of opened centers to at most p. Constraints (3) ensure that each vertex must be assigned to one reference and one backup center. Constraints (4) and (5) restrict that only the opened centers can be assigned to the vertices. Constraints (6) force that the reference center of each vertex must be its closest center. Constraints (7) ensure that the cost z is not less than the length of any path from each vertex to its reference center and then to its backup center.

We can transform pNCP to a series of decision subproblems called z-cost pNCP, where the maximum length of serving paths is fixed as z, and we need to locate p centers under the constraint that the actual cost is no more than z. Once the z-cost pNCP is solved, we can find a feasible solution for pNCP with cost equal to z. Obviously, the value of z must be one of dijk. Let Z = {z1, z2, . . . , zm} be the set containing all distinct values of dijk, where z1 < z2 < zm. Then, we can optimize the pNCP by decreasing z from zm to z1 and solving each corresponding z-cost pNCP.

3 VertexWeighting-Based Tabu Search

3.1 General Framework

The main framework of the WTS algorithm is described in Algorithm 1. First, we employ a constructive heuristic to obtain an initial feasible solution and the upper bound of the cost. Then, WTS converts pNCP into a series of z-cost pNCP decision subproblems and solves them by a solution-based tabu search procedure combined with a weighting technique.

Algorithm 1 The main framework of the WTS algorithm

Input: A graph G, the center number p

Output: The best solution found so far X∗ 1: Solution tabu list ST ←∅ 2: iter ← 1, vertex weights wi ← 1,∀i ∈ V 3: Age aj ← 0, ∀j ∈ C 4: Current solution X, initial cost zq ← Initialize(G, p) 5: Previous solution X′ ←X, X∗ ←X 6: The ordered list of z-costs Z = {z1, z2, ..., zq} 7: while Termination condition is not met do 8: q ← q - 1, update set of illegal vertices I(X) 9: while Termination condition is not met do 10: (i, j) ← FindPair(X′ , ST, iter) // Algorithm 4 11: X ← X ∪ {i} \ {j} 12: if |I(X)| = 0 then 13: X∗ ← X 14: break 15: else if |I(X)| > |I(X′ )| and |I(X)| ≤ 5 then 16: wk ← wk + 1, ∀k ∈ I(X) // Section 3.2 17: end if 18: ST ← ST ∪ {X}, X′ ← X, iter ← iter + 1 19: end while 20: end while 21: return X∗

Algorithm 2 Closing a center

1: function CLOSINGCENTER(j) 2: for all l ∈ X do 3: if FY 1 l = j then 4: FY 1 l ←FY 2 l ,DY 1 l ←DY 2 l 5: (FY 2 l ,DY 2 l ) ←FindNextFY(j, l) 6: else if FY 2 l = j then 7: (FY 2 l ,DY 2 l ) ←FindNextFY(j, l) 8: end if 9: end for 10: for all v ∈ Sj do 11: if FN 1 v = j then 12: FN 1 v ←FN 2 v , DN 1 v ←DN 2 v 13: (F N 2 v , DN 2 v ) ← FindNextFN(j, v) 14: else if FN 2 v = j then 15: (F N 2 v , DN 2 v ) ← FindNextFN(j, v) 16: end if 17: t ← FN 1 v 18: Dv ← DN 1 v +DY 1 t 19: update I(X) and δj 20: end for 21: end function

As soon as the zq-cost pNCP is solved, we decrease zq by setting zq to zq-1 and solve the resulting subproblem again.

Specifically, we first generate an initial solution X with an initial cost zq by a greedy construction algorithm (line 4). The constructive heuristic first picks a center randomly, and then it opens centers one by one under the principle of minimizing objective (1) until |X| = p. According to the initial cost zq, we initialize the ordered cost list Z with all distinct path lengths (line 6), where z1 < z2 < ... < zq. Then, we solve the zq-cost subproblems in turn by decreasing q to q - 1. When solving a subproblem, we relax constraints (7) where the assignment cost of each vertex is no more than zq and record every vertex that violates the constraint as an illegal vertex. For each subproblem, we first record I(X) by the illegal vertices (line 8), and then optimize it iteratively by tabu search until |I(X)| = 0 or the time limit is met (lines 9-20). At each iteration of WTS, it evaluates the neighborhood of the current solution X and selects the best neighborhood move in non-tabu status (line 10, Algorithm 4). If the number of illegal vertices in the current solution X is 0, i.e., |I(X)| = 0, it means that a legal solution for the current decision subproblem is found, and the best solution X* will be updated with X (lines 12-14). Otherwise, when the current solution is trapped into a local optimum, i.e., the best move returned by FindPair() cannot reduce the number of illegal vertices, WTS will increase the weight of each illegal vertex to encourage eliminating their constraint violations in the next iterations (lines 15-16). At the end of each iteration, the associated data structures will be updated (line 18). Finally, once the specified termination condition (e.g., time limit or maximum number of iterations) is met, WTS terminates and returns the best solution X* (line 21).

3.2 Reformulation andWeighting Technique

The weighting technique helps the search to escape from the local optima by introducing an adaptive objective function. It has been successfully applied to many problems, such as unicost set covering problem [Gao et al., 2015], classical p-center problem [Zhang et al., 2020], and optimal camera placement problem [Su et al., 2021]. We transform pNCP into a series of decision subproblems, which determine whether all vertices can assign a reference center and a backup center within specific z-cost. In order to effectively solve each subproblem, we relax the constraints that each vertex must be allocated within z-cost and impose a penalty for each illegal vertex on the objective function. Specifically, let wi be the weight of vertex i ∈ I(X) ⊂ V , the objective function f(X) can be defined as Eq. (10).

The initial weight of each vertex is 1. In the subsequent search, if WTS fails to reduce the number of illegal vertices, i.e., the tabu search is trapped in the local optima and the objective value is below the threshold, WTS will increase the weight wi of each illegal vertex i ∈ I(X) by one unit (Algorithm 1, lines 15-16). Specifically, at the early stage of the neighborhood search, tabu strategy is sufficient for jumping out of the local optima trap, so the threshold is set to 5 in WTS. In the later stage of the search, once stagnation is encountered, the weight of the corresponding illegal vertices will be increased, and the more times a vertex appears in I(X), the greater its weight will be. The weighting strategy changes the landscape of the solution space, and is thus able to effectively prevent vertices from being illegal for a long time.

Algorithm 3 Opening a center

1: function OPEINGCENTER(i) 2: for all l ∈ X do 3: if DY 1 l > dli then 4: FY 2 l ←FY 1 l ,DY 2 l ←DY 1 l 5: FY 1 l ← i,DY 1 l ← dli 6: else if DY 2 l > dli then 7: FY 2 l ← i,DY 2 l ← dli 8: end if 9: end for 10: for all v ∈ Si do 11: if DN 1 v > dli then 12: FN 2 v ←FN 1 v , DN 2 v ←DN 1 v 13: FN 1 v ← i,DN 1 v ← dvi 14: else if DN 2 v > dli then 15: FN 2 v ← i,DN 2 v ← dvi 16: end if 17: t ← FN 1 v 18: Dv ← DN 1 v +DY 1 t 19: update I(X) and δi 20: end for 21: end function

Algorithm 4 Find the best swap pair

1: function FINDPAIR(X, ST, iter) 2: Best swap move (i∗ , j ∗ )←∅, δ ←+∞ 3: for all j ∈ X do // evaluate closing center j 4: ClosingCenter(j) // for evaluation only 5: if (δj < δ) ∨ ((δj = δ) ∧ (aj > a∗j )) then 6: δ ← δj , j ∗ ← j 7: end if 8: end for 9: δ ← +∞ 10: for all i ∈ C \ X do // evaluate opening center i 11: if X ⊕m(i, j ∗ ) /∈ ST then 12: OpeningCenter(i) // for evaluation only 13: if δi < δ then 14: δ ← δi, i∗ ← i 15: end if 16: end if 17: end for 18: return (i∗ , j ∗ ) 19: end function

3.3 Neighborhood Structure and Evaluation

In order to solve each z-cost pNCP subproblem, WTS adopts a swap-based neighborhood structure. Specifically, a swap move consists of closing a center j ∈ X and opening a center i ∈ C \ X, denoted as m(i, j). Based on the incumbent solution X, performing a move produces a new neighborhood solution X ⊕m(i, j). At each neighborhood evaluation, we perform the best closing operation out of O(p) ones and the best opening operation out of O(n - p) ones in turn. So the size of the complete neighborhood is O(p+(n-p)) = O(n).

In local search procedure, the neighborhood evaluation is an extremely important ingredient. To obtain the best neighboring solution and improve the incumbent solution X, the proposed WTS algorithm uses the best-improvement policy, i.e., for a swap move based on closing center and opening center operations, it evaluates all neighborhood moves and performs the closing and opening operations that lead to the best neighboring solution. To evaluate a neighborhood move of the current solution, we need to allocate a reference and a backup center for each client to obtain its assignment cost for updating I(X), which means that the time complexity of the naive neighborhood evaluation is O(2np). However, it is quite time-consuming for large instances. In order to improve efficiency without dramatically sacrificing the solution quality, WTS follows a hierarchical penalty strategy for neighborhood evaluation and an incremental evaluation technique.

On the one hand, in order to accelerate the neighborhood evaluation procedure, WTS uses a set of dedicated data structures to store the k-nearest centers to each vertex and each center, and the corresponding distances, as shown in Table 1, where k is set to 1 and 2. With this set of data structures, the status of vertices can be quickly calculated each time closing or opening a center. Specifically, the assignment cost of each vertex i ∈ V is Di = DN 1 i + DY 1 j (j = FN 1 i ), and the assignment of vertex i is legal iff Di ≤ zq. In addition, WTS adopts an incremental evaluation mechanism to efficiently evaluate all neighborhood moves instead of naively calculating the objective by Eq. (10). Specifically, there are only small changes in the neighboring solutions, thereby calculating their objective by naively evaluating function f is unnecessary. The objective of a neighboring solution can be calculated incrementally by the variation of closing a center i and opening a center j, i.e., f(X ⊕m(i, j)) = f(X) + (δi + δj ), where δ represents the variation of the objective value after closing a center or opening a center. Algorithms 2 and 3 update the data structure and calculate the improvement δ (line 19) after closing center j and opening center i, respectively, where S will be introduced in Section 3.4. The worstcase complexities of ClosingCenter() and OpeningCenter() are O(np) and O(n), respectively.

On the other hand, WTS uses a hierarchical penalty strategy to reformulate the neighborhood evaluation function. We divide the illegality of vertices into two parts, and use different penalty-factors for them. Specifically, when the distance from the illegal vertex i ∈ I(X) to its reference center is not greater than zq, i.e., DN 1 i ≤ zq, the penalty factor is 1. Otherwise, the penalty factor is 2. Combining the hierarchical penalty strategy described in this section and Eq. (10) in Section 3.2, the new evaluation function becomes Eq. (11).

In the neighborhood evaluation process, when there exist multiple best neighboring solutions according to the reformulated evaluation function Eq. (11), we employ an age-based strategy to break ties. Specifically, the age of a center denotes the number of iterations for which it remains opening. The larger the age of a center is, the more iterations the center has been opened for. Thus, we break ties by favoring the neighboring solution which closes an older center. In addition, to avoid re-closing the just-opened center, the center which is just opened will not be evaluated at the next iteration.

Algorithm 4 shows the neighborhood evaluation procedure. In the worst-case, the time complexity is O(n2 ).

3.4 Neighborhood Reduction

Obviously, if the assignment cost of vertex i satisfies Di ≤ zq, then dij ≤ zq and djk ≤ zq, where j and k are the reference and backup of vertex i, respectively, i.e., if the total cost of a vertex i is not more than zq, the distance between i and its reference center, and the distance between the reference center and its backup center are not more than zq. So we adopt a neighborhood reduction rule to accelerate the neighborhood search procedure. Specifically, we ignore the centers which cannot reduce the cost of vertex i, and save the promising centers with the set Si = {j|dij ≤ zq, j ∈ C}.

3.5 Solution-Based Tabu Search Strategy

Tabu search is an optimization algorithm based on local search framework [Glover, 1989]. In recent years, the solution-based tabu search strategy has been widely used to solve many combinatorial problems, such as minimum differential dispersion problem [Wang et al., 2017], 0-1 multidimensional knapsack problem [Lai et al., 2018], and inventory routing problem [Su et al., 2020]. WTS algorithm adopts a solution-based tabu search strategy. Different from the classic attribute-based tabu strategy, the solution-based tabu strategy records complete information of each visited solution, and these solutions will never be evaluated again. This policy avoids the drawbacks of the attribute-based tabu strategy, which may falsely accept already visited solutions or forbid some promising moves.

At each iteration, making a move will produce a new current solution X to be stored into the tabu list ST, which is a Boolean vector of length L = 108 . Specifically, WTS defines a hash function h to map the solution X to an integer in the range [0, L), denoted as h(X), where ST (h(X)) = 1 means that the solution with the hash value h(X) is in tabu status. As there may be hash collisions that lead to incorrect identification of the tabu status of some unvisited solutions, we use multiple hash functions ht (t = 1, 2, 3) and corresponding ST t (ht (X)) to reduce the collision rate. Similarly, the solution X is in tabu state iff ∧3t=1 ST t (ht (X)) = 1. The hash function ht is defined as Eq. (12).

4 Experiments and Analysis

In order to evaluate the performance of the WTS algorithm, we conduct extensive experiments on 413 public benchmark instances, and compare WTS with the state-of-the-art algorithms in the literature, including GRASP-VNS [L´opezS´anchez et al., 2019], ILS [Londe et al., 2021], and BRKGA [Londe et al., 2021].

4.1 Experimental Protocol

There are 413 instances for the pNCP generated from 40 pmedian instances from OR-Library [Londe et al., 2021]. The number of vertices n ranges from 10 to 900, and the number of centers p distributes between 5 and 450. The shortest path between each pair of vertices is computed by Floyd algorithm [Floyd, 1962] for all instances.

All our experiments are run on Windows x64 with an Intel Core i7-10700 2.90 GHz CPU. We performed 20 independent runs on each instance under a 60-second time limit. The computational platform for the reference algorithm GRASPVNS is Intel Core 2 Duo E8400 (3 GHz) with 4 GB RAM. The reference algorithms ILS and BRKGA are tested on an Intel Xeon E5530 CPU at 2.40 GHz and 120 GB of RAM running CentOS Linux and the time limit for each instance is 1800 seconds. The MIP models proposed by [AlbaredaSambola et al., 2015] were solved with IBM ILOG CPLEX 12.10 solver and the cutoff time for each instance was 7 days.

4.2 Computational Results

In order to test the performance of our WTS algorithm, we report the detailed performance metrics and computational results in Tables 2 and 3.

"Closed Instances" and "Open Instances" represent the group of instances whose optimality has been proven and the group of instances that have not yet been proven to be optimal, respectively. Column "Count" describes the number of instances in which the algorithm found the optimal or best solutions. Column "Instance" gives the names of instances. Columns "n" and "p" report the numbers of vertices and centers, respectively. Column "Hit" shows the success rate of the algorithm for reaching the best result it has obtained under the given time limit. Column "CPU" indicates the normalized average total CPU execution time in seconds. If an algorithm did not report the result, the corresponding itemwill be marked with a hyphen "-". Column "LB" gives the lower bound of each instance calculated by the MIP solver CPLEX in 7 days [Londe et al., 2021]. Column "UB" represents the best upper bound obtained by these algorithms, and the numbers in bold stand for the best known results.

Table 2 compares the overall results on all instances obtained by GRASP-VNS [L´opez-S´anchez et al., 2019], CPLEX [Londe et al., 2021], BRKGA (BRKGA-NLS, BRKGA-FI, and BRKFA-BI) [Londe et al., 2021], ILS [Londe et al., 2021] and our WTS. In particular, GRASPVNS only reports 132 instances from pmed1 to pmed8 (except pmed5). We can observe that WTS obtains the optimal results and the best known results for all the instances. Specifically, on the closed instances, our WTS solves 10 more instances to optimality comparing to the previous best known results, which are obtained for the first time, while on the open instances, it improves the previous best known results on 3 instances. Moreover, the average CPU time of WTS is less than 0.2 seconds. WTS is very stable since it obtains the reported results with 100% success rate in all 20 independent runs, while no reference algorithm can achieve this overall performance.

Table 3 shows the details of the 12 improved instances. Note that BRKGA and ILS did not report detailed results for these instances. Except for instances pmed6 200 30 and pmed16 350 40, the upper bound obtained by our WTS algorithm is equal to the lower bound obtained by CPLEX, i.e., LB = UB. This demonstrates that the lower bounds of these instances are their optimal solutions, while CPLEX cannot prove the optimality for these 10 instances within the time limit of 7 days. In summary, these statistics indicate that WTS is highly effective, efficient, and stable for solving the pNCP.

4.3 Importance of WTS Strategies

In order to evaluate the merits of the solution-based tabu search strategy, the weighting technique, and the hierarchical penalty strategy for neighborhood evaluation, we compare WTS with three alternative versions.

1. WLS: Disable the solution-based tabu strategy.
2. TS: Deactivate the weighting technique.
3. WTS′ : Disable the hierarchical penalty strategy, i.e., the objective function uses Eq. (10).

We conduct experiments on four representative instances (pmed16 350 40, pmed21 450 90, pmed29 550 60, and pmed32 650 70). Figure 2 depicts the evolution of the objective value gaps of WTS, WLS, TS, and WTS′ as the search proceeds. Each point (x, y) on the curves represents that the gap between the cost of the current solution and the best known one is y at x iterations.

From Figure 2, one can observe that WTS is the best one of the four versions. Specifically, WLS can obtain similar solution quality to WTS in a short time. However, the convergence rate of WLS decreases constantly in the later search stage. The reason for this phenomenon might be that WLS keeps searching in a loop for previously visited solutions. Furthermore, we can observe that both TS and WTS′ optimize the solution much slower than WTS. In addition, WLS, TS, andWTS′ cannot reach the best solution in 105 iterations, while WTS can obtain the best solution within about 104 iterations. These observations indicate that the solution-based tabu strategy, the weighting technique, and the hierarchical penalty strategy all play important roles, and they are essential for the effectiveness and efficiency of WTS.

5 Conclusion

In this paper, we study a variant of the p-center problem, named the p-next center problem (pNCP), and propose a weighting-based tabu search (WTS) algorithm for solving this challenging NP-hard problem. We solve the original pNCP by decomposing and solving a series of decision subproblems. Tested on 413 widely used benchmark instances and compared with several state-of-the-art algorithms in the literature, WTS shows its advantages in solving the pNCP. It improves the best known results on 12 instances, while matching the best records in the literature for all the remaining ones. Thus, we can conclude that WTS is highly competitive in terms of the effectiveness and the efficiency.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grant 72101094.