# A rounding and clustering-based exact algorithm for the p-center problem

Zacharie Alès, Cristian Duran-Mateluna, Sourour Elloumi

► To cite this version:

Zacharie Alès, Cristian Duran-Mateluna, Sourour Elloumi. A rounding and clustering-based exact algorithm for the p-center problem. Computers and Operations Research, 2024, 183, pp.107185. 10.1016/j.cor.2025.107185. hal-04821165

HAL Id: hal-04821165

https://hal.science/hal-04821165v1

Submitted on 17 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
A rounding and clustering-based exact algorithm for the
p-center problem

Zacharie ALES$^{a,b}$, Cristian DURAN MATELUNA$^{c,d}$, Sourour ELLOUMI$^{a,b}$

<sup>a</sup>: UMA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France.

<sup>b</sup>: CEDRIC, Conservatoire National des Arts et Métiers, 75003 Paris, France.

<sup>c</sup>: University of Santiago of Chile (USACH), Faculty of Engineering, Industrial Engineering Department, Chile

<sup>d</sup>: University of Santiago of Chile (USACH), Faculty of Engineering, Program for the Development of Sustainable Production Systems (PDSPS), Chile

**Abstract**

The *p*-center problem consists of selecting *p* facilities from a set of possible sites and allocating a set of clients to them in such a way that the maximum distance between a client and the facility to which it is allocated is minimized. This paper proposes a new scalable exact solution algorithm based on client clustering and an iterative distance rounding procedure. The client clustering enables to initialize and update a subset of clients for which the *p*-center problem is iteratively solved. The rounding drastically reduces the number of distinct distances considered at each iteration. Our algorithm is tested on 396 benchmark instances with up to 1.9 million clients and facilities. Our results show that our approach outperforms existing methods run on the same computer except when *p* is smaller than 5. In this case, however, we optimally solve all instances in less than 2 minutes on average.

*Keywords:* discrete location; *p*-center problem; clustering algorithm; distance rounding

# 1 Introduction

The *p*-center problem (PCP) is a fundamental problem in location science ([16]). It consists of choosing *p* facility sites to minimize the maximum distance from a client to its nearest open facility. The *p* open facilities are called *centers*. Formally the problem is defined as follows. We consider a set of *N* clients, a set of *M* potential facility sites to open, and their corresponding index set *N* = {1, ..., *N*} and *M* = {1, ..., *M*}. Let *d*<sub>*ij*</sub> be the distance between client *i* and site *j*, and *p* be an integer. The objective is to find a set *S* ⊆ *M* such that |*S*| ≤ *p* and the value *z* = max<sub>*i*∈*N*</sub> (min<sub>*j*∈*S*</sub> *d*<sub>*ij*</sub>) is minimized. This distance *z* is called the *radius* of solution *S*.
The *p*-center is a NP-hard problem ([15]) and is applied in various contexts, including facility location, emergency service planning, telecommunications network design, healthcare planning, and supply chain management. It allows to decide the best places for centers, such as warehouses or hospitals, to minimize distances or response times ([3]). There are a variety of exact and approximate methods for this classical problem and its variants.

In this paper, we focus on the exact solution of large-scale instances. We combine existing ideas from the literature and introduce original ones in order to design an efficient solution algorithm.

Since the optimal radius is equal to one of the distances $d_{ij}$, the first idea consists in considering a subset of the distances and deducing bounds from the result. This idea already appears in the binary search algorithm introduced in [18] and [9], where for a chosen distance $\delta$ one can solve a set-covering problem to deduce whether the optimal radius is more than $\delta$ or not. This allows to iteratively divide the number of possible distances by two. This idea was for example extended by [4] by only considering a subset of distances. We extend this idea by iteratively rounding down the distances with increasing precision, thus allowing us to drastically divide the number of distances at each iteration.

The second idea is related to the management of the distances {$d_{ij}$}$_{i \in \mathcal{N}, j \in \mathcal{M}}$. For medium-scale instances, storing all distances is impractical, while for large instances, even computing them is prohibitive. Consequently, these instances are generally represented by the coordinates of the clients and sites. State-of-the-art methods have developed strategies to obtain an optimal solution while computing only a subset of the $N \times M$ distances. One of these strategies consists in only considering a subset of clients $\mathcal{R}$, that we denote by *representatives* in this article. Optimally solving the (PCP) for $\mathcal{R}$ enables to obtain a lower bound and an upper bound on the optimal radius. In [6] and [5] the authors start with an initial set of representatives and keep increasing it until the lower bound and the upper bound coincide. This strategy was also recently used and improved by [8] to develop one of the most effective methods of solving large-scale instances. In our approach, we use a clustering method to efficiently initialize and update the set of representatives. For this purpose, we leverage the clients coordinates to compute distances between pairs of clients and the barycenter of a set of clients. Note that in the experiments, we use the Euclidean distance, but any other norm could be considered. Furthermore, following previous works, we also round the distances to the next integer. Nevertheless, our approach can compute an optimal solution even for fractional distances.
## Contributions

We introduce a new exact method for the (PCP) which combines the two aforementioned ideas of considering a subset of distances and clients. Our main contributions are the following:

1. We obtain strong bounds on the optimal solution value by rounding down the distances to their nearest multiple of $\beta^\alpha$ and solving the associated subproblem, where $\beta > 1$ is a predefined base. The integer $\alpha$ is decremented at each iteration until an optimal solution of the problem is obtained. This process drastically reduces the number of distinct distances considered at any step. In A we assess the sensitivity of the method to the value of $\beta$ and observe that our method remains robust across the different choices tested. In the remainder of this manuscript we assume that $\beta = 10$.

2. We use a clustering method to initialize the subset $\mathcal{R}$ of representatives of the clients. We additionally take advantage of the partition obtained by introducing the concept of cluster quadrants which enables us to update set $\mathcal{R}$ by only adding a limited number of additional representatives.

3. We design an efficient procedure to manage the sites dominations. The list of dominated sites must be updated either when the representatives set $\mathcal{R}$ is modified or when the bounds on the optimal solution value are improved. In both cases, our procedure significantly reduces the number of sites to be compared.

4. We conduct an extensive computational analysis where we first assess the impact of different original features of our method. Then, we compare our algorithm to two state-of-the art methods ([8, 13]), and we show that our method outperforms the other two methods when $p$ is not very small (i.e., greater than 5). We consider substantially harder instances than those previously used to evaluate these two exact methods by increasing the maximal value of $p$ from 30 to 200 and the maximal number of clients and sites from one million to 1.9 million, all methods being run on the same computer.

The remainder of this article is organized as follows. Section 2 reviews the related literature. Section 3 presents our exact solution method. Illustrations and extensive computational results are presented in Sections 4 and 5. Finally, some conclusions are drawn in Section 6.
## Notations

A solution of the *p*-center problem is represented by a vector $y \in \{0, 1\}^M$ such that $y_j$ is equal to 1 if and only if site $j \in \mathcal{M}$ is opened in solution $y$. The set of opened sites in solution $y$ is denoted by $S(y) \subset \mathcal{M}$. For a given subset of clients $\mathcal{R} \subseteq \mathcal{N}$, $rad_{\mathcal{R}}(y)$ corresponds to the radius of $y$ when only the clients in $\mathcal{R}$ are considered (i.e., $rad_{\mathcal{R}}(y) = \max_{i \in \mathcal{R}} \min_{j \in S(y)} d_{ij}$). Thus, the radius of solution $y$ is $rad_{\mathcal{N}}(y)$. The *allocation distance* of a client $i \in \mathcal{N}$ in solution $y$ is $\min_{j \in S(y)} d_{ij}$. A client is said to be *covered within distance* $d \in \mathbb{R}$ by solution $y$ if its allocation distance is at most $d$. Let $K$ be the number of distinct distances in $\{d_{ij}\}_{i \in \mathcal{N}, j \in \mathcal{M}}$ and let $\mathcal{D} = \{D^1, \dots, D^K\}$ be these distances in increasing order.

For a given $\alpha \in \mathbb{Z}$, let $d_{ij}^{\alpha}$ be distance $d_{ij}$ rounded down to its nearest multiple of $10^{\alpha}$. We denote by $(PCP_{\alpha})$ the *p*-center problem in which distances $d_{ij}^{\alpha}$ are considered instead of $d_{ij}$. Let $rad_{\mathcal{R}}^{\alpha}(y)$ be the radius of $y$ for clients $\mathcal{R}$ in $(PCP_{\alpha})$ (i.e., $rad_{\mathcal{R}}^{\alpha}(y) = \max_{i \in \mathcal{R}} \min_{j \in S(y)} d_{ij}^{\alpha}$).

# 2 Literature review

Several MILP formulations of the *p*-center problem have been presented and discussed in the literature. The first and classical MILP formulation for the (PCP) was proposed in [10] that consider integer variables to represent the location and allocation decisions, and one real variable to represent the radius. [12] introduce an alternative formulation which relaxation bound dominates that of [10]. Other formulations are presented in [4], [2], and [13]. We refer the reader to [11] for an extensive computational comparison of these formulations. The direct solution of any of these formulation in a standard solver fails at optimally solving instances with several thousands of clients and sites in reasonable time.

