Effective Approaches to Solve P-Center Problem via Set Covering and SAT

The classic <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-center problem consists of choosing a set of <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula> vertices in an undirected graph as facilities in order to minimize the maximum distance between each client vertex and its closest facility. The problem is equivalent to covering all vertices by no more than <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula> circles with the smallest possible radius, which can be tackled by solving a series of the decision version of set covering subproblems with the same cardinality constraint (<inline-formula> <tex-math notation="LaTeX">$\le p$ </tex-math></inline-formula>) and gradually decreasing the covering radius. In this paper, we solve the <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-center problem via set covering and SAT. We first transform the <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-center problem into a series of set covering subproblems and simplify them by some reduction rules. Then, we present two kinds of encoding methods to convert them into CNF format and solve them with several state-of-the-art SAT solvers. Tested on three sets of totally 70 benchmark instances, our proposed approach can improve the previous best known results for 3 instances using the heuristic SAT solvers while proving the optimality for 59 instances using the exact SAT solvers. The computational results demonstrate the effectiveness of the proposed approach in terms of both solution quality and computational efficiency. In addition, the main advantage of our approach is twofold: The independence of the subproblems allows the problem to be solved in parallel; The approach to transform the original problem into SAT is flexible such that various state-of-the-art SAT solvers can be used.


I. INTRODUCTION
The p-center problem is a classical facility location problem which has important applications in the fields of telecommunication industry, transportation science, public services and so on. It consists of seeking the location of p facilities and assigning clients to them in order to minimize the maximum distance between a client and its nearest facility.
Hakimi [1] first introduced the absolute center problem to find a center C in a graph G such that the maximum distance from C is minimized. The absolute center problem can be considered as the p-center problem when p = 1. However, there may be multiple centers that need to locate in real applications. Therefore, Hakimi [2] extended the absolute center problem to the p-center problem for more general applications. The p-center problem has been proven to be NP-hard [3] and many algorithms including exact and heuristic methods have been proposed for solving this problem.
The associate editor coordinating the review of this manuscript and approving it for publication was Juan Liu .
Most exact methods are based on solving a finite series of related subproblems, for example, the set covering problem. Minieka [4] proposed the first set covering based approach to solve the p-center problem by decreasing the covering radius step by step until the optimal value of the set covering subproblem within the current radius is greater than p. Garfinkel et al. [5] improved the approach proposed by Minieka [4] by using a binary search technique for the selection of the covering radius to reduce the search space. Elloumi et al. [6] proposed a new integer programming formulation (ELP) which is also based on the set covering subproblem and is better than the classical formulation proposed by Daskin [7]. Calik and Tansel [8] developed new integer programming formulations and an exact algorithm based on decomposition to their models (IP) for solving the p-center problem, which provides a tighter lower bound than the approach proposed by Elloumi et al. [6]. Daskin [9] presented an algorithm where the covering radius is also searched by a binary search, but the subproblem is replaced by the maximal set covering problem which consists of maximizing the number of covered clients by no more than p facilities within the given radius. Ilhan and Pinar [10] proposed a two-phase method based on a feasibility subproblem where it checks whether it is feasible to cover all the clients by no more than p facilities within a given radius. In this method, it first solves the linear programming (LP) feasibility problems to obtain a suitable lower bound and then solves a series of IP feasibility problems starting from the lower bound. Al-Khedhairi and Salhi [11] proposed modifications (IP* and Daskin*) to the previous algorithms: They enhanced the algorithm proposed by Ilhan and Pinar [10] by reducing the number of ILP iterations needed to find the optimal solution, and modified the approach proposed by Daskin [7] with tighter initial lower and upper bounds and a more appropriate binary search method which can reduce the number of subproblems to be solved.
There are also various effective heuristic and metaheuristic algorithms for the p-center problem. Mladenovic et al. [12] presented a basic variable neighborhood search (VNS) and two tabu search (TS) heuristics for this problem. Their experiments showed that the proposed VNS and two TS algorithms are superior to previous classical heuristics, and VNS performs the best on average in terms of both the solution quality and computational time. However, TS performs slightly better for the instances with smaller p values. Pullan [13] proposed a memetic genetic algorithm (PBS), which is a population-based meta-heuristic that uses phenotype crossover and directed mutation operators to generate new starting points for a local search. For larger p-center instances, PBS is able to effectively utilize a number of computer processors. It can get high quality solutions for the classical benchmarks. Yin et al. [14] proposed a greedy randomized adaptive search procedure with path-relinking (GRASP/PR) algorithm, which combines both GRASP and path-relinking. In their algorithm, each iteration of GRASP/PR consists of the construction of a randomized greedy solution, followed by a tabu search procedure. GRASP/PR was considered to be the best heuristic for it can obtain the optimal solutions for most of the public benchmarks.
Satisfiability (SAT) is the problem of determining if there exists an assignment of boolean variables that satisfies a given Boolean formula and is the first problem proven to be NP-complete [15]. There are many literatures about solving combinatorial optimization problems based on SAT. Van [16] proposed a method for solving graph coloring problem by encoding it into SAT. They improved the previous lower bounds for the classical benchmarks such as DSJC125.5. Soh et al. [17] proposed a SAT based exact approach for solving the two-dimensional strip packing problem (2SPP). They showed that their method is competitive with the previous state-of-the-art 2SPP methods. In particular, they found better solution for instance HT08 than the best solution found by the previous exact and heuristic methods.
In this paper, we present a method for solving the p-center problem via set covering and SAT, since the SAT problem has been intensively studied in the academic society over the last decades and a lot of progress has been made. Thus, various state-of-the-art SAT solvers can be used, which is an advantage of our method. To our best knowledge, our paper establishes a bridge between the p-center problem and SAT for the first time. We first transform the problem to a finite series of the decision version of set covering subproblems with the same cardinality constraint (≤ p) and gradually decreasing covering radius. Then, we simplify these set covering subproblems by some data reduction rules which were proposed by Alber et al. [18] and encode them into CNF format with two different encoding modes. We then solve the p-center problem by solving these encoded subproblems using the state-of-the-art SAT solvers. Note that these encoded subproblems are independent and can be solved in parallel. Our method is tested on three sets of totally 70 classical benchmarks of p-center problem and it can obtain the optimal or best known solutions for all the 70 instances, where the previous best known results are improved for 3 instances. Furthermore, we can prove the optimality for 59 out of 70 instances. To the best of our knowledge, there are some instances which are proven to optimality for the first time.
The rest of this paper is organized as follows. Section II describes the formal definition of the p-center problem and its decision version problem. Section III presents the method for solving the p-center problem via set covering and SAT. Section IV reports the computational results and comparison with the state-of-the-art algorithms in the literature. Section V analyzes and discusses the importance of data reduction, how to prove the optimality for the p-center problem by SAT and the performance differences of the two encoding modes used in our algorithm. Section VI summarizes the main contribution of this work and concludes the paper.