A key idea to solve the *p*-center was introduced by [18] and presented in [9]. For a given distance $D^k \in \mathcal{D}$, determining if the optimal radius of the (PCP) is at most equal to $D^k$ can be achieved through the solution of a minimum-cardinality set cover problem denoted by $(SC_k)$. In $(SC_k)$, one constraint is associated to each client $i \in \mathcal{N}$. This constraint imposes that at least one of the sites $\{j \in \mathcal{M} \mid d_{ij} \le D^k\}$ must be opened. Consequently, if there exists a solution with at most *p* open sites, then the optimal solution value of (PCP) is at most $D^k$. From this observation, [18] proposes to solve the (PCP) by incrementing *k* until $(SC_k)$ leads to a solution of value at most *p*. As presented in Algorithm 1, this can be improved by solving the set cover problems through a binary search on $\mathcal{D}$ ([14]). To reduce the size of $\mathcal{D}$, an initial solution is additionally computed at
Step 3.

[10] also considered a binary search approach, but solved maximum coverage problems instead of set cover problems. [1] propose implementation enhancements that reduce the number of iterations. [12] improved Algorithm 1 by first relaxing the integrity of the set cover variables in Step 8. This enables to quickly reach a very strong lower bound called $LB^*$. Different heuristics are also introduced to obtain an upper bound called $UB^*$. These bounds significantly reduce the possible values of the optimal radius index in $\mathcal{D}$ and thus the number of binary search iterations.

**Algorithm 1: BinSearch:** solve (PCP) by binary search [14]

**Data:** clients $\mathcal{N}$, sites $\mathcal{M}$, $p \in \mathbb{N}^*$

**Result:** An optimal solution $y$ of (PCP) and its radius

1 Compute $\{d_{ij}\}_{i \in \mathcal{N}, j \in \mathcal{M}}$

2 Compute $\mathcal{D} = \{D^1, ..., D^K\}$ the sorted distinct distances in $\{d_{ij}\}_{i \in \mathcal{N}, j \in \mathcal{M}}$

3 $y \leftarrow$ heuristic solution

4 $\bar{k} \leftarrow$ distance index such that $rad_{\mathcal{N}}(y) = D^{\bar{k}}$

5 $lb, ub \leftarrow 1, \bar{k}$

6 **while** $lb < ub$ **do**

7 $k \leftarrow \lfloor \frac{lb+ub}{2} \rfloor$

8 $\hat{y}, opt \leftarrow$ solve $(SC_k)$

9 **if** $(SC_k)$ is feasible and $opt \le p$ **then**

10 $ub \leftarrow k$

11 $y \leftarrow \hat{y}$

12 **else**

13 $lb \leftarrow k + 1$

14 **return** $y, rad_{\mathcal{N}}(y)$

In another approach, [6] and [5] presented iterative algorithms in which only a subset of clients $\mathcal{R}$ is considered. Let $y_{\mathcal{R}}$ be the optimal solution of the (PCP) in which only clients $\mathcal{R}$ are considered. They showed that $y_{\mathcal{R}}$ is optimal for the original problem if and only if $rad_{\mathcal{R}}(y) = rad_{\mathcal{N}}(y)$. Otherwise, they add a constant number of clients to subset $\mathcal{R}$.

Finally, to the best of our knowledge, the two best exact methods also consider a progressively increasing subset of clients to solve the (PCP). [13] explore a Benders decomposition of the classic formulation of [10]. They propose a method based on a specialized branch-and-cut that considers a lifting procedure for the Benders cuts. On the other hand, [8] leverage domination rules and a local search heuristic. These two methods can both solve instances with up to one million clients for
small values of $p$ (i.e., $p \le 30$).

We now provide further details on the approach from [8] as it shares similarities with the algorithm that we propose.

## An existing row generation algorithm

The method proposed in [8] is summarized in Algorithm 2. In this algorithm, the initial subset of clients $\mathcal{R}$ is obtained in Step 1 by randomly selecting $p + 2$ clients, solving the (PCP) on these clients and repeating the procedure until the radius of the worst solution obtained is not increased for a given number of iterations. The $p + 2$ clients leading to the worst radius then constitute the initial set $\mathcal{R}$.

**Algorithm 2:** Row generation algorithm from [8]

**Data:** clients $\mathcal{N}$, sites $\mathcal{M}$, $p \in \mathbb{N}^*$

**Result:** An optimal solution $y$ of the (PCP)

1. $\mathcal{R} \leftarrow$ Initial subset of $\mathcal{N}$

2. $LB, UB \leftarrow -\infty, +\infty$

3. **while** $LB < UB$ **do**

4. &nbsp;&nbsp;&nbsp;&nbsp; $\hat{y}, rad_{\mathcal{R}}(\hat{y}) \leftarrow \text{BinSearch}(\mathcal{R})$

5. &nbsp;&nbsp;&nbsp;&nbsp; $LB \leftarrow rad_{\mathcal{R}}(\hat{y})$

6. &nbsp;&nbsp;&nbsp;&nbsp; **if** $rad_{\mathcal{N}}(\hat{y}) < UB$ **then**

7. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $y, UB \leftarrow \hat{y}, rad_{\mathcal{N}}(\hat{y})$

8. &nbsp;&nbsp;&nbsp;&nbsp; **if** $LB < UB$ **then**

9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; /* Increase $\mathcal{R}$ */
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\{y^1, ..., y^S\} \leftarrow$ alternative solutions obtained by local search around $\hat{y}$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Update $\mathcal{R}$ by deducing from $\{\hat{y}, y^1, ..., y^S\}$ new clients to add

11. **return** $y, rad_{\mathcal{N}}(y)$

At each iteration of the while loop, the $p$-center problem $(PCP_R)$ which only considers the clients in $\mathcal{R}$ in solved (Step 4). To that end, the authors propose a variation of Algorithm 1. In this variant, starting with $t = 0$, the set covers $(SC_{2t})$ are iteratively solved by incrementing $t$ until a solution of value at most $p$ is obtained. The last tested distance thus constitutes an upper bound on the optimal radius that enables to reduce the number of distance indices considered in the classical binary search algorithm. This provides a solution $\hat{y} \in \{0, 1\}^M$. Since only a subset $\mathcal{R} \subset \mathcal{N}$ is considered, the radius $rad_R(\hat{y})$ constitutes a lower bound $LB$ on the optimal radius (Step 5). The radius $rad_N(\hat{y})$,
constitutes an upper bound which may enable to decrease *UB* (Step 7).

To speed up the resolution of the set covers within the binary search, the authors also identified the
non-dominated sites.

**Definition 1 ([7])** A site $j_1 \in \mathcal{M}$ is said to be dominated by another site $j_2 \in \mathcal{M}$ relatively to a client subset $\mathcal{R}$ if $d_{i,j_1} \ge d_{i,j_2} \quad \forall i \in \mathcal{R}$ and either $\exists i \in \mathcal{R} \ d_{i,j_1} > d_{i,j_2}$ or $j_1 < j_2$.

We can observe that with this definition, two different sites cannot dominate each other.

**Proposition 1 ([7])** There exists an optimal solution of $(PCP_{\mathcal{R}})$ in which no dominated site is opened.

From this proposition, we deduce that solving $(PCP_{\mathcal{R}})$ for the non-dominated sites is the same as
solving $(PCP_{\mathcal{R}})$ for sites $\mathcal{M}$.

If the optimality is not reached (i.e., $LB < UB$), the set of considered clients $\mathcal{R}$ is increased. This step can have a significant impact on the solution time. Indeed, on the one hand, adding too much clients will slow down the solution of the ($PCP$) at each subsequent iteration. On the other hand, adding to few or non-relevant clients may increase the number of remaining iterations required to reach an optimal solution.

To add a relevant set of clients, the authors first solve maximal coverage problems through a local search heuristic in order to generate a set of alternative solutions $Y = \{y^1, ..., y^S\}$. In each of these problems, the objective is to find a solution $y^h$ which enables to maximize the number of clients in $\mathcal{N} \setminus \mathcal{R}$ covered within distance $LB$ (i.e., the optimal radius of $(PCP_{\mathcal{R}})$) while ensuring that $rad_{\mathcal{R}}(y^h) = rad_{\mathcal{R}}(\hat{y})$. The heuristic starts by randomly selecting a small set of clients $C \subset \mathcal{N} \setminus \mathcal{R}$. Solution $y^h$, which is initially equal to $\hat{y}$, is then iteratively modified by:

1. randomly replacing a given number of sites in $y^h$; and then

2. iteratively applying the site replacement that maximizes the number of clients in $C$ covered within distance $LB$ until a local minima is reached.

If solution $y^h$ is the best currently known solution of the maximum coverage problem, it is added to
$Y$. Moreover, if $y^h$ enables to cover all clients in set $C$ within distance $LB$, new clients are added
to $C$. Throughout the algorithm, it is ensured that $rad_{\mathcal{R}}(y^h) = rad_{\mathcal{R}}(\hat{y})$ by preventing any move
that would increase the maximal radius of the clients in $\mathcal{R}$. The algorithm stops when the size of $Y$
remains the same for a given number of consecutive iterations. Finally, set $\mathcal{R}$ is updated through
a completion heuristic that identifies a subset of clients of $N \setminus R$ that together would lead to an increase of the radius for all the solutions $\{\hat{y}, y^1, ..., y^S\}$ (Step 10).

Throughout the algorithm, it is ensured that $rad_R(y^h) = rad_R(\hat{y})$ by preventing any move that would increase the maximal radius of the clients in $R$.