II. PROBLEM DEFINITION
Given a complete undirected graph G = (N , E), where N is the set of nodes and E is the set of edges. The distance between any two nodes u and v is denoted by d(u, v). The p-center problem is to find a set of facilities S ⊆ N , such that |S| = p and the objective function: is minimized. The p-center problem can be transformed from an optimization problem to a series of decision subproblems. Assuming that all distinguishing lengths of edges in G are sorted in a decreasing order, i.e., d 1 > d 2 > . . . > d K . Then, we choose d i (i = 1, 2, . . . , K ) in decreasing order as the coverage radius and remove all the edges whose lengths are larger than d i . For every d i , we judge if there is a feasible solution for the p-center problem. We repeatedly choose next d i until d i is invalid, i.e., there does not exist a set S ⊆ N which satisfies both |S| = p and max v∈N min u∈S d(u, v) ≤ d i . Then, d i−1 is the best solution for the original p-center problem.
The methodology of solving an optimization problem by transforming it into a series of decision problems is a popular VOLUME 8, 2020 technique for solving challenging combinatorial optimization problems [19]- [21]. The reason lies in the fact that some optimization problems are much more difficult to solve than their corresponding decision problems. By transforming an optimization problem to a decision one, the search space of the problem becomes greatly narrowed and restricted. In addition, it also becomes easy to identify the reason to cause the infeasibility of a decision problem, which can guide the search effectively and efficiently. For these reasons, we employ this methodology to solve the p-center problem.

III. SOLVING P-CENTER PROBLEM VIA SET COVERING AND SAT
The subproblem in the decision version of the p-center problem can be transformed to the decision version of set covering problem, which can be encoded into CNF format and then solved with SAT solvers. Based on this idea, the p-center problem can be solved via set covering and SAT.

A. TRANSFORMING THE SUBPROBLEM TO SET COVERING
Ilhan and Pinar [10] first transformed the subproblem of p-center problem to the decision version of the set covering problem, which was formulated as an IP model. Based on the model, we make a small modification such that it can be encoded into CNF format. Given a complete graph G = (N , E) and a coverage radius d, used for the set covering subproblem can be obtained by removing the edges that cannot cover a node, i.e., whose length is greater than d. For each node u ∈ N , we can get a set S u ⊆ N such that for each node v ∈ S u the distance between u and v is less than or equal to d. Then, we can obtain the set F = {S 1 , S 2 , . . . , S u , . . . , S n }, where S u = {u}∪{v ∈ N |(u, v) ∈ E }. Thus, judging if there exists a set S ⊆ N , such that |S| = p and max v∈N min u∈S d(u, v) ≤ d is equivalent to judging if there exists p sets in F which can cover all the nodes, i.e., if there exists a subset F ⊆ F which satisfies N = ∪ S u ∈F S u and |F | = p. From another perspective, for each node u ∈ N , F must have at least one set that contains u. Furthermore, if F can cover all the nodes by less than p sets, we can add arbitrarily p − |F | sets to F and F can definitely cover all the nodes. The cardinality constraint |F | = p can be relaxed to |F | ≤ p. Let x c be a binary variable which means that S c belongs to F if x c = 1, we can get the following model.
In this model, constraint (1) denotes that for each node u ∈ N , there exists at least one set S c which contains u belonging to F , and constraint (2) restricts that we can choose at most p sets from F. Figure 1 shows an example of transforming a p-center problem into set covering problems with different radii.

B. ENCODING SET COVERING INTO CNF
In order to solve the subproblem with SAT solvers, we need to encode the model into CNF format where the formula is conjunction of clauses and each clause is a disjunction of literals; each literal is a propositional variable x or its negation ¬x.
For constraint (1), we can encode it with: Since the graphs of different set covering subproblems of a p-center problem are different, constraint (1) may be different for a vertex in two different subproblems. For example, there are three sets (S 1 , S 2 and S 6 ) which contain vertex 1 in Figure 1-(b). Thus, we can obtain a clause (x 1 ∨ x 2 ∨ x 6 ) for vertex 1 to ensure it to be covered, while the clause for vertex 1 in Figure 1 (2) is the cardinality constraint, which is independent of the specific graphs of both the original p-center problem and the transformed set covering problems. That is to say, it is only related to the values of n and p. We encoded the cardinality constraint in sequential counter and parallel counter encoding modes to maintain the cardinality constraint.
The sequential counter encoding mode was introduced in Sinz [22] for the n c=1 x c ≤ p cardinality constraint. The encoding mode works by encoding a sequential counter circuit ( Figure 2) that sequentially calculates the sums s i = i j=1 x j which is represented as unary number with p bits (s i,1 , . . . , s i,p ), and sets the overflow bits v i to be true if s i exceeds p (i = 1, 2, . . . , n). In addition, all of the overflow bits v i should be zero to ensure that s i ≤ p. Then, the cardinality constraint n c=1 x c ≤ p can be encoded with 2np + n − 3p − 1 clauses and n + (n − 1)p variables. For example, the problem in Figure 1 with n = 6 and p = 2 will introduce 23 clauses and 16 variables. For the sake of space limit, we show a small example with n = 3  and p = 2, where the constraint can be encoded with 8 clauses and 7 variables as follows: The clause set A s calculates the sum s 1 which is decided only by x 1 , and B s calculates the sum s 2 which is decided by both x 2 and s 1 . In addition, the clause set C s is to set the overflow bits v 2 and v 3 to be zero.
Sinz [22] also introduced another encoding mode based on a parallel counter circuit. The counter ( Figure 3) recursively splits the n-bit input bits x i (i = 1, . . . , n) into two halves where the number of true variables is respectively counted. The results of the two halves are represented by m-bit (m = log 2 n ) binary numbers and added by a standard binary adder to get the result of the current recursive layer. In addition, the (m + 1)-bit output of the uppermost layer is the result of the whole circuit. The counter is completely implemented by full-adder and half-adder, both of which can be encoded based on the well-known equations, where each half-adder (a ⊕ b) can be encoded by three clauses: and each full-adder (a ⊕ b ⊕ c) can be encoded by seven clauses: In addition, there is a comparator circuit which forces the (m + 1)-bit output value of the parallel counter to be not greater than p, and can be encoded by at most m + 1 clauses. The encoding of the n c=1 x c ≤ p constraint based on the circuit below requires at most 7n − 3 log 2 n − 6 clauses and 3n − 2 variables. For example, when n = 3 and p = 2, the constraint can be encoded as follows: The clause sets A p and B p are generated from the parallel counter circuit which has just one full-adder. The clause set C p is to ensure the cardinality to be not greater than 2 which comes from the comparator circuit.