We now present a new efficient exact algorithm based on a clustering of the client that takes also advantage of an iterative distance rounding procedure.

# 3 An exact rounding and clustering-based algorithm for (PCP)

This section is dedicated to the presentation of our exact solution method. A high-level view of it is described in Algorithm 3. We start by the clustering phase that partitions the set of clients into a given number $k \ge p$ of clusters. Then, the subset $R$ of representatives is initialized with the clusters medoids. The medoid of a cluster corresponds to the element of the cluster which is the closest to its barycenter. The lower bound $LB$ is initialized 0 and a feasible solution $y$ is chosen using $R$ to initialize the upper bound $UB$. These initialization Steps 1–4 are detailed in Subsection 3.1.

Then in the remaining Steps 5–11, we perform our iterative rounding. The decreasing parameter $\alpha$ gives the precision of the rounding. It is initialized from the number of digits in $UB$. At each iteration, the crucial step is the exact solution of problem $(PCP_\alpha)$, i.e., the $p$-center problem with the rounding precision $10^\alpha$. This is detailed in Algorithm 4 and in Subsection 3.2. At each iteration of this algorithm, the set of representatives $R$ is enlarged by adding, within each cluster, uncovered clients for different $(PCP)$ solutions. These solutions include the current solution as well as additional ones obtained by local search. We also introduce the notion of cluster quadrants in order to efficiently select the clients added to $R$. The update of the subset of clients $R$ is detailed in Section 3.3. Another important aspect of our method is the elimination of dominated sites for the set of representatives (see Definition 1) within the solution procedure of each rounded problem $(PCP_\alpha)$. We put a focus on our novel ideas on this aspect in Subsection 3.4.

We now present our new algorithm and its main features in Subsection 3.1 to 3.4.

## 3.1 The initialization steps

In Algorithm 3, the initial subset of clients is not selected randomly, as in Algorithm 2, but through the $k$-means algorithm ([17]) (Step 1). This clustering heuristic provides a partition of the clients into $k$ clusters. The first advantage of using this clustering approach is that it enables to consider a
subset of representative clients $\mathcal{R}$ of $\mathcal{N}$ by initializing $\mathcal{R}$ with the $k$-medoids of the partition (Step 2).

The second advantage is that each non-representative client in $\mathcal{N}\setminus\mathcal{R}$ is now allocated to a cluster which we will leverage to update $\mathcal{R}$. The number of clusters $k$ considered is a parameter of our method. Similarly to [8], we observed empirically that the value $p+2$ leads to good results in average.

**Algorithm 3:** Function `solveByRounding` that optimally solves the (PCP) with an increasingly precise rounding of the distances

**Data:** clients $\mathcal{N}$, sites $\mathcal{M}$, $p \in \mathbb{N}^*$, $k \ge p$

**Result:** An optimal solution $y$ of the (PCP)

1. Clustering: use the $k$-means heuristic to obtain a partition $\mathcal{C}$ of $\mathcal{N}$ into $k$ clusters

2. $\mathcal{R} \leftarrow$ the $k$ medoids of clustering $\mathcal{C}$

3. $y \leftarrow$ build a feasible solution from $\mathcal{R}$

4. $LB, UB \leftarrow 0, \text{rad}_{\mathcal{N}}(y)$

5. $\alpha \leftarrow$ (number of digits in $UB$) - 1 // Used to round the distances following Definition 2

6. **while** $LB < UB$ **do**

7. &nbsp;&nbsp;&nbsp;&nbsp;`ŷ`, `rad_{\mathcal{N}}^{\alpha}(ŷ)` $\leftarrow$ `solvePCP_{\alpha}(\mathcal{N}, \mathcal{M}, p, \mathcal{R}, \mathcal{C}, LB, UB, \alpha)$ by Algorithm ??

8. &nbsp;&nbsp;&nbsp;&nbsp;$LB \leftarrow rad_{\mathcal{N}}^{\alpha}(ŷ)$

9. &nbsp;&nbsp;&nbsp;&nbsp;**if** $rad_{\mathcal{N}}(ŷ) < UB$ **then**

10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y, UB \leftarrow ŷ, rad_{\mathcal{N}}(ŷ)$

11. &nbsp;&nbsp;&nbsp;&nbsp;`α \leftarrow α - 1`

12. **return** Solution $y$

## 3.2 Exact solution of the $p$-center problem for the rounding precision $10^{\alpha}$

The bottleneck of Algorithm 3 is the solution of the (PCP) subproblems through binary search (Step ??). Since the number of binary search iteration is in $O(\log(|\mathcal{D}|))$, we reduce the size of $\mathcal{D}$ by iteratively considering an increasingly precise rounding of the distances.

More precisely, for a given rounding factor $\alpha \in \mathbb{Z}$, each distance $d_{ij}$ is rounded down to the nearest multiple of $10^{\alpha}$.

**Definition 2 (Distance rounding)** For a given client $i \in \mathcal{N}$, a given site $j \in \mathcal{M}$, a rounding factor $\alpha \in \mathbb{N}$, together with a lower bound $LB$ and an upper bound $UB$ on the optimal value of (PCP)

$$d_{ij}^{\alpha} \triangleq \min \left( \max \left( LB, 10^{\alpha} \left\lfloor \frac{d_{ij}}{10^{\alpha}} \right\rfloor \right), UB + 1 \right).$$
The distance $d_{ij}^{\alpha}$ is imposed to be in the interval $[LB, UB + 1]$ in order to further reduce $|D|$.

The rounding factor $\alpha$ is initially equal to the number of digits in the radius of a heuristic solution $y$ minus one (Step 5). It is then decremented at each iteration thus enabling the (PCP) to be solved with a more precise rounding at the next iteration (Step 11). Note that when only integer distances are considered we have $d_{ij}^{0} = d_{ij}$ which ensures that an optimal solution is obtained for $\alpha \ge 0$.

At Step ??, we use function `solvePCP_α` to compute an optimal solution of the (PCP) problem ($PCP_α$) in which each distance $d_{ij}$ is replaced by $d_{ij}^{\alpha}$. Note that this function could be replaced by any (PCP) solution method (e.g. Algorithm 2) provided that the distances can be rounded down whenever they are computed. Function `solvePCP_α` takes as an input both clients sets $\mathcal{R}$ and $\mathcal{N}$. It considers $\mathcal{R}$ as the initial set of clients and update it by adding clients from $\mathcal{N} \setminus \mathcal{R}$ until an optimal solution $\hat{y}$ of ($PCP_α$) is obtained.

The following property ensures that optimally solving ($PCP_α$) enables to obtain a lower bound on the optimal radius of the original problem (Step 8).

**Property 1** For any solution $\hat{y}$, $rad_{\mathcal{N}}^{\alpha}(\hat{y}) \le rad_{\mathcal{N}}(\hat{y})$.

**Proof.** Since $d_{(ij)}^{\alpha} \le d_{ij}$ for all $ij \in \mathcal{N} \times \mathcal{M}$, then $rad_{\mathcal{N}}^{\alpha} = \max_{i \in \mathcal{N}} \min_{j \in S(\hat{y})} d_{ij}^{\alpha} \le \max_{i \in \mathcal{N}} \min_{j \in S(\hat{y})} d_{ij} = rad_{\mathcal{N}}(\hat{y})$. $\square$

<figure><img src="image_2.png" /></figure>

After obtaining the lower bound $rad_{\mathcal{N}}^{\alpha}(\hat{y})$, we compute the radius of solution $\hat{y}$ without rounding down the distances (Step 9) which constitutes an upper bound. If it is currently best known, we update $UB$ and $\hat{y}$ (Step 10). Note that since the radius of $\hat{y}$ on the rounded down distances is equal to $LB$, the new value of $UB$ is at most $LB + 10^{\alpha}$. Thus, throughout its iterations, this algorithm provides an increasingly tight guarantee on the quality of its bounds. Moreover, at each iteration, there can only be at most 10 distinct distances $d_{ij}^{\alpha}$.

Note that we only cluster the clients once at Step 1 of Algorithm 3. Since the clients added to $\mathcal{R}$ within function `solvePCP_α` are likely to be relevant for the subsequent iterations, they remain in $\mathcal{R}$ until the end of the algorithm.

We now present function `solvePCP_α`, detailed in Algorithm 4, that takes advantage of the initial client clustering to update the subset of clients $\mathcal{R}$. This function starts by computing for each client $i \in \mathcal{R}$ and each site $j \in \mathcal{M}$ the distance $d_{ij}^{\alpha}$.
(PCP) relaxation

In Step 12 of Algorithm 4, an optimal solution $\hat{y}$ of an integer ($PCP_\alpha$) for clients $\mathcal{R}$ is computed. Since this is the most time consuming step of our algorithm, we instead solve the linear relaxation of ($PCP_\alpha$) whenever possible. The integer ($PCP_\alpha$) is only considered if the linear relaxation solution $y^r$ does not enable to deduce any new client to add to $\mathcal{R}$.
**Algorithm 4:** Function `solvePCP<sub>α</sub>` that optimally solves (*PCP*<sub>α</sub>).

**Data:** clients $\mathcal{N}$, sites $\mathcal{M}$, $p \in \mathbb{N}^*$, $\mathcal{R} \subset \mathcal{N}$, $\mathcal{C}$ a partition of $\mathcal{N}$, $LB \in \mathbb{R}$, $UB \in \mathbb{R}$, $\alpha \in \mathbb{N}$

**Result:** An optimal solution $y$ of ($PCP_\alpha$)

1 Compute $\{d_{ij}^\alpha\}_{i \in \mathcal{R}, j \in \mathcal{N}}$

2 $LB^\alpha, UB^\alpha \leftarrow LB, UB$

3 $J_R^\alpha \leftarrow$ non-dominated sites for $\mathcal{R}$

4 `addedClients` $\leftarrow$ true

5 $y \leftarrow \emptyset$

6 **while** $LB^\alpha < UB^\alpha$ **do**

7 &nbsp;&nbsp;&nbsp;**if** `addedClients` **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`LB`<sup>r</sup>, `y`<sup>r</sup> $\leftarrow$ solve the relaxed (PCP<sub>α</sub>) for $\mathcal{R}$ and $J_R^\alpha$ by relaxed binary search
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`LB`<sup>α</sup> $\leftarrow$ `LB`<sup>r</sup>
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ŷ $\leftarrow$ compute from `y`<sup>r</sup> an integer solution
11 &nbsp;&nbsp;&nbsp;**else**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ŷ, rad<sub>\*R</sub>(ŷ) $\leftarrow$ solve the integer (PCP<sub>α</sub>) for $\mathcal{R}$ and $J_R^\alpha$ by binary search
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`LB`<sup>α</sup> $\leftarrow$ rad<sub>\*R</sub>(ŷ)
14 &nbsp;&nbsp;&nbsp;**if** $LB^\alpha < UB^\alpha$ **then**
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;{$y^1, ..., y^S$} $\leftarrow$ alternative solutions obtained by local search around ŷ
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $y_0 \in \{ŷ, y^1, ..., y^S\}$ **do**
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** rad<sub>\*N</sub>(y<sub>0</sub>) < $UB^\alpha$ **then**
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;y, $UB^\alpha$ $\leftarrow$ $y_0$, rad<sub>\*N</sub>(y<sub>0</sub>)
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`addedClients` $\leftarrow$ Update $\mathcal{R}$ by adding at most one client for each quadrant of each cluster in $\mathcal{C}$ and each solution $\{ŷ, y^1, ..., y^S\}$
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $LB^\alpha$ or $UB^\alpha$ has been improved **then**
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Clip the distances in the interval [$LB^\alpha, UB^\alpha + 1$] and update the non-dominated sites $J_R^\alpha$
22 &nbsp;&nbsp;&nbsp;**return** y, rad<sub>\*N</sub>(y)
To that end, we define a boolean parameter *addedClients* equal to true if and only if clients have been added to *R* at the last iteration. The variable *addedClients* is initially set to true (Step 4). Thus, in the first iteration the linear relaxation of (*PCP*) for clients *R* and sites *J*<sub>*R*</sub><sup>*α*</sup> is solved (Step 8). This variable is then set to false only if no client is added to *R* at Step 19. In that case, the integer (*PCP*<sub>*α*</sub>) is solved at the next iteration (Step 12).

At Step 15 of Algorithm 4, we use the same local search approach as in Algorithm 2 from [8] to deduce from an integer solution $\hat{y}$, a set of alternative solutions $\{y^1, \dots, y^S\}$. However, when we solve the linear relaxation of a ($PCP$) at Step 10, we compute a fractional solution $y^r$. To build an integer solution $y$ from $y^r$, we open the $p$ sites which values are maximal in $y^r$ (Step 10).

We observed empirically that the value $LB$ obtained when solving ($PCP_\alpha$) at Steps 10 or 12 is often the same for two consecutive iterations. Indeed, adding clients often does not lead to an increase of the optimal radius. Consequently, at the first iteration of the binary search Algorithm 1, we set $k$ to the value 1 instead of $\lfloor \frac{lb+ub}{2} \rfloor$ (Step 7). As a consequence, $D^k$ is equal to $LB$ and if the solution of the set cover at Step 8 does not require more than $p$ sites, we directly know that the optimal radius of ($PCP_\alpha$) is not increased since the last iteration. This proved to be more efficient than the variation of the binary search considered in [8] and described in Section 14.

### 3.3 Updating set $\mathcal{R}$ through the clusters quadrants and local search

The choice of the clients added to $\mathcal{R}$ at each iteration has a crucial impact on our algorithm. Indeed, increasing $\mathcal{R}$ will tend to slow down the solution of the ($PCP$) which would be an incentive not to add too many clients. However, adding too few clients will make the lower bound $LB$ increase more slowly, thus increasing the number of iterations.

To find a relevant compromise, we can first observe that only adding to $\mathcal{R}$ clients which are covered within $LB$ in the current solution $y$ cannot lead to an improvement of $LB$ at the next iteration (since the same solution $y$ would still lead to a radius of value $LB$). Thus, we focus on the clients of $\mathcal{N} \setminus \mathcal{R}$ which are not covered in $y$ within distance $LB$. Let $\mathcal{N}'(y)$ be this set of clients.

Since the $k$-means algorithm provides a coherent partition of the clients, we independently search in each cluster for clients to add to set $\mathcal{R}$. For each cluster, we want to add clients which are as far away from their medoid as possible (so that they are more likely to increase the radius) and which are also as far away from each other as possible (to make it harder to cover them without covering other points in the cluster). Consequently, as represented in Figure 1a, we divide the clients of each cluster in four *quadrants* according to their position in relation to that of their medoid. Figure 1b
represents the solution $y$ in which three sites are opened. The clients represented by a white circle correspond to the set $\mathcal{N}'(y)$. We can see in the left cluster that a client can be covered within distance $LB$ by a different site than the medoid of its cluster. For each quadrant, we add to $\mathcal{R}$ the client of $\mathcal{N}'(y)$ which is the furthest away from the representative of the cluster, if any. For example, in Figure 1c, three clients are removed from two of the clusters and only one client is removed for the leftmost cluster. We apply this process to all solutions $\hat{y}, y^1, ..., y^S$ (Step 19).

Note that the construction of quadrants as presented above is only possible when the customers are represented by two-dimensional coordinates. Nevertheless, the quadrants of a cluster $C$ can still be built in a more general context by performing a clustering of the clients in $C$. Each cluster thus obtained constitutes a quadrant of $C$.

## 3.4 Improving the computation of the sites dominations

In Algorithm 4, whenever the lower bound $LB$ or the upper bound $UB$ is improved, each distance $d_{ij}^{\alpha}$ is updated according to Definition 2. This enables to both reduce the number of distances considered in the binary search algorithm and potentially increase the number of dominated sites.

To reduce the computational time, we first observe that between two consecutive iterations, the set of dominated sites $J_{\mathcal{R}}^{\alpha}$ only needs to be updated if either:

1. one of the bounds $LB$ and/or $UB$ has been improved;

2. the number of considered clients $\mathcal{R}$ has been increased.

Moreover, in both cases, we do not need to test the domination of all possible couples of sites. More precisely, when the first case occurs, we only need to check if the non-dominated sites become dominated. Indeed, improving the bounds reduce the number of distinct distances, thus a dominated site remains dominated. Furthermore, we do not need to check if a non-dominated site is dominated by a dominated site.

In the second case, we only need to check if dominated sites become non-dominated. Indeed, adding clients to $\mathcal{R}$ cannot make a non-dominated site dominated – except in the case where two sites had the same distance to all clients before $\mathcal{R}$ was updated. To improve the computation of the dominances in that case, we define the vector $dom \in \mathbb{N}^M$ such that $dom_j$ is equal to 0 if site $j$ is not dominated and to the lowest index of a site that dominates $j$ otherwise. If $dom_j \neq 0$, we can avoid testing if site $j$ is dominated by sites $\{1, ..., dom_j - 1\}$. Indeed, if all these sites were not dominating $j$ before the addition of the clients, they can not dominate it after either.
<figure><img src="image_3.png" /><figcaption>(a) Three clusters and their medoids (▲) which are included in $\mathcal{R}$. The other clients in the cluster (●) are not included in $\mathcal{R}$. The clusters are represented by the convex hull of their clients and their quadrants are delimited by dotted lines.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) The squares represent the opened sites in solution $y$. The white circles correspond to the clients $\mathcal{N}'(y)$ which are not covered within distance $LB$ in solution $y$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) The three clusters after their update. Set $\mathcal{R}$ now contains 10 clients.</figcaption></figure>

<caption>Figure 1: Example of an update of a set $\mathcal{R}$ and three clusters.</caption>
Following the suggestions of a reviewer, we tested three variants of our algorithm which did not enable us to significantly improve the performances. We hereafter detail the three variants. At Step 17 of Algorithm 4 the radius of each solution is computed according to distances $d^\alpha$ which may enable to improve $UB^\alpha$ at Step 18. The first variant was to additionally compute the radius of the solutions for the original distances $d$ and update the global bound $UB$ of Algorithm 3 if it is improved. The second variant was to limit the solution time of $\mathbf{solvePCP}_\alpha$ to a quarter of the remaining time for all iterations except the last one. If the time limit is reached, $\mathbf{solvePCP}_\alpha$ returns non-tight bounds leading to more distinct distances in the next iteration. Finally, we tried to modify Step 10 of Algorithm 4 which initially only computes a feasible integer solution from the fractional solution $y^r$ by setting to 1 the $p$ largest components of $y^r$ and to 0 the others. In this variant, we compute an additional feasible solution by considering the sites in decreasing order of their component in vector $y^r$ and opening the $p$ first ones that reduce the radius of the solution.

# 4 An illustrative example

We illustrate our algorithm by considering the TSPLib ([19]) instance nu3496 which contains the coordinates of 3496 cities in Nicaragua. To obtain a (PCP) instance, each city both constitutes a client and a candidate site. The number of sites to open is $p = 5$.

Figure 2 represents the main steps of our algorithm. In Figure 2a, the partitioning of the clients in 7 clusters is depicted. Each cluster is represented by a distinct colour and by its star-shaped medoid. For this instance the parameter $\alpha$ is initially set to 3. Consequently, in the first iteration we solve the problem ($PCP_{\alpha=3}$) in which all distances are rounded down to their closest multiple of 1000. During the resolution of ($PCP_3$), 55 new clients are added to the set of representatives. These clients are represented by triangles in Figure 2b and the red squares correspond to the opened sites in the optimal solution of ($PCP_3$). The optimal radius is 1000, which provides a lower bound on the radius of the original (PCP) problem. Moreover, when considering the true distances (i.e., without rounding them down), this same solution leads to a radius of 1603 which thus constitutes an upper bound.

In the next iteration, we solve ($PCP_{\alpha=2}$) where the distances are rounded down to their closest multiple of 100. Thanks to the bounds [1000, 1603] obtained in the previous iteration, we know that there are only 7 possible values for the optimal radius of ($PCP_2$): 1000, 1100, 1200, 1300, 1400, 1500, or 1600. More generally, in our algorithm, we necessarily solve ($PCP_\alpha$) problems containing
at most 10 distances. As represented in Figure 2c, we only need to add 13 additional clients to the set of representatives to compute the optimal solution of (PCP₂). We deduce from this solution that the optimal radius of the original problem is in the interval [1100, 1197].

The third iteration, represented in Figure 2d, adds 11 clients to the set of representatives to solve (PCP₁) and ensures that the optimal radius is between 1120 and 1123. Finally, an optimal radius of 1123 is obtained during the last iteration in which no new client is added to the set of representatives. In the end, our algorithm optimally solves the problem by only considering 86 clients out of the 3496 contained in the instance.

<figure><img src="image_6.png" /><figcaption>(a) The 7 clusters and their medoids (initial representatives).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) It. 1: optimal solution when distances are rounded down to multiples of 10³. This iterations adds 55 new representatives and provides the bounds [1000, 1603].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) It. 2: optimal solution when distances are rounded down to multiples of 10². This iterations adds 13 new representatives and provide the bounds [1100, 1197].</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) It. 3: optimal solution when distances are rounded down to multiples of 10¹. This iterations adds 11 new representatives and provide the bounds [1120, 1123].</figcaption></figure>

<caption>Figure 2: Illustration of our algorithm on instance **nu3496** with $p = 5$. The medoids of the clusters are represented by stars, the other clients in the representative set by triangles and the sites opened in the optimal solution of each iteration by red squares.</caption>
# 5 Computational experiments

Following previous works of [8] and [13], we consider in this experiment 44 instances from the TSPLib ([19]) containing the coordinates of 1600 to 1904711 cities. For each TSP instance, a p-center instance is obtained by considering a client and a potential site at the location of each city (i.e. $|N| = |M|$). The values of $p$ considered are $p \in \{2, 3, 5, 10, 15, 20, 50, 100, 200\}$. Note that $p$ is at most 30 in the results presented in [8] and [13]. A time limit of 3 hours is considered. For all instances, the distances between the clients and the sites are obtained by computing the Euclidean distance between their coordinates rounded to its closest integer. For our experiments, we consider an Intel XEON E5-2643 processor 3,4 GHz, with 24 threads and 252 GB of RAM.

First, in Section 5.1 we evaluate the efficiency of three of our algorithmic features: the site dominations, the local search, and the iterative rounding of the distances. This enables to identify that the use of site dominations on the largest instance is too time consuming and should be deactivated. Then, in Section 5.2 we compare our approach to the two best approaches from the litterature: the iterative algorithm from [8] and the branch-and-cut from [13]. We kindly thank the authors of these articles for providing us with the resources to execute their methods on our computer, thus allowing a fair comparison. All MILPs are solved using IBM ILOG CPLEX 22.1.

## 5.1 Features comparison of our method

We first compare the performances of our method in four different configurations:

1. *All features activated:* We use Algorithm 3 to solve the (PCP) which calls Algorithm 4 at each iteration to solve (PCP$_{α}$).

2. *Without the local search:* We do not use the local search algorithm at Step 15 of Algorithm 4. Consequently, at Step 19 only solution $\hat{y}$ is considered to add clients to $\mathcal{R}$ instead of $\{\hat{y}, y^1, \dots, y^S\}$.

3. *Without dominations:* The non-dominated sites $J_{\mathcal{R}}^{\alpha}$ are not computed in Step 3 of Algorithm 4, nor updated in Step 3. Thus, (PCP$_{α}$) is always solved by considering all the sites instead of only the non-dominated ones.

4. *Without rounding the distances:* We directly solve (PCP$_{α=0}$) through Algorithm 4 without using Algorithm 3. As a consequence, the number of distinct distances between clients and sites may be high, which is likely to increase significantly the number of set cover iterations.
The results are summarized in Figure 3 through performance profiles of CPU times and gap. Each curve of the performance profile corresponds to the results of one configuration. Each point (x, y) of a curve indicates that for x percents of the instances, the configuration was no more than y times worse than the best method. In particular, y = 1 corresponds to the proportion of instances on which the configuration is the best.

<figure><img src="image_10.png" /><figcaption>Figure 3: Performance profiles of our algorithm with different features deactivated. A higher curve on the graph corresponds to a more effective method.</figcaption></figure>

The results in Figure 3 are split in three categories depending on the size of the instance: $N \le 10000$, $N \in [10000, 230000]$, and $N \ge 230000$. For each size, we represent two performance profiles: one for the solution time, one for the final gap obtained. Note that for a given configuration, the resolution time of an instance is only accounted for if it is below the time limit.

We can first see in Figure 3 that the configuration in which the distance are not rounded is always the worst one, both in terms of gap and of resolution time. This is even more significant for the largest instances. The use of Algorithm 3 which iteratively rounds down the distances more and more precisely, thus seems to be the most efficient feature of our algorithm. We can then see that the computation of the site dominations enables to speed up the resolution except for the largest instances. For these instances, most of the available time is spent computing the dominations instead
of solving the set cover problems which eventually leads to poor upper and lower bounds. As a consequence, in the next section we do not compute the dominations for these instances. Finally, we can see that removing the local search tends to slow down the resolution but enables to improve the number of instances solved to optimality. Since our objective is to improve the solution of hard instances, we keep this feature in the following.

In our method, the distances at each iteration are rounded down to their nearest multiple of $\beta^\alpha$
with $\beta = 10$. Following the recommendation of a reviewer, we experimentally showed that other
values of $\beta$ do not significantly impact the quality of the results. This experiment is detailed in A.

## 5.2 Comparison to state-of-the-art methods

In this section, we compare our algorithm with the iterative algorithm from [8] and the branch-and-cut
from [13]. As mentioned before, thanks to their authors, we were able to run the three methods on
the same computer and execute their code for larger values of $p$ with the same time limit than us.
The results are summarized by the performance profiles in Figure 4 which splits the instances in
three categories depending on the number of sites opened: $p \in \{2, 3, 5, 10\}$, $p \in \{20, 30, 50\}$, and
$p \in \{100, 200\}$. To give a more detailed view of the results, we also present Tables 1 and 2 in which
the results are averaged for each value of $p$ and $N$, respectively. In these tables, a value is in bold if
it is strictly better than the one obtained by the two other methods.

From Table 1, we can see that the three methods solve all instances to optimality within 3 hours for
the smallest values of $p \in \{2, 3, 5\}$. For these instances, the branch-and-cut from [13] is the fastest.
However, for all the other values of $p$, our algorithm provides the best results in terms of number
of instances solved to optimality, final gap and resolution time. In particular, the resolution time
reduction when applying our method for $p \ge 10$ is very significant: from 3 times faster for $p = 10$ to
60 times faster for $p = 200$. This is clearly highlighted by the two top left profiles of Figure 4.
<figure><img src="image_11.png" /><figcaption>Figure 4: Performance profiles of our algorithm compared with two state-of-the-art methods. A higher curve on the graph corresponds to a more effective method.</figcaption></figure>

<table><caption>Table 1: Average results obtained for each method and each value of parameter *p*. A value is in **bold** if it is strictly better than that of the other methods. All methods run on the same computer with a time limit of 3 hours. *: only computed on instances optimally solved by the three methods.</caption><thead><tr><th rowspan="2">p</th><th rowspan="2">#inst.</th><th colspan="3">Our algorithm</th><th colspan="3">[8]</th><th colspan="3">[13]</th></tr><tr><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th></tr></thead><tbody><tr><td>2</td><td>44</td><td>44</td><td>0.0</td><td>63</td><td>44</td><td>0.0</td><td>22</td><td>44</td><td>0.0</td><td><b>9</b></td></tr><tr><td>3</td><td>44</td><td>44</td><td>0.0</td><td>74</td><td>44</td><td>0.0</td><td>53</td><td>44</td><td>0.0</td><td><b>22</b></td></tr><tr><td>5</td><td>44</td><td>44</td><td>0.0</td><td>144</td><td>44</td><td>0.0</td><td>252</td><td>44</td><td>0.0</td><td><b>124</b></td></tr><tr><td>10</td><td>44</td><td><b>43</b></td><td><b>0.3</b></td><td><b>153</b></td><td>41</td><td>0.6</td><td>465</td><td>33</td><td>1.2</td><td>1098</td></tr><tr><td>15</td><td>44</td><td><b>35</b></td><td><b>0.9</b></td><td><b>21</b></td><td>30</td><td>2.5</td><td>182</td><td>19</td><td>6.0</td><td>509</td></tr><tr><td>20</td><td>44</td><td><b>30</b></td><td><b>1.9</b></td><td><b>107</b></td><td>20</td><td>5.7</td><td>542</td><td>13</td><td>10.9</td><td>1384</td></tr><tr><td>50</td><td>44</td><td><b>13</b></td><td><b>10.7</b></td><td><b>19</b></td><td>7</td><td>16.6</td><td>812</td><td>5</td><td>22.0</td><td>1251</td></tr><tr><td>100</td><td>44</td><td><b>9</b></td><td><b>24.4</b></td><td><b>19</b></td><td>6</td><td>25.6</td><td>660</td><td>3</td><td>25.8</td><td>677</td></tr><tr><td>200</td><td>44</td><td><b>7</b></td><td><b>27.4</b></td><td><b>18</b></td><td>5</td><td>32.4</td><td>1094</td><td>4</td><td>31.5</td><td>1779</td></tr></tbody><tfoot><tr><td>sum/avg.</td><td>396</td><td><b>269</b></td><td><b>7.3</b></td><td><b>94</b></td><td>241</td><td>9.2</td><td>242</td><td>209</td><td>10.8</td><td>406</td></tr></tfoot></table>

In Table 2, the bold values confirm that for most of the instances, we outperform the two state-of-the-art methods. The method from [8] is able to obtain better gaps on 9 instances and the branch-and-cut from [13] is faster than our method on 5 instances. Four of these five instances are the largest ones. This coincides with the observation in Section 5.1 that computing the domination for these instances was too time consuming. Nevertheless, the number of instances solved to optimality
and the final gaps obtained by our method remains better in average for these large instances than
that of the branch-and-cut.
<table><caption>Table 2: Average results obtained for each method and each instance size. A value is in **bold** if it is strictly better than that of the other methods. All methods run on the same computer with a time limit of 3 hours. *: only computed on instances optimally solved by the three methods.</caption><thead><tr><th rowspan="2">Name</th><th rowspan="2">N</th><th rowspan="2">#inst.</th><th colspan="3">Our algorithm</th><th colspan="3">[8]</th><th colspan="3">[13]</th></tr><tr><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th><th>#opt. gap[%]</th><th>time[s]*</th></tr></thead><tbody><tr><td>rw</td><td>1621</td><td>9</td><td>9</td><td>0.0</td><td><b>3</b></td><td>9</td><td>0.0</td><td>49</td><td>9</td><td>0.0</td><td>14</td></tr><tr><td>u</td><td>1817</td><td>9</td><td><b>9</b></td><td><b>0.0</b></td><td><b>25</b></td><td>8</td><td>1.2</td><td>770</td><td>8</td><td>0.2</td><td>1694</td></tr><tr><td>rl</td><td>1889</td><td>9</td><td><b>9</b></td><td><b>0.0</b></td><td><b>7</b></td><td>8</td><td>0.3</td><td>142</td><td>6</td><td>0.6</td><td>1045</td></tr><tr><td>mu</td><td>1979</td><td>9</td><td>9</td><td>0.0</td><td><b>4</b></td><td>9</td><td>0.0</td><td>122</td><td>9</td><td>0.0</td><td>227</td></tr><tr><td>pr</td><td>2392</td><td>9</td><td><b>8</b></td><td><b>0.6</b></td><td><b>4</b></td><td>6</td><td>4.1</td><td>25</td><td>4</td><td>4.5</td><td>196</td></tr><tr><td>pcb</td><td>3038</td><td>9</td><td><b>6</b></td><td><b>8.2</b></td><td><b>4</b></td><td>5</td><td>8.3</td><td>20</td><td>4</td><td>9.5</td><td>31</td></tr><tr><td>nu</td><td>3496</td><td>9</td><td><b>9</b></td><td><b>0.0</b></td><td><b>7</b></td><td>7</td><td>2.2</td><td>103</td><td>7</td><td>2.7</td><td>778</td></tr><tr><td>ca</td><td>4663</td><td>9</td><td><b>9</b></td><td><b>0.0</b></td><td><b>1</b></td><td>8</td><td>1.0</td><td>14</td><td>6</td><td>1.9</td><td>4</td></tr><tr><td>rl</td><td>5915</td><td>9</td><td><b>6</b></td><td><b>6.6</b></td><td><b>1</b></td><td>5</td><td>9.0</td><td>5</td><td>4</td><td>12.2</td><td>3</td></tr><tr><td>rl</td><td>5934</td><td>9</td><td>6</td><td><b>5.6</b></td><td><b>17</b></td><td>6</td><td>7.5</td><td>360</td><td>4</td><td>9.5</td><td>658</td></tr><tr><td>tz</td><td>6117</td><td>9</td><td><b>7</b></td><td><b>1.5</b></td><td><b>6</b></td><td>6</td><td>6.1</td><td>31</td><td>4</td><td>8.9</td><td>19</td></tr><tr><td>eg</td><td>7146</td><td>9</td><td>8</td><td><b>0.3</b></td><td>4</td><td>8</td><td>2.4</td><td>28</td><td>7</td><td>5.3</td><td>4</td></tr><tr><td>pla</td><td>7397</td><td>9</td><td><b>7</b></td><td><b>2.2</b></td><td><b>4</b></td><td>6</td><td>5.5</td><td>51</td><td>5</td><td>8.8</td><td>35</td></tr><tr><td>ym</td><td>7663</td><td>9</td><td><b>7</b></td><td>6.5</td><td><b>4</b></td><td>6</td><td><b>4.3</b></td><td>31</td><td>6</td><td>8.0</td><td>12</td></tr><tr><td>pm</td><td>8079</td><td>9</td><td><b>8</b></td><td><b>0.8</b></td><td><b>10</b></td><td>7</td><td>3.0</td><td>478</td><td>7</td><td>5.3</td><td>152</td></tr><tr><td>ei</td><td>8246</td><td>9</td><td><b>6</b></td><td>10.7</td><td><b>7</b></td><td>5</td><td><b>9.7</b></td><td>46</td><td>4</td><td>10.6</td><td>30</td></tr><tr><td>ar</td><td>9152</td><td>9</td><td>6</td><td><b>2.7</b></td><td><b>210</b></td><td>6</td><td>8.8</td><td>691</td><td>6</td><td>8.6</td><td>2377</td></tr><tr><td>ja</td><td>9847</td><td>9</td><td>7</td><td><b>3.7</b></td><td><b>3</b></td><td>7</td><td>4.1</td><td>16</td><td>6</td><td>8.0</td><td>4</td></tr><tr><td>gr</td><td>9882</td><td>9</td><td><b>7</b></td><td>7.1</td><td><b>6</b></td><td>6</td><td><b>4.6</b></td><td>62</td><td>6</td><td>8.7</td><td>23</td></tr><tr><td>kz</td><td>9976</td><td>9</td><td>6</td><td><b>2.1</b></td><td><b>11</b></td><td>6</td><td>7.6</td><td>96</td><td>5</td><td>10.6</td><td>31</td></tr><tr><td>fi</td><td>10639</td><td>9</td><td><b>6</b></td><td><b>5.3</b></td><td>9</td><td>5</td><td>10.6</td><td>21</td><td>4</td><td>10.7</td><td><b>8</b></td></tr><tr><td>rl</td><td>11849</td><td>9</td><td><b>5</b></td><td><b>6.1</b></td><td><b>17</b></td><td>4</td><td>10.6</td><td>66</td><td>4</td><td>12.0</td><td>767</td></tr><tr><td>usa</td><td>13509</td><td>9</td><td><b>6</b></td><td><b>6.8</b></td><td><b>8</b></td><td>5</td><td>11.6</td><td>87</td><td>4</td><td>12.1</td><td>158</td></tr><tr><td>brd</td><td>14051</td><td>9</td><td><b>6</b></td><td>11.5</td><td><b>14</b></td><td>5</td><td><b>10.1</b></td><td>39</td><td>4</td><td>11.7</td><td>2353</td></tr><tr><td>mo</td><td>14185</td><td>9</td><td>6</td><td><b>4.9</b></td><td><b>6</b></td><td>6</td><td>8.0</td><td>38</td><td>4</td><td>11.8</td><td>14</td></tr><tr><td>ho</td><td>14473</td><td>9</td><td>6</td><td>8.6</td><td><b>6</b></td><td>6</td><td><b>7.2</b></td><td>23</td><td>4</td><td>11.3</td><td>29</td></tr><tr><td>d</td><td>15112</td><td>9</td><td><b>5</b></td><td>18.4</td><td><b>7</b></td><td>4</td><td><b>11.8</b></td><td>20</td><td>4</td><td>14.0</td><td>18</td></tr><tr><td>it</td><td>16862</td><td>9</td><td>6</td><td><b>5.6</b></td><td><b>14</b></td><td>6</td><td>7.2</td><td>106</td><td>6</td><td>10.1</td><td>50</td></tr><tr><td>d</td><td>18512</td><td>9</td><td><b>5</b></td><td><b>8.8</b></td><td><b>7</b></td><td>4</td><td>12.6</td><td>21</td><td>3</td><td>13.1</td><td>17</td></tr><tr><td>vm</td><td>22775</td><td>9</td><td>6</td><td><b>5.8</b></td><td><b>8</b></td><td>6</td><td>8.0</td><td>61</td><td>5</td><td>10.8</td><td>12</td></tr><tr><td>sw</td><td>24978</td><td>9</td><td><b>6</b></td><td><b>5.7</b></td><td><b>13</b></td><td>5</td><td>10.5</td><td>36</td><td>4</td><td>16.4</td><td>22</td></tr><tr><td>fyg</td><td>28534</td><td>9</td><td>4</td><td><b>11.8</b></td><td><b>100</b></td><td>4</td><td>13.6</td><td>867</td><td>4</td><td>15.8</td><td>1156</td></tr><tr><td>bm</td><td>33708</td><td>9</td><td><b>6</b></td><td><b>6.2</b></td><td><b>32</b></td><td>5</td><td>10.0</td><td>134</td><td>4</td><td>12.2</td><td>238</td></tr><tr><td>pla</td><td>33810</td><td>9</td><td>4</td><td><b>9.3</b></td><td><b>25</b></td><td>4</td><td>15.3</td><td>283</td><td>3</td><td>16.1</td><td>70</td></tr><tr><td>bby</td><td>34656</td><td>9</td><td>4</td><td><b>11.5</b></td><td><b>130</b></td><td>4</td><td>12.8</td><td>737</td><td>4</td><td>15.4</td><td>484</td></tr><tr><td>pba</td><td>38478</td><td>9</td><td>4</td><td><b>12.4</b></td><td><b>16</b></td><td>4</td><td>13.5</td><td>51</td><td>3</td><td>16.5</td><td>24</td></tr><tr><td>ch</td><td>71009</td><td>9</td><td>5</td><td>14.5</td><td><b>94</b></td><td>5</td><td><b>12.9</b></td><td>327</td><td>4</td><td>14.2</td><td>509</td></tr><tr><td>pla</td><td>85900</td><td>9</td><td><b>4</b></td><td><b>9.3</b></td><td><b>76</b></td><td>3</td><td>18.9</td><td>840</td><td>3</td><td>18.0</td><td>306</td></tr><tr><td>sra</td><td>104815</td><td>9</td><td><b>5</b></td><td>16.4</td><td>61</td><td>4</td><td><b>12.8</b></td><td>96</td><td>3</td><td>17.1</td><td><b>46</b></td></tr><tr><td>usa</td><td>115475</td><td>9</td><td>4</td><td><b>13.9</b></td><td><b>240</b></td><td>4</td><td>15.8</td><td>1048</td><td>4</td><td>15.8</td><td>1678</td></tr><tr><td>ara</td><td>238025</td><td>9</td><td>4</td><td>19.5</td><td>288</td><td>4</td><td><b>16.1</b></td><td>230</td><td>3</td><td>16.3</td><td><b>124</b></td></tr><tr><td>lra</td><td>498378</td><td>9</td><td><b>4</b></td><td><b>12.5</b></td><td>579</td><td>3</td><td>19.4</td><td>922</td><td>3</td><td>12.6</td><td><b>230</b></td></tr><tr><td>lrb</td><td>744710</td><td>9</td><td>3</td><td>18.2</td><td>990</td><td>3</td><td>22.3</td><td>1370</td><td>3</td><td>18.2</td><td><b>984</b></td></tr><tr><td>world</td><td>1904711</td><td>9</td><td><b>6</b></td><td><b>19.0</b></td><td>2076</td><td>4</td><td>25.3</td><td>812</td><td>5</td><td>30.3</td><td><b>272</b></td></tr></tbody><tfoot><tr><td>sum/avg.</td><td></td><td>396</td><td><b>269</b></td><td><b>7.3</b></td><td><b>94</b></td><td>241</td><td>9.2</td><td>242</td><td>209</td><td>10.8</td><td>406</td></tr></tfoot></table>
# 6 Conclusions

We introduce a new exact method to solve the *p*-center problem based on an iterative rounding of the distances which greatly limits the number of possible value for the optimal solution at each iteration. Our experimental analysis highlights that this new feature is by far the one that most improves the performances of our algorithm and allows it to scale up. Moreover, our method includes a clustering of the clients which enables to efficiently initialize the subset of representatives $\mathcal{R}$. We also introduce the notion of cluster quadrants which enables to update $\mathcal{R}$ by adding to it relevant clients which are not close to one another and not too many either. Finally, we also efficiently manage the computation of the site dominations and limit the number of MILP subproblems to solve by considering their linear relaxation as often as possible. Through extensive experiments, we compare our approaches to the iterative algorithm from [8] and the branch-and-cut from [13]. As soon as parameter *p* is greater than 5, our method significantly outperforms these two state-of-the-art algorithms in terms of resolution time, final gap and number of instances solved to optimality. In the future we will try to better handle site dominations of the largest instances and to speed up the update of the subset of clients.

# Acknowledgments

The authors acknowledge the financial support from the Programa de Cooperación Científica ECOS-ANID project number ECOS230013 and the Programa Regional STIC-AMSUD project number AMSUD240008.

We would like to thank the authors of [13] and [8] for providing us with the executable and source code of their methods, respectively.
References

[1] Salhi Al-Khedhairi, A. Enhancements to two exact algorithms for solving the vertex p-center problem. *Journal of Mathematical Modelling and Algorithms*, 4:129–147, 2005.

[2] Zacharie Ales and Sourour Elloumi. Compact MILP formulations for the p-center problem. In Jon Lee, Giovanni Rinaldi, and A. Ridha Mahjoub, editors, *Combinatorial Optimization*, pages 14–25, Cham, 2018. Springer International Publishing.

[3] Hatice Çalık, Martine Labbé, and Hande Yaman. p-Center problems. In *Location Science*, pages 51–65. Springer, 2019.

[4] Hatice Calik and Barbaros C. Tansel. Double bound method for solving the p-center location problem. *Computers & Operations Research*, 40(12):2991–2999, 2013.

[5] Doron Chen and Reuven Chen. New relaxation-based algorithms for the optimal solution of the continuous and discrete p-center problems. *Computers & Operations Research*, 36(5):1646–1655, 2009.

[6] R. Chen and G. Y. Handler. Relaxation method for the solution of the minimax location-allocation problem in euclidean space. *Naval Research Logistics*, 34(6):775–788, 1987.

[7] Richard L. Church. Symposium on location problems: In memory of Leon Cooper. *Journal of Regional Science*, 24(2):185–201, 1984.

[8] Claudio Contardo, Manuel Iori, and Raphael Kramer. A scalable exact algorithm for the vertex p-center problem. *Computers & Operations Research*, 103:211–220, 2019.

[9] M. Daskin. Center problems. In *Network and Discrete Location: Models, Algorithms, and Applications*, pages 154–197. Wiley, New York, 1995.

[10] Mark Daskin. Network and discrete location: Models, algorithms and applications. *Journal of the Operational Research Society*, 48, 01 1996.

[11] Cristian Durán Mateluna. *Exact solution methods for large-scale discrete p-facility location problems*. Theses, Institut Polytechnique de Paris, January 2024.

[12] Sourour Elloumi, Martine Labbé, and Yves Pochet. A New Formulation and Resolution Method for the p-Center Problem. *INFORMS Journal on Computing*, 16(1):84–94, 2004.
[13] Elisabeth Gaar and Markus Sinnl. A scaleable projection-based branch-and-cut algorithm for the p-center problem. *European Journal of Operational Research*, 303(1):78–98, 2022.

[14] R S Garfinkel, A W Neebe, and M R Rao. The m-Center Problem: Minimax Facility Location. *Management Science*, 23(10):1133–1142, 1977.

[15] O. Kariv and S. L. Hakimi. An algorithmic approach to network location problems. ii: The p-medians. *SIAM Journal on Applied Mathematics*, 37(3):539–560, 1979.

[16] Gilbert Laporte, Stefan Nickel, and Francisco Saldanha-da Gama. *Location Science* (2nd ed). Springer International Publishing, 2019.

[17] James MacQueen et al. Some methods for classification and analysis of multivariate observations. In *Proceedings of the fifth Berkeley symposium on mathematical statistics and probability*, 1967.

[18] Edward Minieka. The m-center problem. *SIAM Review*, 12(1):138–139, 1970.

[19] Gerhard Reinelt. TSPLIB—A Traveling Salesman Problem Library. *ORSA Journal on Computing*, 3(4):376–384, 1991.
A Rounding the distances using different bases

In our method, distances $d^\alpha$ are obtained by rounding down each distance $d$ to its closest multiple
of $\beta^\alpha$ with $\beta = 10$. In this section, we assess the impact of the choice of $\beta$ in relation to the
order of magnitude of the distances contained in the instances. For this purpose, we created
instances in which the coordinates of the clients are randomly generated in the interval $[0, M]$
with $M \in \{10^1, 10^3, 10^5, 10^7\}$ representing the *magnitude* of an instance. For each magnitude, we
generated one instance for each of the following number of clients $n \in \{10^3, 5 \times 10^3, 10^4, 1.5 \times 10^4, 2 \times$
$10^4, 5 \times 10^4, 10^5, 2 \times 10^5, 5 \times 10^5\}$. Each instance is solved for each value of $\beta \in \{2, 5, 10, 50, 100\}$, and
each number of opened sites $p \in \{2, 3, 5, 10, 15, 20\}$. Table 3 presents the average results obtained.
We observe that all the bases lead to similar results. There does not seem to be a clear relation
between the couple $(\beta, M)$ and the quality of the results.
<table><caption>Table 3: Average solution times and gaps obtained. The columns BX correspond to the results when the base $\beta$ is equal to X. A value is in **bold** if it is strictly better than that of the other methods. All bases were run with a time limit of 1 hour. *: only computed on instances optimally solved by all the methods.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">n</th><th colspan="5">Average gap (%)</th><th colspan="5">Average time (s)*</th></tr><tr><th>B2</th><th>B5</th><th>B10</th><th>B50</th><th>B100</th><th>B2</th><th>B5</th><th>B10</th><th>B50</th><th>B100</th></tr></thead><tbody><tr><td></td><td>10<sup>3</sup></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1</td><td>1</td><td>1</td><td><b>0</b></td><td>1</td></tr><tr><td></td><td>5 × 10<sup>3</sup></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td><b>42</b></td><td>74</td><td>72</td><td>56</td><td>53</td></tr><tr><td></td><td>10<sup>4</sup></td><td>8.3</td><td>8.3</td><td><b>0.0</b></td><td>8.3</td><td>8.3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td></tr><tr><td></td><td>1.5 × 10<sup>4</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td><b>0.0</b></td><td>8.3</td><td>2</td><td>1</td><td>1</td><td>1</td><td>2</td></tr><tr><td rowspan="5">10</td><td>2 × 10<sup>4</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>6</td><td>4</td><td>4</td><td><b>2</b></td><td>3</td></tr><tr><td>5 × 10<sup>4</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>14</td><td>7</td><td><b>5</b></td><td>6</td><td>6</td></tr><tr><td>10<sup>5</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>34</td><td>20</td><td>14</td><td><b>13</b></td><td>14</td></tr><tr><td>2 × 10<sup>5</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>170</td><td>94</td><td>72</td><td>57</td><td><b>52</b></td></tr><tr><td>5 × 10<sup>5</sup></td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>8.3</td><td>517</td><td>372</td><td>238</td><td>231</td><td>231</td></tr><tr><td></td><td>10<sup>3</sup></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>4</td><td>4</td><td>4</td><td>6</td><td>6</td></tr><tr><td></td><td>5 × 10<sup>3</sup></td><td>0.4</td><td>0.2</td><td><b>0.1</b></td><td>0.2</td><td>0.8</td><td><b>79</b></td><td>82</td><td>131</td><td>141</td><td>121</td></tr><tr><td></td><td>10<sup>4</sup></td><td>2.3</td><td>2.8</td><td><b>1.1</b></td><td>1.2</td><td>2.0</td><td>24</td><td>27</td><td>24</td><td>34</td><td>42</td></tr><tr><td></td><td>1.5 × 10<sup>4</sup></td><td>2.3</td><td>4.9</td><td><b>1.7</b></td><td>6.5</td><td>2.3</td><td>65</td><td><b>56</b></td><td>77</td><td>91</td><td>96</td></tr><tr><td rowspan="5">10<sup>3</sup></td><td>2 × 10<sup>4</sup></td><td>2.2</td><td>4.7</td><td><b>1.7</b></td><td>6.6</td><td>2.9</td><td>156</td><td>197</td><td><b>147</b></td><td>168</td><td>149</td></tr><tr><td>5 × 10<sup>4</sup></td><td>2.7</td><td>5.0</td><td><b>2.2</b></td><td>7.3</td><td>4.8</td><td>106</td><td><b>66</b></td><td>89</td><td>129</td><td>178</td></tr><tr><td>10<sup>5</sup></td><td>3.6</td><td>5.4</td><td><b>2.3</b></td><td>7.8</td><td>4.6</td><td><b>363</b></td><td>382</td><td>515</td><td>712</td><td>1156</td></tr><tr><td>2 × 10<sup>5</sup></td><td><b>4.0</b></td><td>5.6</td><td>4.3</td><td>12.1</td><td>12.8</td><td>442</td><td>244</td><td><b>204</b></td><td>237</td><td>321</td></tr><tr><td>5 × 10<sup>5</sup></td><td><b>4.1</b></td><td>6.4</td><td>11.2</td><td>13.9</td><td>17.4</td><td>1782</td><td>1046</td><td>1053</td><td><b>1003</b></td><td>1505</td></tr><tr><td></td><td>10<sup>3</sup></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td><b>5</b></td><td>7</td><td>9</td><td>9</td><td>14</td></tr><tr><td></td><td>5 × 10<sup>3</sup></td><td><b>0.1</b></td><td>0.4</td><td>0.4</td><td>0.3</td><td>1.0</td><td>16</td><td>15</td><td>15</td><td>25</td><td>25</td></tr><tr><td></td><td>10<sup>4</sup></td><td>2.1</td><td><b>1.3</b></td><td>1.4</td><td>3.2</td><td>1.6</td><td><b>43</b></td><td>58</td><td>47</td><td>78</td><td>106</td></tr><tr><td></td><td>1.5 × 10<sup>4</sup></td><td>1.9</td><td><b>1.3</b></td><td>1.7</td><td>5.2</td><td>2.6</td><td><b>178</b></td><td>258</td><td>225</td><td>241</td><td>368</td></tr><tr><td rowspan="5">10<sup>5</sup></td><td>2 × 10<sup>4</sup></td><td><b>1.9</b></td><td>17.8</td><td>2.0</td><td>5.2</td><td>2.7</td><td>406</td><td>361</td><td><b>242</b></td><td>346</td><td>820</td></tr><tr><td>5 × 10<sup>4</sup></td><td>3.7</td><td>17.8</td><td><b>3.0</b></td><td>5.2</td><td>3.4</td><td>189</td><td>139</td><td>149</td><td><b>97</b></td><td>160</td></tr><tr><td>10<sup>5</sup></td><td>4.1</td><td>20.7</td><td><b>3.3</b></td><td>5.5</td><td>6.4</td><td>451</td><td>381</td><td>418</td><td><b>276</b></td><td>1010</td></tr><tr><td>2 × 10<sup>5</sup></td><td>7.2</td><td>21.1</td><td><b>5.4</b></td><td>6.5</td><td>13.8</td><td>1162</td><td>626</td><td>762</td><td><b>482</b></td><td>680</td></tr><tr><td>5 × 10<sup>5</sup></td><td><b>7.9</b></td><td>21.3</td><td>10.6</td><td>9.5</td><td>17.1</td><td>3542</td><td>1632</td><td>1198</td><td><b>780</b></td><td>1475</td></tr><tr><td></td><td>10<sup>3</sup></td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>8</td><td>8</td><td>11</td><td>12</td><td>12</td></tr><tr><td></td><td>5 × 10<sup>3</sup></td><td>0.7</td><td><b>0.1</b></td><td>0.4</td><td>0.4</td><td>0.9</td><td>18</td><td>14</td><td><b>12</b></td><td>14</td><td>16</td></tr><tr><td></td><td>10<sup>4</sup></td><td>1.0</td><td>1.4</td><td>1.6</td><td>2.6</td><td>1.7</td><td><b>56</b></td><td>57</td><td>74</td><td>92</td><td>106</td></tr><tr><td></td><td>1.5 × 10<sup>4</sup></td><td>1.3</td><td><b>1.1</b></td><td>1.8</td><td>2.6</td><td>2.1</td><td>208</td><td><b>165</b></td><td>175</td><td>177</td><td>190</td></tr><tr><td rowspan="5">10<sup>7</sup></td><td>2 × 10<sup>4</sup></td><td>5.4</td><td>4.3</td><td><b>1.9</b></td><td>2.6</td><td>2.5</td><td>310</td><td>227</td><td><b>226</b></td><td>283</td><td>337</td></tr><tr><td>5 × 10<sup>4</sup></td><td>4.7</td><td>4.9</td><td><b>2.2</b></td><td>2.8</td><td>3.0</td><td>300</td><td><b>164</b></td><td>198</td><td>205</td><td>246</td></tr><tr><td>10<sup>5</sup></td><td>7.1</td><td>4.6</td><td>3.3</td><td><b>3.2</b></td><td>6.5</td><td>618</td><td>436</td><td>387</td><td><b>262</b></td><td>363</td></tr><tr><td>2 × 10<sup>5</sup></td><td>6.6</td><td>6.5</td><td>6.7</td><td><b>5.7</b></td><td>12.5</td><td>1338</td><td>1024</td><td>934</td><td>809</td><td><b>741</b></td></tr><tr><td>5 × 10<sup>5</sup></td><td>10.2</td><td>8.0</td><td>8.6</td><td><b>5.1</b></td><td>16.7</td><td>3152</td><td>2092</td><td>1451</td><td><b>794</b></td><td>948</td></tr><tr><td></td><td>avg</td><td>4.0</td><td>6.3</td><td><b>3.6</b></td><td>4.8</td><td>5.6</td><td>439</td><td>287</td><td>255</td><td><b>218</b></td><td>321</td></tr></tbody></table>