C. DATA REDUCTION FOR SET COVERING
Before encoding the model of set covering into CNF, we can employ two reduction rules, which were introduced by Alber et al. [18], to reduce the graph in the preprocessing stage. These two rules are based on exploring local structures of the graph and try to replace them by simpler structures. It is worth mentioning that these rules are implemented with polynomial-time complexity and can ensure that the solution domain of the reduced graph is the same as the original graph for the set covering problem. So, these rules do not change the optimal solution for the set covering problem. Based on this, we can use these two rules to reduce the graph for each set covering subproblem. Both rules can determine some vertices to be centers or not, and remove these determined vertices. Given a graph G = (N , E) and a covering radius d, the induced graph G = (N , E ) used for the set covering subproblem can be obtained by removing the edges that cannot cover a node, i.e., whose length is greater than d. 2 (v) and can be fixed or removed in a more complex way. Interested readers are referred to [18] for more details.
Let N c and N d respectively denote the set of vertices which are determined to be centers, and the set of vertices which are determined not to be centers. The solution of the set covering subproblem with graph G and covering radius d contains all the vertices in N c and does not contain any vertex in N d . Note that each vertex in N d can be covered by one vertex in N c , so we need to find at most p − |N c | centers from N \(N c ∪N d ) to cover all the vertices of N \(N c ∪N d ). By doing this, all the vertices of N can be covered by a center and the number of centers is no more than p. Let G = (N , E ) and and E is obtained by removing all the edges related to the vertices in N c ∪ N d , we will solve the following set covering model via SAT solvers: Figure 4 gives an example of rule 1 based on the graph in Figure 1-(b), which shows that the neighborhood of vertex 3 is partitioned into 1 Therefore, after employing the reduction rules for the graph in Figure 4, vertex 3 is added into N c , while vertices 2, 4, 5 and 6 are added into N d . Then, we can obtain a reduced graph which only contains vertex 1 by removing N c and N d and the set covering model will become x 1 ≥ 1 and x 1 ≤ 1, which can be encoded with only one clause (x 1 ). Thus, {x 3 , x 1 } is a feasible solution for this example.

IV. COMPUTATIONAL RESULTS
In this section we report the experimental results of different SAT solvers on the classical standard benchmark instances of the p-center problem and compare the performance of these SAT solvers with the state-of-the-art algorithms.

A. SAT SOLVERS
The famous SAT Competition has been organized for 21 times to solve the satisfiability problem since 1992.
There are many state-of-the-art SAT solvers submitted to these SAT competitions. One advantage of our method is that it can utilize different SAT solvers in parallel. To solve our problem, we test 7 SAT solvers on the instances from OR-lib [23] and TSP-Lib [24] and eventually choose 3 representative solvers including Maple_CM solver [25], MapleLCMDistChronoBT solver [26] and Sparrow2Riss-2018 solver [27], [28]. The Maple_CM solver is a new conflict-driven clause learning (CDCL) solver which was developed by extending clause minimization to original clauses based on Maple_LCM [29]. The Maple_LCM solver got the first place of the main track in SAT Competition 2017 which was obtained by implementing the learnt clause minimization approach by using unit propagation based on the solver MapleCOMSPS_DRUP [30], [31].
The MapleLCMDistChronoBT is the winner solver of the main track of SAT Competition 2018. The solver was based on the SAT Competition 2017 winner, Maple_LCM_Dist [29], and was updated with chronological backtracking (configuration {T = 100, C = 4000}) based on the results in [32].
The Sparrow2Riss-2018 solver is the first prize of the random track of SAT Competition 2018 which is a combination of the solvers Sparrow [33] and Riss. Sparrow is a Stochastic Local Search (SLS) solver that uses promising variables and probability distribution based selection heuristics, while Riss is a CDCL solver which is based on MINISAT [34] search engine and GLUCOSE 2.2 [35], [36]. As SLS solvers cannot prove unsatisfiability, they combined Sparrow and Riss by first trying to solve the SAT problem with Sparrow, which is limited for the execution with 5 · 10 8 flips, and then running Riss which starts from the solution obtained by Sparrow. By doing this, the solver's overall behavior can still be deterministic. However, Gaussian Elimination and Cardinality Constraint reasoning are applied in the solver for our experiments, which means that it cannot emit proof for unsatisfiable instances.

B. PROBLEM INSTANCES AND EXPERIMENTAL PROTOCOL
We carry out computational experiments on three representative sets of instances from OR-Lib [23] and TSP-Lib [24], respectively, which were widely used by previous state-of-the-art reference algorithms.

1) OR-LIB INSTANCES
The first set of instances consists of 40 randomly generated instances with |N | ranging from 100 to 900 and p ranging from 5 to 90. The graphs in this set are not complete, so we need to calculate the shortest path between each pair of nodes to get a complete graph.

2) TSP-LIB INSTANCES
Both the second and third sets of instances are from TSP-Lib, which are real world application instances from the task of drilling holes in printed circuit boards and are usually used as the benchmarks for various routing and location problems. The two sets respectively consist of 15 instances from u1060 (|N | = 1060) and 15 instances from u1817 (|N | = 1817). For these 30 instances, p ranges from 10 to 150 and the nodes are given with two-dimensional coordinates. Thus, we need to calculate the Euclidean distance between each pair of nodes to get the whole distance matrix, i.e., a complete graph.
Our algorithm was programmed in C++ and the experiments were conducted on a Linux PC with Intel Xeon CPU E5-2609 v2 2.50 GHZ processor and 32 GB RAM with 8 cores. The time limit for each given radius of each instance was set to 4 hours. Besides, we test each instance for five times to get the average running time.

C. EXPERIMENTAL DESIGN
For each instance, we try all possible covering radii from the previous best radius reported in the literature. Due to the independence of set covering subproblems, we solve 8 (the number of cores) subproblems with different covering radii of an instance at the same time. Once a feasible solution is obtained for a certain radius, all the subproblems with larger radii are stopped and the subsequent smaller radii are tried. The computation of an instance stops if a subproblem with a radius d i−1 outputs 'SAT' and a subproblem with a radius d i outputs 'UNSAT', or when the time limit is met. For each given radius, we can get a graph where all the edges whose lengths are larger than the given radius are deleted. First, we employ the data reduction rules mentioned in Section III-C for the graph to get a simpler one. Next, we encode the instance for the given radius with a simplified graph into CNF format by sequential counter and parallel counter. So, there are two different SAT instances to solve with SAT solvers for each given radius. We finally solve these SAT problems transformed from set covering subproblems to judge if the given radius is feasible for the p-center problem by using Maple_CM, MapleLCMDistChronoBT and Sparrow2Riss-2018.
Note that a SAT instance encoding a p-center instance consists of a relatively small subset of very long clauses only containing positive literals (constraint 1) and a large subset of clauses encoding the cardinality constraints (constraint 2). This fact is not exploited by the general-purpose SAT solvers we test in this paper, but could be exploited by a SAT solver specialized to solve the p-center problem.

D. RESULTS OF SAT SOLVERS
In this section, we conduct experiments for solving the feasibility subproblems by the three SAT solvers on the three sets of 70 instances. Table 1 shows the computational results of the three solvers on instances from OR-Lib for the given radius, while Table 2 and Table 3 respectively show the results on the sets of u1060 and u1817 instances from TSP-Lib. In these three tables, columns SEQ and PAR in Maple_CM, MapleLCMDistChronnoBT and Sparrow2Riss-2018 show the average CPU time ('−' means timeout) for each given radius in the sequential counter encoding and the parallel counter encoding modes. Note that there is just one row for one instance if all the solvers can get the same minimum radius for both encoding modes in the given time limit. Otherwise, there will be multi-row results for one instance with different covering radius. From Tables 1-3, we can get the minimum radius for all the instances in different encoding modes that each solver can obtain within the time limit. The minimum radius of one instance in certain encoding mode for one solver is the best solution of the p-center problem that the solver can obtain for the instance in the encoding mode. For example, the best solution of pmed33 obtained by Maple_CM or MapleLCMDistChronnoBT in any encoding mode is 29, while Sparrow2Riss-2018 can get 28 in the SEQ encoding mode. Row Num. shows the number of instances for which the solver can solve the instances to optimality or reach the best known results and row Avg. gives the average time of solving the instances to their corresponding best results.
From Table 1, one observes that in the same encoding mode, for both SEQ and PAR, Maple_CM and MapleLCMDistChronnoBT can get the same best solution for all the 40 instances, which means that Maple_CM is consistent with MapleLCMDistChronnoBT in terms of the solution quality on these small instances in OR-Lib. Sparrow2Riss-2018 has a little worse performance in terms of the solution quality and computational time than other two solvers in PAR encoding mode, but outperforms them in SEQ encoding mode for it can obtain better results for 9 instances in the SEQ mode. Besides, Sparrow2Riss-2018 in the SEQ mode outperforms other two solvers for it has the smallest average computational time. This experiment demonstrates that both encoding mode and SAT solver have important affects on the search performance of the algorithm.
From Table 2 and Table 3, one observes that in the same encoding mode, whether SEQ or PAR, all the three solvers perform the same in terms of the solution quality. It is worth mentioning that they can obtain the optimal solutions for all the 30 large instances from TSP-Lib in the SEQ encoding mode. In addition, Maple_CM outperforms MapleLCMDistChronnoBT in terms of the computational time to reach the best results on the set of u1060 instances, while it is opposite on the set of u1817 instances. Both of them are better than Sparrow2Riss-2018 in terms of the computational time.

E. COMPARISON WITH OTHER REFERENCE ALGORITHMS
In this subsection, we compare the best results of the three SAT solvers with three exact algorithms (ELP, IP* and Daskin*) and two metaheuristic algorithms (PBS and GRASP/PR) in the literature. The total computational time of each instance using SAT solvers can be taken as the computational time of the set covering subproblem with the smallest radius that can be solved for the following reasons: 1) The set covering subproblems with different radii of an instance can be solved in parallel. 2) The initial coverage radius is set as the previous best known result for each instance, which is a common setting in solving many combinatorial optimization problems, such as graph coloring [19], maximum clique [37], etc.
3) The number of radii between the previous best known radius and the best radius that our algorithm can obtain is less than 8 for each instance. Because the reduction rules are designed for the set covering subproblems, the reference algorithms cannot directly use these reduction rules since these reference algorithms solve the original p-center problem where the objective is to minimize the coverage radius. This is also one of the advantages of our approach that existing reduction rules for set covering can be used. For this reason, we just cited the existing results in the corresponding references. Table 4 presents the best results of SAT solvers on the 40 OR-Lib instances and makes comparison with ELP, IP*, Daskin*, PBS and GRASP/PR. Table 5 and Table 6 respectively show the best results of SAT solvers on the u1060 and u1817 instances and make comparison with ELP, PBS and GRASP/PR. Columns f best and cpu respectively show the best values obtained by an algorithm and their corresponding computational time ( = 0.001). Row Num. presents the number of instances for which the optimal or the best known results can be obtained for the corresponding solver. As PBS, GRASP/PR and SAT solvers use floating point data but the reference exact algorithms use integer data, a floating point result is considered to be identical to the integer result if the rounded floating point result obtained by PBS, GRASP/PR  or SAT solvers equals the integer result obtained by the exact algorithms. There may be dozens of different edges in one unit of solution, so the solutions obtained by the reference exact algorithms are a little rough to some extent. Besides, the f best column of ELP shows the optimal values obtained by ELP, where the optimal solution that is not obtained by ELP is marked with''?''.
From Table 4, one observes that all the algorithms can obtain the optimal solutions for all the 40 instances. We can also observe that the average time of SAT solvers for the 40 OR-Lib instances is faster than IP* and Daskin* but slower than ELP, PBS and GRASP/PR. The reason might lie in the fact that these instances are relatively small and easy and transforming these instances into SAT may miss the original problem structure such that algorithms specially designed for the p-center problem, such as PBS and GRASP/PR, can solve them more easily. Table 5 shows that the results of SAT solvers can match the previous best known results for all the u1060 instances. From Table 6, one observes that solving the p-center problem via SAT solvers can improve the previous best known results for 3 cases compared with GRASP/PR (u1817 with p = 80, 130, 150). In addition, for the instances of TSP-Lib, the SAT solvers can obtain the results showed in Tables 5 and 6 with less average computational time, and the results of TSP-Lib instances obtained by SAT solvers can be proven to be optimal solutions (See Section V-B for the proof).
In sum, we present the summarized computational statistics of all the solvers on the two datasets (OR-Lib and TSP-Lib) in Table 7, where column Sparrow gives the statistics of the Sparrow2Riss-2018 solver, column Exact shows the average statistics of the two exact SAT solvers (with the SEQ encoding), and column Ref presents the statistics of the best reference algorithm GRASP/PR since it outperforms  other reference algorithms in terms of both solution quality and run time. Rows Num. and Avg. present the number of instances for which the optimal or the best known results can be obtained and the average computational time to reach the best known results for the corresponding solvers, respectively. From Table 7, one observes that for the random instances in OR-Lib Sparrow2Riss-2018, which is from the random track and is a specialized solver for random instances, can hit the lower bounds for all the instances in OR-Lib, while other exact SAT solvers hit the lower bounds for only 31 instances.  For the real world application instances in TSP-Lib, all the SAT solvers reach the best known results for all the 30 instances where the exact SAT solvers need much less time than Sparrow2Riss-2018. This phenomenon shows that the exact SAT solvers are more suitable for solving the real world application instances while it is relatively easy for the random solver Sparrow2Riss 2018 for solving the random instances. In addition, although the reference algorithms in the literature can obtain equal or better results with less computational time than the SAT solvers on the random instances, the exact SAT solvers can obtain better results with less computational time than the reference algorithms on the real world application instances, which is of practical significance.

A. IMPORTANCE OF DATA REDUCTION
In order to evaluate the effectiveness of the data reduction rules to the p-center problem approach via SAT solvers, we conduct experiments to compare the two versions of the algorithm with and without the data reduction rules. Note that the data reduction can determine some centers. Let k be the number of centers determined by the data reduction. Figure 5 shows that the ratio k/p is less than 0.2 for most instances, which means the impact of data reduction is small for these instances. But for some instances, we can improve the computational efficiency dramatically. Table 8 presents the computation results of the two versions of the algorithm in the SEQ mode for some instances when using the Maple_CM solver, where column t 1 (s) and t 2 (s) respectively show the average computational time of the listed instances with and without the data reduction rules, and Ratio shows the ratio of the number of centers determined by data reduction to p. From Table 8, one observes that pmed29 can determine 30% centers by data reduction and the time for solving pmed29 by Maple_CM with and without data reduction is 19.09s and 3764.20s when the given radius is 13 and the encoding mode is SEQ. Furthermore, data reduction can improve the solution quality for 4 instances: pemd20, pmed24, pmed30 and pmed34, which shows the importance of the data reduction rule.

B. THE PROOF OF OPTIMALITY
The p-center problem is transformed to a series of feasibility set covering subproblems with different covering radius which are solved by SAT solvers. Note that if an exact SAT solver tackles an unsatisfiable instance, the solver can prove the unsatisfiability of the instance if it outputs 'UNSAT'. Thus, it is guaranteed that there are no feasible solutions for a feasibility set covering subproblem if the output of the SAT solver is 'UNSAT'. Furthermore, we can prove that d i−1 is the optimal solution of the p-center problem if the output of SAT solver is 'SAT' for d i−1 and 'UNSAT' for d i , because this means that there does not exist a feasible solution for the set covering subproblems with the radii which are shorter than d i−1 , i.e., there does not exist a feasible solution X with f (X ) < d i−1 for the p-center problem..
Although Sparrow2Riss-2018 can get the optimal solutions for all the 70 instances, it needs to apply the Gaussian Elimination and Cardinality Constraint reasoning which cause the solver to be unable to emit proof for unsatisfiable instances. So, we have to use Maple_CM or MapleLCMDistChron-noBT, i.e., the exact SAT solvers to prove the optimality for these 70 instances. We can prove the optimality for 29 out of 40 random instances in OR-Lib, and all the 30 real application instances in TSP-Lib, thus closing these instances (the optimal solution of the 30 instances are listed by f opt in Table 2 and Table 3). Although the previous exact method ELP claimed to prove the optimality for 22 of the 30 instances in TSP-Lib, it can only give integer lower bound for the p-center problem. That is to say, the exact solution obtained by ELP is the rounded integer value of the real exact solution. To our best knowledge, our proposed SAT-based method proves the optimality for all the 30 TSP-Lib instances by giving the real exact solutions for the first time, showing its high performance for real world application instances. VOLUME 8, 2020  Figure 6 shows the scatter plots comparing SEQ and PAR encoding modes for the three SAT solvers on the sets of pmed, u1060 and u1817 instances. A point (x, y) in the plots corresponds to an instance, where x(y) represents the solving time in seconds of SEQ (PAR). A point (x, y) where x = 14400 (y = 14400) means that the instance was not solved within the time limit.

C. THE COMPARISONS OF ENCODING MODES
Although all the 70 instances have a much smaller number of clauses and variables when using the PAR encoding mode, which is very efficient in terms of used memory, the solver's performance may be worse due to the lack of propagation strength in the PAR encoding mode. From Figure 6a -Figure 6c, one observes that the PAR encoding mode outperforms SEQ on the pmed instances for most instances when using Maple_CM and MapleLCMDistChron-noBT, while it is opposite for Sparrow2Riss-2018. Figure 6d -Figure 6f show that the two encoding modes have similar performance in terms of computational efficiency on the set of u1060 instances for all the three solvers.  Figure 6i show that the SEQ encoding mode outperforms PAR in terms of computational efficiency on the set of u1817 instances for all the three solvers.

VI. CONCLUSION
In this paper, we present a new method for solving the p-center problem via set covering and SAT, which is the first SAT-based method for solving the p-center problem to our knowledge. Our proposed approach has the advantage of allowing the problem to be solved in parallel due to the independence of the set covering subproblems. The method consists of transforming the problem to a finite series of feasibility set covering subproblems, the simplification of these feasibility set covering subproblems, encoding them into CNF format with different encoding modes and then solving the encoded subproblems using state-of-the-art SAT solvers. The encoding modes that we use include sequential counter and parallel counter, and the SAT solvers consist of Maple_CM, MapleLCMDistChronoBT and Sparrow2Riss-2018.
We compare the computational results of 70 well-known benchmark instances for different encoding modes when using different SAT solvers. We find that Sparrow2 Riss-2018 is the best solver for solving these three sets of 70 instances when using the sequential counter encoding mode, for it can obtain the best solutions for all the instances. However, it also has the worst time performance when using the parallel counter mode, which indicates the importance of choosing an appropriate encoding mode.
This work demonstrates that SAT-based approaches provide a competitive alternative for solving the p-center problem. In fact, our SAT based approach improves the previous best known results for 3 instances and equals the best known ones for the remaining instances, meaning that our SAT-based approach is better than the previous approaches in terms of the best solution quality. Considering the notorious NP-hardness of the p-center problem, we believe this is significant.
In addition, this work establishes, for the first time to our best knowledge, a bridge between the SAT and p-center problems, opening promising perspectives. In fact, the results presented in the paper are obtained by only using general-purpose SAT solvers. On the one hand, SAT solving is a very fast developing field and our results make it clear that any progress in that field can be beneficial to solve the p-center problem. On the other hand, SAT solvers could be specialized for the p-center problem, by taking its particular features into account, which would greatly improve the performance of SAT solvers for the p-center problem.
It is noteworthy that our proposed SAT-based method proves the optimality for all the 30 TSP-Lib instances by giving the real exact solutions for the first time and thus closes these instances, showing its high performance for real world application instances. Combined with the proof by the the exact algorithms (ELP and IP*) in the literature on the optimality of solutions for the instances from OR-Lib, we can conclude that our proposed algorithm can obtain the optimal solutions for all the 70 instances. The study in this work inspires us that similar method can be applied for solving other challenging combinatorial optimization problems.
[37] Q. Wu  CHUMIN LI received the B.S. and Ph.D. degrees in computer science from the University of Technology of Compiegne, France, in 1985 and 1990, respectively. He is currently a Professor of computer science with the University of Picardie Jules Verne. His research interests include the practical resolution of NP-hard problems, including SAT, CSP, MaxSAT, MinSAT, MaxClique, and GCP. He is particularly interested in the intrinsical relationships between these problems. One of his research directions is to find and exploit these relationships to solve them. A recent example is the exploitation of the relationships between MaxSAT and MaxClique to solve MaxClique.