Michele Barbato

Francisco Canas
email: francisco.m.p.canas@gmail.com

Luís Gouveia

Pierre Pesneau

Node based compact formulations for the Hamiltonian p-Median Problem

Keywords: Combinatorial optimization, Integer linear programming, Polyhedral theory, Valid inequalities, Hamiltonian p-median problem, Location Routing

In this paper, we introduce, study and analyse several classes of compact formulations for the symmetric Hamiltonian p-Median Problem (HpMP). Given a positive integer p and a weighted complete undirected graph G = (V, E) with weights on the edges, the HpMP on G is to find a minimum weight set of p elementary cycles partitioning the vertices of G. The advantage of developing compact formulations is that they can be readily used in combination with off-the-shelf optimization software, unlike other types of formulations possibly involving the use of exponentially sized sets of variables or constraints. The main part of the paper focuses on compact formulations for eliminating solutions with less than p cycles. Such formulations are less well known and studied than formulations which prevent solutions with more than p cycles. The proposed formulations are based on a common motivation, that is, the formulations contain variables that assign labels to nodes, and prevent less than p cycles by stating that different depots must have different labels and that nodes in the same cycle must have the same label. We introduce and study aggregated formulations (which consider integer variables that represent the label of the node) and disaggregated formulations (which consider binary variables that assign each node to a given label). The aggregated models are new. The disaggregated formulations are not, although in all of them new enhancements have been included to make them more competitive with the aggregated models. The two main conclusions of this study are: i) in the context of compact formulations, it is worth looking at the models with integer node variables, which have a smaller size. Despite their weaker LP relaxation bounds, the fewer variables and constraints lead to faster integer resolution, especially when solving instances with more than 50 nodes; ii) the best of our compact models exhibit a performance that, overall, is comparable to that of the best methods known for the HpMP (including branch-and-cut algorithms), solving to optimality instances with up to 226 nodes within 1 hour. This corroborates our message that the knowledge of the inequalities for preventing less than p cycles is much less well understood.

Introduction

In this paper, we propose, study and analyse a several classes of compact formulations for the symmetric Hamiltonian p-Median Problem (HpMP). Given a positive integer p and a weighted complete undirected graph G = (V, E) with a weight d ij associated to each edge {i, j} ∈ E, the HpMP on G is to find a minimum-weight partition in p elementary cycles of the vertices of G (the weight of the p cycles is the sum of the weights of the edges in those cycles). If p = 1, the HpMP is the classical symmetric Traveling Salesman Problem (TSP) showing that the HpMP is NP-hard if p is part of the input. If p is not fixed, the Content and Contribution of the Paper. Formulations that are polynomial in the number of variables and constraints have the advantage that they can be readily used in combination with off-the-shelf optimization software, unlike other types of formulations possibly involving the use of exponentially sized sets of variables or constraints. While typically requiring less computational time, algorithms that are based on such formulations require the use of specialized methods, such as constraint separation, which may not always be easy to understand, implement or use. Moreover, compact formulations are convenient to solve hard combinatorial problems appearing in subroutines of more complex algorithms (e.g., separation problems in branch-and-cut algorithms).

As pointed out by Gollowitzer et al. [START_REF] Gollowitzer | A comparison of several models for the hamiltonian p-median problem[END_REF], formulations for this problem can be viewed as composed by three sets of constraints: i) the degree constraints on the nodes; ii) one set that, together with the first, prevents more than p cycles; iii) one set that, together with the first, prevents less than p cycles. More recently, Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] contextualized this framework in terms of directed models using node variables y i that indicate whether node i plays the role of a "depot" in a cycle in the solution together with the constraint i∈V y i = p. The main argument from Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] to justify the use of these variables in the formulations is that in solutions satisfying only the degree constraints and having more than p cycles, at least one of the cycles will have no depots assigned to it. Thus, the cycle(s) with no depots assigned to can be "eliminated" by inequalities similar to the standard sub-tour elimination constraints known from the TSP (or ATSP). Alternatively, in solutions satisfying only the degree constraints and having less than p cycles, at least one cycle will have more than one depot assigned to it. Thus, the cycle(s) with more than one depot assigned to can be "eliminated" by inequalities similar to the "path-elimination" inequalities known from models of variants with multi-depots.

With respect to inequalities of type ii), the literature on compact formulations for the ATSP is vast (see, for instance, Roberti and Toth [START_REF] Roberti | Models and algorithms for the asymmetric traveling salesman problem: an experimental comparison[END_REF] and Öncan et al. [START_REF] Öncan | A comparative analysis of several asymmetric traveling salesman problem formulations[END_REF]) and thus, there are plenty of possibilities. However, one observation about this adaptation is needed, namely that the adapted constraints will need to include the y i variables. In this paper, we study two classes of models for modelling the generic constraints ii). We use an adaptation of the Desrochers and Laporte (DL) model [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF] and an adaptation of the well known multicommodity flow (MCF) model by Claus [START_REF] Claus | A new formulation for the travelling salesman problem[END_REF]. The first model is known to have a weak linear programming (LP) relaxation bound but has a rather good performance to obtain the optimal integer solution when used within a ILP package. The second model has a LP relaxation bound that is considerably stronger than the one given by the DL model. However, due to the large number of variables and constraints, it is less efficient to solve. The reader is referred to Roberti and Toth [START_REF] Roberti | Models and algorithms for the asymmetric traveling salesman problem: an experimental comparison[END_REF] for the evaluation of the performance of these two models, and other models, to obtain the optimal integer solution when used within a ILP package, for the case of the ATSP. The reason these two models are included in this study is to have an idea of how the LP relaxation bounds for the whole model change when formulations for modelling part ii) change. Although we do not expect the MCF based model to be competitive, a (substantial) improvement of the reported LP relaxation bounds might suggest equivalent and alternative solution methods that might later be explored.

The main part of the paper focuses on compact formulations for part iii). The main reason is that inequalities which prevent solutions with less than p cycles are less well known than inequalities which prevent solutions with more than p cycles. Every formulation presented in this work for part iii) stems from a common motivation, that is, the formulations contain variables that assign labels to nodes, and prevent less than p cycles by stating that different depots must have different labels and that nodes in the same cycle must have the same label. The formulations described and proposed in this study will be described in Section 3.

Due to the focus on inequalities for modelling part iii), we will also address a variant of the HpMP. This variant is motivated by some empirical observations taken from the literature, namely that the instances which are more difficult to solve are obtained with large values of p and that, in many cases, inequalities from part ii) are not needed in the models to solve instances with large values of p. Thus, to focus our study on inequalities of part iii), we also address the variant where the number of cycles is restricted to be greater than or equal to p. We denote this variant by HpMP ≥ and observe that it is modelled by considering only the two subsets of constraints i) and iii).

Our study will start by addressing models for the HpMP ≥ , which will be differentiated by the constraints used for part iii). For the HpMP, the models will be obtained by considering the models for the HpMP ≥ either augmented with the adapted DL constraints, or the adapted MCF constraints, for part ii). The formulations will be compared in terms of LP relaxation bound values (both theoretically and computationally) and on the basis of computational time requirements for attaining optimal solutions for a number of benchmark instances. One of the main conclusions of this study is that, in the context of compact formulations, it is worth looking at the more compact models with the node variables. Despite the weaker LP relaxation bounds, the fewer variables and constraints lead to faster computational times, especially when solving instances with more than 50 nodes. This might result from the use of current ILP packages that benefit strongly from built-in enhancements. redcheck again One other conclusion is that, although not outperforming more sophisticated methods known from the literature, the results provided by our compact models are often comparable, especially if p is not very close to 1 or to n 3 . Section 2 describes a base model that will be used in most of the formulations studied in this paper. In Section 3, we address formulations for the HpMP ≥ , or alternatively, formulations for preventing less than p cycles. We also make a brief discussion on the comparison of the LP relaxations of the models of the different classes discussed before. In Section 4, we consider models for the HpMP that are obtained from models described in Section 3, by augmenting them either with the adapted DL constraints or the adapted MCF constraints, for part ii). We also compare these two approaches. In Section 5, we describe the computational experiments, including a description of the instances tested and an evaluation of the results obtained with the models discussed in the previous sections. A final section concludes the paper.

Formulations -The BASE model

In this paper we will consider several compact formulations for the HpMP and HpMP ≥ . They are expressed in extended spaces involving several sets of variables in addition to the edge variables describing the HpMP in the "natural" space. The integer hull of such a formulation is the polytope coinciding with the convex hull of its solutions. Projecting the integer hull of a formulation for the HpMP onto the subspace of edge variables we obtain the convex hull of the incidence vectors of the HpMP solutions. For ease of notation, from now on, we say that an inequality is valid for the HpMP if it is valid for one of the integer hulls of the formulations presented below. To simplify the notation in the pages that follow, we also denote by δ(i) = {{i, j} ∈ E, ∀j ∈ V } the set of all edges adjacent to node i. Finally, for a given formulation F, we refer to its LP relaxation as F L , and to its optimal value as v(F).

In this section, we present a base formulation which is contained in most of the formulations for the HpMP that are presented in the pages that follow. This base formulation includes two sets of binary variables: i) variables u ij will take value 1 if the edge {i, j} ∈ E is included in one of the p cycles and 0 otherwise, and ii) variables y j will indicate whether node j ∈ V is a depot (y j = 1) or not (y j = 0) of the cycle node j is in. Consider, also, the following set of constraints: Min.

{i,j}∈E

d ij u ij (1a) s.t.
{i,j}∈δ(i)

u ij = 2, ∀i ∈ V (1b) j∈V y j = p, (1c)
u ij ∈ {0, 1}, ∀{i, j} ∈ E (1d) y j ∈ {0, 1}, ∀j ∈ V (1e)
The objective function corresponds to minimizing the sum of the weights of the edges which define the p cycles. Equalities (1b) are the usual assignment constraints stating that each node is included in one and only one cycle. Constraints (1d) define the u ij variables as binary. A solution to the formulation described by the u ij variables alone is composed of several disjoint cycles (with at least three nodes each) covering all nodes of the graph. Equality (1c) indicates the number of nodes that play the role of depots of the cycles. Finally, constraints (1e) define the y j variables as binary. We denote by "BASE" the formulation described above. Observe that for the moment there are no inequalities relating the two different sets of variables. In the following sections we describe several sets of such constraints, including sets of constraints that prevent solutions from having less than p cycles and sets of constraints that prevent solutions from having more than p cycles. Besides the two sets of variables previously defined, in some of the models to be presented later on, other sets of variables will be used to define the new sets of constraints. Also, in one class of formulations the additional variables will make unnecessary the use of these depot variables.

There are many equivalent solutions on the u ij variables that differ on the choice of the y i variable to denote the depot of each cycle. In order to reduce the number of equivalent solutions, most of the formulations we study in this work will include inequalities that enforce the following symmetry breaking (SB) strategy: a node can be a depot only if it is the node with the lowest index in the cycle it belongs to. In fact, for some of the formulations presented in this work, it is far from clear how to write valid formulations with the same sets of variables without imposing this property. The SB strategy allows us to tighten the redundant inequality u ij ≤ 1 into:

u ij + y j ≤ 1, ∀{i, j} ∈ E : i < j (2)
Proposition 1. Constraints (2) are valid for the HpMP ≥ .

Proof. Consider a pair of nodes i, j such that i < j. We show that we cannot have both y j = 1 and u ij = 1. If y j = 1, then, j must be the node with the lowest index in the cycle it belongs to. Therefore, since i < j, nodes i and j cannot be in the same cycle, and thus, these two nodes cannot be adjacent (that is,

u ij = 0).
The previous inequality will be added to most of the formulations presented later on. The computational results will also show that the addition of these inequalities is relevant for the computational times.

Preventing less than p cycles

As mentioned in the introduction, all the formulations presented in this paper for preventing solutions with less than p cycles are based on a similar motivation: i) assigning different labels to the depots and ii) guaranteeing that two consecutive nodes in a cycle have the same label. We will use the designation "continuity constraints" for constraints guaranteeing this last condition. These formulations can be viewed in two different ways. The first one depends on the type of variables used to represent the label of the node, either by using integer variables to represent the label of a node (aggregated formulations) or binary variables that assign each node to a given label (disaggregated formulations). The second view depends on the meaning of the label. It can either be the index of the node serving as a depot for a cycle (node-depot assignment) or the index, between 1 and p, of the cycle (node-cycle assignment). The aggregated models are new. The disaggregated formulations are not, although in all of them new enhancements have been included in the models. Since we are modelling the HpMP ≥ , we may have feasible solutions with more than p cycles. In such situations, we can have nodes from different cycles with the same labels. These situations will be detailed in the proofs of the validity of the models introduced in the next sections.

Aggregated formulations

In this subsection, we discuss two classes of aggregated models that include integer variables indicating the value of the label of a node and constraints guaranteeing that nodes in the same cycle must have the same label.

We start by presenting a relaxed formulation, denoted by Rel, that can be viewed as a base to build the models presented in this subsection. The formulation Rel is built from the BASE formulation by adding the generic label integer variables z i , ∀i ∈ V and the following constraints:

z i ≤ z j + (i -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (3a) z j ≤ z i + (j -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (3b) 1 ≤ z i ≤ i and integer, ∀i ∈ V (3c)
Constraints (3a) and (3b) guarantee z i = z j for any pair of adjacent nodes i and j. The Rel formulation relaxes the condition i) stated at the introduction of Section 3, namely that depots must have different labels. Thus, it allows infeasible solutions -for instance, for any value of p, any solution featuring a single Hamiltonian cycle is feasible, since z i = 1, ∀i ∈ V is feasible.

The two classes differ in the interpretation of the label of a node and are obtained by considering different sets of constraints linking the label variables z i with the depot variables y i .

Node-depot assignment formulations

In this section we consider the label of a node as the index of the node serving as depot for the cycle. More precisely, we consider the integer variables s i ∈ N, ∀i ∈ V , which indicate the index of the depot of the cycle containing node i (more precisely, if node i belongs to a cycle with depot d, then, s i = d). Models using a similar set of variables have been proposed by Burger et al. [START_REF] Burger | A node current-based 2-index formulation for the fixeddestination multi-depot travelling salesman problem[END_REF] (and later by Bektaş et al. [START_REF] Bektaş | Compact formulations for multi-depot routing problems: Theoretical and computational comparisons[END_REF]) for a multi-depot routing problem. In this section, we show how to adapt these models for the HpMP ≥ by adding depot variables. We also introduce a new (as far as we know) enhancement of the "continuity" inequalities which implies some of the constraints of the original model leading to a valid model with fewer constraints. Consider the following set of constraints:

s i ≤ s j + (i -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (4a)
s j ≤ s i + (j -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (4b) s i ≥ (i -1)y i + 1, ∀i ∈ V (4c) 1 ≤ s i ≤ i and integer, ∀i ∈ V (4d)
Constraints (4a), (4b) and (4d) are constraints (3a), (3b) and (3c) after setting s i = z i , ∀i ∈ V . Constraints (4c) and (4d) guarantee that s i = i if node i is a depot (this follows from considering y i = 1 in constraint (4c) for node i). When y i = 0, these two constraints guarantee that 1 ≤ s i ≤ i. This is consistent with the SB strategy mentioned earlier, as it implies only the node with the lowest index in a cycle can be its depot. Observe that although the label variables s i are interpreted as integer, they can be defined as continuous. This also applies to all the formulations described in this subsection.

We denote by ANDA (for "Aggregated Node to Depot Assignment") the formulation which results from augmenting BASE with constraints (4a)-(4d). Observe that since this formulation is an augmentation of Rel, we have that v(ANDA L ) ≥ v(Rel L ). Additionally, our computational results show that there exist instances for which v(ANDA L ) > v(Rel L ).

We observe that in the two cited works on the multi-depot routing problem, more elaborate upper bounding constraints on the s i variables were also necessary for the validity of the model. In the case of the HpMP ≥ , similar constraints, such as s i ≤ i -1 + y i , ∀i ∈ V , could have been added to the model. However, as shown in the next result, these constraints are not needed to obtain a valid model for the HpMP ≥ (this follows from the SB strategy imposing that s i ≤ i). Also, although not redundant in the LP relaxation of the model, these constraints only contribute for very small improvements of the LP relaxation bounds. For this reason, they were not included in the ANDA model and will not be further considered in the remainder of the text.

Proposition 2. The formulation ANDA is valid for the HpMP ≥ .

Proof. Consider a feasible solution for the HpMP ≥ and consider an assignment of values, 0 or 1, to the u ij variables corresponding to this solution. It is easy to choose an assignment of values for the remaining two sets of variables that satisfy (4a)-(4d). To see this, consider the cycles in the feasible solution to be sorted by ascending order of the indices of their nodes with lowest indices. Considering d k the node with the lowest index in cycle k such that 1 ≤ k ≤ p, set y d k = 1. For every node i, set s i = d k if node i is in some cycle k such that 1 ≤ k ≤ p, and set s i = 1 if node i is in some cycle k such that k > p. The resulting solution is feasible for the ANDA model.

To see the converse (that is, that any partition of the graph in less than p cycles violates some of these constraints), observe that a feasible solution for BASE is composed of k disjoint cycles that cover all nodes of the graph. Constraints (4a) and (4b) imply that s i j = s i l , ∀j, l ∈ {1, . . . , m} for the nodes in any cycle C = {i 1 , . . . , i m } of this solution.

Assume that k < p. In such a solution, there will be at least one cycle with at least two depots (this follows from (1c)), e.g., nodes a and b such that a < b. As mentioned before, s a = s b must hold. However, since nodes a and b are depots, we also have that y a = y b = 1. Constraints (4c) and (4d) imply s b = b > a = s a , leading to s b > s a . But this contradicts what was observed before, namely that s a = s b . Thus, any solution with less than p cycles is not feasible for the ANDA model.

The validity of constraints (2) raises the question of knowing whether the term (1 -u ij ) in constraints (4a) and/or constraints (4b) can be replaced by (1 -u ij -y j ). In fact, while it is not apparent how to enhance constraints (4b), we can provide the following stronger inequality in the case of (4a):

s i ≤ s j + (i -1)(1 -u ij ) -(j -1)y j , ∀{i, j} ∈ E : i < j (4e)
The next result proves the validity of constraints (4e) and that they are a strengthened version of (4a). The result also states that in the presence of constraints (4d), constraints (4c) are a particular case of (4e): Proposition 3. Constraints (4e) are valid for the HpMP ≥ and:

• Constraints y j ≥ 0, ∀j ∈ V and (4e) imply (4a);

• Constraints (4d) and (4e) imply (4c).

Proof. To prove the validity of constraints (4e), consider some edge {i, j} ∈ E. We assume y j = 1 (otherwise the inequality would correspond to the original valid inequality). The SB strategy implies u ij = 0 if y j = 1 and i < j, and constraints (4d) and (4c) imply s j = j if y j = 1. Therefore, the corresponding constraint (4e) becomes s i ≤ i, which is implied by constraints (4d). Thus, constraints (4e) are valid for the HpMP ≥ .

Regarding the two implications, the first is trivial. To prove the second result, consider the constraint (4e) for a given node j and node i = 1. Since s 1 = 1 (which is implied by constraints (4d)), the resulting constraint becomes (4c) for the same node j.

We denote by SANDA the formulation which results from replacing constraints (4a) and (4c) in ANDA with constraints (4e). Proposition 3 leads to: Proposition 4. The formulation SANDA is valid for the HpMP ≥ and v(SANDA L ) ≥ v(ANDA L ).

Finally, we also study the effect of adding constraints (2) to the previously defined models. We denote by ANDA + and SANDA + the models which result from augmenting the models ANDA and SANDA, respectively, with these constraints. The following proposition is an obvious consequence of the definition of these models:

Proposition 5. • v(ANDA + L ) ≥ v(ANDA L ); • v(SANDA + L ) ≥ v(SANDA L ).
The computational results reported in Section 5 indicate that the addition of constraints (2) also seems to have a greater impact on the LP relaxation bounds and on the computational times than the addition of constraints (4e). The main conclusion is that the most competitive model from this class of models is the SANDA + model.

Node-cycle assignment formulations

In this section, we consider the label of a node as the index of the cycle that node is in. More precisely, we consider the integer variables v i ∈ N, ∀i ∈ V such that, if node i belongs to the k-th cycle, then,

v i = k, ∀k = 1, . . . , p.
As noted before, to the best of our knowledge, formulations using such variables are new. Consider the following set of constraints:

v i ≤ v j + (i -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (5a) v j ≤ v i + (j -1)(1 -u ij ), ∀{i, j} ∈ E : i < j (5b) v i ≥ i j=1 y j -(i -2)(1 -y i ), ∀i ∈ V (5c) v i ≤ i j=1 y j , ∀i ∈ V (5d) 1 ≤ v i ≤ i and integer, ∀i ∈ V (5e)
Constraints (5a), (5b) and (5e) are constraints (3a), (3b) and (3c) after setting v i = z i , ∀i ∈ V . Constraints (5c) and (5d) relate the v i variables with the y i variables. When y i = 0, constraints (5c) are redundant. When y i = 1, these constraints, together with (5d), guarantee that v i = i j=1 y j , that is, the cycle number of a given node i equals the number of depot nodes with an index less than or equal to i. Constraints (5e) set lower and upper bounds on the v i variables. Observe that constraints (1c) and (5d) also imply that v i ≤ p, ∀i ∈ V , which is consistent with the definition of the v i variables given before. The equality v i = i j=1 y j when y i = 1 is consistent with the assumption that the binary node variables y i are interpreted as indicating whether node i is the depot of a cycle and is the node of smallest index in the cycle. Observe also that constraints (5c) and (5d) enforce an additional SB strategy, namely that the cycles must be sorted by ascending order of the indices of their depots. Finally, although the new variables v i are defined as integer, they can be defined as continuous without altering the validity of the model. This applies to all the formulations described in this subsection.

We denote by WANCA (for "Weak Aggregated Node to Cycle Assignment") the formulation which results from augmenting BASE with constraints (5a)-(5e). Proposition 6. The formulation WANCA is valid for the HpMP ≥ .

Proof. Consider a feasible solution for the HpMP ≥ and consider an assignment of values, 0 or 1, to the u ij variables corresponding to this solution. It is easy to choose an assignment of values for the remaining two sets of variables that satisfy (5a)-(5e). To see this, consider the cycles in the feasible solution to be sorted by ascending order of the indices of their nodes with lowest indices. Let d k be the node with the lowest index in cycle k, 1 ≤ k ≤ p, and set y d k = 1. For every node i, set v i = k if node i is in some cycle k, 1 ≤ k ≤ p. In the case that the solution has more than p cycles, set, for instance, v i = 1 if node i is in some cycle k such that k > p. The resulting solution is feasible for the WANCA model.

To see the converse (that is, that any partition of the graph in less than p cycles violates some of these constraints), observe that a feasible solution for BASE is composed of k disjoint cycles that cover all nodes of the graph. We begin by observing that for any pair of adjacent nodes i and j, constraints (5a) and (5b) imply v i = v j , which, in turn, implies that

v i k = v i l , ∀k, l ∈ {1, . . . , m}
for the nodes in any cycle C = {i 1 , . . . , i m } of this solution.

Assume that k < p. In such a solution, there will be at least one cycle with at least two depots (this follows from (1c)), e.g., nodes a and b such that a < b. As mentioned before, v a = v b must hold. However, since nodes a and b are depots, we also have that y a = y b = 1. Constraints (5c) and (5d) imply

v b = b j=1 y j ≥ 1 + a j=1 y j = 1 + v a leading to v b ≥ v a + 1.
But this contradicts what was observed before, namely that v a = v b . Thus, some of constraints (5a)-(5e) are violated by any solution featuring less than p cycles.

In contrast to what happens with the ANDA model, the LP relaxation bound of the WANCA model does not improve on the LP relaxation bound of the original BASE model.

Proposition 7. v(BASE L ) = v(WANCA L ).
Proof. As WANCA is an augmentation of BASE, v(WANCA L ) ≥ v(BASE L ) clearly holds. Therefore, it is only necessary to prove v(WANCA L ) ≤ v(BASE L ). This can be proven by showing that for any feasible solution for BASE L , there is also a feasible solution for WANCA L with the same cost. Consider a feasible solution (u b , y b ) for BASE L . Observe that for the remainder of the proof, the value of the variables y b are irrelevant. Consider now a solution (u w , y w , v w ) for WANCA L such that u w = u b , y 1 = 1, y 2 = 0, y i = 0.5, ∀i ∈ V : 3 ≤ i ≤ 2p, y i = 0, ∀i ∈ V : i > 2p and v i = 1, ∀i ∈ V . It is easy to see that this solution is feasible for WANCA L and has the same cost as the original solution.

This result shows that the linking constraints (5c) and (5d) are effective only on the integer solutions. A more detailed comparison between the two classes of aggregated models can be found in Section 3.1.3.

As mentioned earlier, constraints (1c) and (5d) imply v i ≤ p, ∀i ∈ V . These upper bounds indicate that the constraints (5a), (5b), (5c) and (5e) can be "enhanced", as follows:

v i ≤ v j + min{p -1, i -1}(1 -u ij ), ∀{i, j} ∈ E : i < j (6a) v j ≤ v i + min{p -1, j -1}(1 -u ij ), ∀{i, j} ∈ E : i < j (6b) v i ≥ i j=1 y j -min{p -1, i -2}(1 -y i ), ∀i ∈ V (6c) 1 ≤ v i ≤ min{i, p} and integer, ∀i ∈ V (6d)
Constraints (6a), (6b), (6c) and (6d) are enhanced versions of (and can replace) constraints (5a), (5b), (5c) and (5e), respectively. We denote by ANCA (for "Aggregated Node to Cycle Assignment") the formulation which results from augmenting the BASE formulation with constraints (5d) and (6a)-(6d), and omit the proof of the validity of ANCA as it is very similar to that of the validity of WANCA. Proposition 8. The formulation ANCA is valid for the HpMP ≥ .

The enhancement obtained by using the "min" coefficients in every constraint as well as by restricting the value of the variables would suggest that the ANCA model would lead to promising results with respect to LP relaxation bounds, at least when compared with WANCA. Quite surprisingly, there are no improvements for the benchmark instances considered in the computational experiment. On the other hand, there are improvements for other randomly generated instances. This, combined with the fact that ANCA is a strengthening of WANCA, proves that v(ANCA L ) ≥ v(WANCA L ) and that there exist instances for which the inequality is strict.

The validity of constraints (2) also leads to a similar question made in the context of the ANDA model, namely of knowing whether constraints (6a) and/or (6b) can be enhanced. And once again, while it is not apparent how to enhance constraints (6b), we can present an enhanced version of constraints (6a):

v i ≤ v j + min{p -1, i -1}(1 -u ij ) -min{p, i}y j , ∀{i, j} ∈ E : i < j (6e)
Proposition 9. Constraints (6e) are valid for the HpMP ≥ . Constraints y j ≥ 0, ∀j ∈ V and (6e) imply (6a).

Proof. To prove the validity of constraints (6e), consider a pair of nodes i, j such that i < j and we assume that y j = 1 (otherwise the inequality corresponds to the original valid inequality). It has already been shown (in the proof of Proposition 1) that y j = 1 implies u ij = 0 if i < j. Therefore, the corresponding constraint (6e) becomes v i ≤ v j -1. But this inequality is implied by constraints (6c) and (5d) which state that v j = j l=1 y l ≥ i l=1 y l + 1 ≥ v i + 1 when y j = 1, proving these constraints are valid. The proof of the implication is trivial.

We observe, however, that constraints (6a) cannot be "enhanced" by using the coefficient of variable y j from SANDA (i.e, adding the term "-min{p, j -1}y j " to the previous constraints). The reason for this is that the value of the coefficient (j -1) is too high, conflicting with constraints (6c) and/or (5d), leading to a not valid inequality. To see this, consider i < j -1 ≤ p. Then, the "inequality" becomes v i ≤ v j + (i -1)(1 -u ij ) -(j -1)y j . Consider also that j is a depot, that is, y j = 1. In this case, due to constraints (6c) and (5d), we have that v j = j k=1 y k . Finally, observe that u ij = 0 since i < j and j is a depot. Then the "inequality" reduces to v i + (j -i) ≤ j k=1 y k which is not necessarily valid when 1 < j -i. Observe, however, that when i = j -1, the inequality is valid, since we obtain the corresponding constraint (6e).

We denote by SANCA ("S" for "Strengthened") the formulation which results from replacing constraints (6a) in ANCA with constraints (6e). Proposition (9) leads to: Proposition 10. The formulation SANCA is valid for the HpMP ≥ and v(SANCA L ) ≥ v(ANCA L ).

The two models described in this subsection can also be enhanced by adding constraints (2) as was the case of the ANDA and SANDA models. We denote by ANCA + and SANCA + the models which result from augmenting ANCA and SANCA respectively with constraints (2). The following proposition is an obvious consequence of the definition of these models:

Proposition 11. • v(ANCA + L ) ≥ v(ANCA L ); • v(SANCA + L ) ≥ v(SANCA L ).
The computational results reported in Section 5 indicate that the addition of constraints ( 2) and (6e) has a positive impact on the LP relaxation bounds and on the computational times obtained with these models, and unlike what happened with the aggregated node-depot models, the impact of augmenting the ANCA model with constraints (6e) is comparable to that of augmenting the same model with constraints [START_REF] Bektaş | New path elimination constraints for multi-depot routing problems[END_REF]. As a result of this, the main conclusion to be drawn from the computational testing of these models is that formulation SANCA + is the most competitive aggregated node-cycle assignment model. However, regarding computational times, the computational results indicate that the aggregated node-depot formulations are more competitive and we try to provide an explanation for this behaviour in the next paragraphs.

To conclude this section on aggregated node-cycle assignment formulations, we observe that we can obtain a more compact and still valid model (although with a weaker LP relaxation) by removing constraints (6c) and (5d) from the SANCA model. Although more compact, the fact that it has a weaker LP relaxation may explain why, according to some preliminary computational testing, this model is not as competitive as the SANCA model for solving instances of the HpMP ≥ . Thus, we omit details and the proof of its validity in this paper.

Comparing the Aggregated Node-Depot and Node-Cycle Assignment formulations

In this section we emphasize the differences between the aggregated node-depot and node-cycle assignment formulations, which can explain the computational results reported in Section 5. We begin the comparison with models ANDA and WANCA. The ANDA and WANCA models differ only in the constraints linking the label variables with the depot variables. As pointed before, the set of linking constraints used in the node-cycle model does not lead to an improvement with respect to the LP relaxation of the BASE model.

With respect to comparing the ANCA model with the ANDA model, although ANDA sometimes provides slightly better LP relaxation bounds, the differences are minimal. In fact, we did not find any instance where ANCA has a strictly better LP relaxation bound than ANDA, although, in our computational experiments, we have instances for which the model ANDA augmented with constraints (5d) and (6a)-(6d) of the ANCA model resulted in LP relaxation bounds which were strictly better than those obtained with the original ANDA model. This shows that the LP relaxation of the ANDA model does not dominate that of the ANCA model. Also, for obtaining the optimal solutions, the CPU times produced by ANDA are, in general, better than the CPU times produced by ANCA, as reported in Section 5. To explain this behaviour, we recall one property of each class of models, namely that for the ANDA model we have that s i = i if node i is a depot, and with respect to the ANCA (and also the WANCA) model, we have that v i = i j=1 y j if node i is a depot. One immediate consequence of this is that in the node-depot models, the label of a node depends only on the choice of the depot of the cycle while, in node-cycle models, the label depends on the number of depots previously chosen. This might explain the different behaviour of the two classes of models for solving instances of the problem. In particular, in a branch-and-bound algorithm, when a y i variable is fixed to 1, the value of the corresponding v i variable in the node-cycle models might still remain dependent on other y j variables for which the values have not been fixed, whereas the value of the corresponding s i variable in the node-depot models always becomes a constant (since y i = 1 implies s i = i). In other words, the benefits of using constraints (4c) seem to outweigh the negatives of not being able to lower the coefficients on constraints (3a) and (3b).

With respect to comparing the model SANDA with the model SANCA, we begin by observing that one of the main differences between these models is in the lifted constraints (4e) and (6e). Consider, again, the "generic" label variables z i , ∀i ∈ V introduced at the beginning of this section and the relaxed lifted constraint:

z i ≤ z j + (i -1)(1 -u ij ) -iy j , ∀{i, j} ∈ E : i < j (7)
The lifted inequality from the SANDA model is obtained from (7) by setting s i = z i , ∀i ∈ V and adding "-(j -i -1)y j " to the right-hand side of constraints [START_REF] Burger | A node current-based 2-index formulation for the fixeddestination multi-depot travelling salesman problem[END_REF], whereas the lifted inequality from the SANCA model can be obtained from ( 7) by setting v i = z i , ∀i ∈ V and lowering the coefficients in constraints [START_REF] Burger | A node current-based 2-index formulation for the fixeddestination multi-depot travelling salesman problem[END_REF] to obtain (6e). This argument, where we start with a relaxed inequality and provide different enhancements to obtain two different valid inequalities, suggests that the strengthened ND model and the strengthened NC models have non-related LP relaxations (as confirmed by the computational results). The results also indicate that, empirically, we gain more by adding the term "-(j -i)y j " to the right-hand side of ( 7) and not reducing the coefficients (and thus, being allowed to keep the original interpretation of the s i variables) than by doing the coefficient reduction.

Disaggregated formulations

In this section we discuss two classes of disaggregated models that include binary node variables that assign nodes to cycles. The two classes differ in the interpretation of the node assignment variables. The models discussed in the two next subsections are built from models known from the literature. We emphasize, however, that in both cases a reduced and valid version of the main model are proposed for the first time. The reduced models contain half of the "continuity" constraints of the main models, and while the removal of these constraints results in slightly lower LP relaxation bounds, the lower number of constraints means these models result in lower computational times when compared with their "complete" counterparts. In spite of this, however, these reduced models are still not competitive with the aggregated models presented in Section 3.1. In one class of the models, a new set of inequalities is also proposed for the first time. The new inequalities lead to substantial improvements on the LP relaxation bounds as well as substantial reductions on the CPU times to obtain the optimal solution.

Node-depot assignment formulations

In this section we return to the interpretation of Section 3.1.1 and consider the label of a node as the index of the node serving as depot for the cycle. More precisely, we consider the binary node-depot assignment variables

s i d , ∀i, d ∈ V : d ≤ i, indicating if node d is the depot of the cycle node i is in.
Observe that these variables are defined only for d ≤ i due to the SB strategy, and that variables s d d indicate if node d is a depot. The models discussed in this subsection is built upon models known from the literature (see Gollowitzer et al. [START_REF] Gollowitzer | A comparison of several models for the hamiltonian p-median problem[END_REF], Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF] and Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF], this last reference describes a model for the directed case). As with the ANDA models, for solutions with more than p cycles, we will have nodes from more than one cycle assigned to the same depot d since the constraints of the model will limit the number of depots to p. This will be clarified in the proof of Proposition 12. The following system prevents less than p cycles:

d∈V :d≤i s i d = 1, ∀i ∈ V (8a) s i d ≤ s d d ∀i, d ∈ V : d ≤ i (8b) s i d ≤ s j d + 1 -u ij , ∀{i, j} ∈ E : i < j, ∀d ∈ V : d ≤ i (8c) s j d ≤ s i d + 1 -u ij , ∀{i, j} ∈ E : i < j, ∀d ∈ V : d ≤ i (8d)
s i i = y i , ∀i ∈ V (8e) s i d ∈ {0, 1} ∀i, d ∈ V : d ≤ i (8f)
Constraints (8a) indicate that any node i must be assigned to exactly one depot and (8b) indicates that if node j is in the cycle with node d as the depot, then, node d must be a depot. Constraints (8c) and (8d) are similar to the continuity inequalities presented in the other models and indicate that two adjacent nodes must be assigned to the same depot. Constraints (8e) relate the node-depot assignment variables with the depot variables defined in the BASE model. Constraints (8f) define the new variables as binaryhowever, it is easy to see that they can be defined as continuous without altering the validity of the model.

Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF] present an exponentially sized generalization of the continuity constraints (8c) and (8d) which consider subsets of depots, instead of single node sets, and according to the results in Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], these constraints appear to be effective in solving instances of the problem. However, since the focus of this paper is on compact formulations, we do not elaborate further on these generalized constraints.

We denote by DNDA (for "Disaggregated Node-Depot Assignment") the formulation which results from augmenting BASE with constraints (8a)-(8f). Although the proof of the validity of the DNDA model for the HpMP ≥ can be taken and adapted from the literature, we present a proof here, since it allows an easier proof of the reduced model. Proposition 12. The formulation DNDA is valid for the HpMP ≥ .

Proof. Consider a feasible solution for the HpMP ≥ and consider an assignment of values, 0 or 1, to the u ij variables corresponding to this solution. It is easy to choose an assignment of values for the remaining two sets of variables that satisfy (8a)-(8f). To see this, consider the cycles in the feasible solution to be sorted by ascending order of the indices of their nodes with lowest indices. Considering d k the node with the lowest index in cycle k such that 1

≤ k ≤ p, set s d k d k = 1. For every node i, set s i d k = 1 if node i is in some cycle k such that 1 ≤ k ≤ p,
and for all the nodes i in some cycle k such that k > p we set, for instance, s i 1 = 1. The resulting solution is feasible for the DNDA model. To see that partitions of the graph in less than p cycles are not feasible, we consider a solution with less than p cycles, and show it is not feasible for the aforementioned model. We begin by observing that for any pair of adjacent nodes i and j and any depot d, constraints (8c) and (8d) imply s i d = s j d . This, in turn, means that for any cycle with nodes {i 1 , . . . , i m }, constraints (8c) and (8d) imply the following:

s ir d = s i l d , ∀r, l ∈ {1, . . . , m}, ∀d ∈ V (8g)
To see that this model prevents solutions with less than p cycles, start by considering a solution with less than p cycles. Under this assumption, constraints (1c) and (8e) imply one of the cycles will include at least two depots -for instance, nodes a and b. This means s a a = s b b = 1. However, according to what was previously observed, as a consequence of s a a = 1, s i a = 1 for every node i in the cycle to which node a belongs must hold, including node b. However, if s b a = s b b = 1, constraint (8a) for node b is violated. Therefore, some of constraints (8a)-(8f) are violated by any solution featuring less than p cycles, and the proposition holds.

We present next two classes of valid inequalities that strengthen the LP relaxation of the DNDA model. The first class was presented in [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF] and the second is new. The following constraints state that if some node j is a depot, there must be two other nodes assigned to node j (in other words, each cycle must have at least three nodes):

j∈V :j>i s j i ≥ 2s i i , ∀i ∈ V (8h)
As mentioned in [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], these constraints are relevant for the reported computational results. The second class is an enhancement of constraints (2) and a consequence of the SB strategy:

u ij + j d=i+1 s j d ≤ 1, ∀{i, j} ∈ E : i < j (8i)
These constraints state that if some edge {i, j} belongs to one of the cycles, then, the depot of that cycle cannot be any node with an index greater than i.

Proposition 13. Constraints (8i) are valid for the HpMP ≥ .

Proof. We prove the validity of these constraints by showing that, for any pair of nodes i, j such that i < j, j d=i+1 s j d = 1 implies u ij = 0. Assume that j d=i+1 s j d = 1. In this case, the depot of the cycle to which node j belongs is some node d such that d > i. Therefore, i cannot belong to that cycle, which implies that nodes i and j cannot be adjacent, and therefore, u ij = 0 must hold.

We denote by DNDA

+ the model which results from augmenting DNDA with constraints (8h) and (8i).

The previous propositions show that:

Proposition 14. The formulation DNDA + is valid for the HpMP ≥ and v(DNDA

+ L ) ≥ v(DNDA L ).
Our computational results will show that the addition of these constraints improves the LP relaxation bound of the DNDA model and that, in several cases, the solution times become substantially smaller. However, the disadvantage of the disaggregated models, namely of the "node-depot assignment" models, is the large number of variables and constraints, which makes it difficult for solvers to solve large instances using these models. In order to potentially improve the computational times obtained with this class of models, we show next that we can define a valid model, DNDA-(DNDA-+ ), which is obtained by removing constraints (8d) from DNDA (DNDA + ). We now show DNDA-and DNDA-+ are valid formulations for the HpMP ≥ :

Proposition 15. The formulations DNDA-and DNDA-+ are valid for the HpMP ≥ .

Proof. The proof is done only for the DNDA-model since a similar proof holds for the DNDA-+ model.

We prove the validity of the DNDA-model by showing that for any cycle with nodes {i 1 , . . . , i m } in any feasible solution for this model, the equality (8g) still holds, as the rest of the proof would be the same as that of Proposition 12.

To prove that (8g) still holds, as in the proof of Proposition 12, we only need to prove the following equality holds for any pair of adjacent nodes i, j such that i < j:

s i d = s j d , ∀d ∈ V : d ≤ i, j
To prove the equality above, we consider a feasible solution for DNDA-and a pair of adjacent nodes i and j (we assume i < j). Since these nodes are adjacent, constraints (8c) imply

s i d ≤ s j d , ∀d ∈ V : d ≤ i, j.
Observe that constraints (8a) and (8f) also imply that there exists exactly one node, d , such that s i d = 1. Since i and j are adjacent, constraints (8f) and (8c) imply that s j d = 1. Constraints (8a) for nodes i and j imply a∈V :a =d s i a = a∈V :a =d s j a = 0, which, combined with 0 ≤ s l a , ∀l, a ∈ V , implies

s i a = s j a = 0, ∀a ∈ V : a = d . Therefore, s i d = s j d , ∀d ∈ V : d ≤ i, j holds.
The definition of these models leads to the following result:

Proposition 16. v(DNDA + L ) ≥ v(DNDA- + L ) and v(DNDA- + L ) ≥ v(DNDA-L ).
The results from our computational results indicate that in terms of computational times, the models DNDA-+ and DNDA + are the more competitive models, showing that the addition of constraints (8h) and (8i) is relevant, independently of using the complete or the reduced model.

Node-cycle assignment formulations

Unlike any other model presented in this study, the models discussed in this subsection do not require the y i variables to guarantee the requirement of the p cycles. Thus, they might be viewed as not following the paradigm of the three previous classes of models. The "complete" model discussed in this section is from [START_REF] Gollowitzer | New models for and numerical tests of the hamiltonian p-median problem[END_REF] (see, also, [START_REF] Gollowitzer | A comparison of several models for the hamiltonian p-median problem[END_REF]). However, the variant of the model using only half of the so-called "continuity" constraints is new. The model includes the binary variables v k i , ∀i ∈ V, ∀k = 1, . . . , p : k ≤ i, indicating whether node i is in cycle k. Observe that these variables are only defined for k ≤ i, since, in order to reduce the number of equivalent solutions, we follow Gollowitzer et al. [START_REF] Gollowitzer | New models for and numerical tests of the hamiltonian p-median problem[END_REF] and consider that a node i cannot be in a cycle k such that i < k. However, despite reducing the number of equivalent solutions, this strategy does not eliminate symmetries entirely, as, for some cycle k ∈ {1, . . . , p -1}, the node with the lowest index of cycle k can have an index greater than the index of the node with the lowest index of cycle k + 1. Similarly to what happened in the definition of the variables in the previous class of models, for solutions with more than p cycles, we will have more than one cycle assigned to the same index k. This will be clarified in the proof of Proposition 17. Consider the following model to prevent less than p cycles:

min{i,p} k=1 v k i = 1, ∀i ∈ V (9a) v k i ≤ v k j + 1 -u ij , ∀{i, j} ∈ E : i < j, ∀k = 1, . . . , p : k ≤ i (9b) v k j ≤ v k i + 1 -u ij , ∀{i, j} ∈ E : i < j, ∀k = 1, . . . , p : k ≤ i (9c) i∈V :i≥k v k i ≥ 3, ∀k = 1, . . . , p (9d)
v k i ∈ {0, 1}, ∀i ∈ V, ∀k = 1, . . . , p (9e)
Constraints (9a) guarantee that each node is assigned to exactly one cycle. Constraints (9b) and (9c) guarantee that two adjacent nodes are in the same cycle. Constraints (9d) guarantee that each cycle has at least three nodes. Observe that a right-hand side value of 1 would suffice to obtain a valid formulation. Indeed, our computational results show that the stronger right-hand side leads to improvements on the LP bounds -however, these modifications do not drastically affect the behaviour of these formulations when compared with the other models presented in this paper. These constraints have a similar meaning to constraints (8h) presented in Section 3.2.1 for the DNDA models. Finally, constraints (9e) define the v k i variables as binary. For the DNDA models, we have observed that the s i d variables, ∀i, d ∈ V : d = i, can be defined as continuous. However, for the validity of the DNCA models, the v k i variables need to be defined as binary. As noted above, the previous model does not require the y i variables to guarantee the requirement of the p cycles. This requirement is guaranteed by the range of the variation of the index k in the variables v k i . We denote by DNCA (for "Disaggregated Node-Cycle Assignment") the formulation which results from augmenting BASE with constraints (9a)-(9e) and removing constraints (1c) and (1e) and the y i variables. The proof of the validity of the DNCA model for the HpMP ≥ can be adapted from the one given in Gollowitzer et al. [START_REF] Gollowitzer | New models for and numerical tests of the hamiltonian p-median problem[END_REF]. However, as was the case of the DNDA model, we present a proof here, since it allows an easier proof of the reduced model. Proof. Consider a feasible solution for the HpMP ≥ and consider an assignment of values, 0 or 1, to the u ij variables corresponding to this solution. It is easy to choose an assignment of values for the remaining two sets of variables that satisfy (9a)-(9e). To see this, consider the cycles in the feasible solution to be sorted by ascending order of the indices of their nodes with lowest indices. For every node i, set v k i = 1 if node i is in some cycle k such that 1 ≤ k ≤ p, and for all the nodes i in some cycle k such that k > p we set, for instance, v 1 i = 1. The resulting solution is feasible for the DNCA model. To prove the converse we begin by observing that for any pair of adjacent nodes i and j and any cycle k, constraints (9b) and (9c) imply v k i = v k j . This, in turn, means that for any cycle with nodes {i 1 , . . . , i m }, constraints (9b) and (9c) imply the following:

v k ir = v k i l , ∀r, l ∈ {1, . . . , m}, ∀k = 1, . . . , p (9f)
We consider, now, a solution, S, with fewer than p cycles, and show that it is not feasible for the DNCA model. Constraints (9d) imply that for each k = 1, ..., p there exist some nodes i such that each v k i = 1. But since the index k ranges from 1 to p and the solution S has fewer than p cycles, at least one of these cycles must have one node, say i, such that v r i = 1 for a given r, and another node, say j, such that v s i = 1 for a given s, r = s. But this contradicts (9c) since these constraints imply that s = r. Therefore, some of constraints (9a)-(9e) are violated by any solution featuring less than p cycles, and the proposition holds.

Similarly to what was done for model DNDA, we can also define a new model, DNCA-, which is obtained from DNCA by removing constraints (9c). We now show that DNCA-is a valid formulation for the HpMP ≥ : Proposition 18. The formulation DNCA-is valid for the HpMP ≥ .

Proof. We prove the validity of the DNCA-model by showing that for any cycle with nodes {i 1 , . . . , i m } in any feasible solution for this model, the equality (9f) still holds, as the rest of the proof would be the same as that of Proposition 17.

To prove that (9f) still holds, as in the proof of Proposition 17, we only need to prove the following equality holds for any pair of adjacent nodes i, j such that i < j:

v k i = v k j , ∀k = 1, . . . , p : k ≤ i, j
To prove the equality above, we consider a feasible solution for DNCA-and a pair of adjacent nodes i and j (we assume i < j). Since these nodes are adjacent, constraints (9b) imply v k i ≤ v k j , ∀k = 1, . . . , p : k ≤ i, j. Observe that constraints (9a) and (9e) also imply that there exists exactly one cycle index, k , such that v k i = 1. Since i and j are adjacent, constraints (9e) and (9b) imply that v k j = 1. Constraints (9a) for nodes i and j imply min{i,p} q=1:q =k v q i = min{j,p} q=1:q =k v q j = 0, which, combined with 0 ≤ v q l , ∀l ∈ V, q = 1, . . . , p, implies v q i = v q j = 0, ∀q = 1, . . . , p : q = k. Therefore,

v k i = v k j , ∀k = 1, . . . , p ∈ V : k ≤ i, j holds.
The definition of the DNCA and DNCA-models results in the following:

Proposition 19. v(DNCA L ) ≥ v(DNCA-L ).
The computational results indicate that, although not considerably worse, the node-cycle models are not as good as the node-depot models. In terms of computational times, these results are interesting since the node-cycle models have much fewer variables than the node-depot models. We give an explanation in the computational results section.

Comparing aggregated formulations with disaggregated formulations

In this subsection we make a brief comparison on the behaviour of the models in each of the four classes just described. As pointed out (and confirmed in the computational results section) the aggregated models have worse LP relaxation bounds than the disaggregated models, but are strongly preferred for obtaining the optimal integer solutions. In order to give some insight in the relationship between the models, observe that the variables of the node-depot (node-cycle) aggregated models can be related with the variables of the node-depot (node-cycle) disaggregated models as follows:

s i = i d=1 ds i d , ∀i ∈ V ( 10
)
v i = min{i,p} k=1 kv k i , ∀i ∈ V (11)
Observe that by using [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF] (respectively ( 11)), we can add any constraint of the aggregated nodedepot (respectively node-cycle) assignment models to the disaggregated node-depot (respectively nodecycle) models. We can, for instance, consider two augmented disaggregated models -the first is obtained from DNDA by adding constraints (4a) and (4b) (where the s i variables are replaced according to [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF]), and the second is obtained from DNCA by adding constraints (6a) and (6b) (where the v i variables are replaced according to [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF]). Some preliminary computational testing shows there are instances where these augmented disaggregated models have strictly better LP relaxation bounds than their non-augmented counterparts, showing some of the constraints found in the aggregated models are not implied by the disaggregated models. However, while, in terms of their LP relaxations, these augmented disaggregated models are always at least as good as the corresponding aggregated models, they are not competitive with their aggregated counterparts regarding computational times (due to having many more variables and constraints). A similar situation, although with better LP bounds, arises when the disaggregated models are augmented with the lifted constraints (4e) or (6e).

Finally, the relationship between the two sets of variables, aggregated and disaggregated, may also suggest new inequalities to enhance the LP relaxations of the aggregated models. For instance, in Sections 3.2.1 and 3.2.2, it was stated that the addition of constraints (8h) helps in speeding up computational times and that the enhanced right-hand side in constraints (9d) results in improvements to the LP relaxation bounds. The following inequalities can be derived from those constraints:

i∈V s i ≥ i∈V 3iy i + n -3p ( 12
) i∈V v i ≥ p k=1 3k + n -3p (13)
These constraints are a consequence of the fact that it is assumed that each cycle must include at least three nodes, and are valid for the aggregated node-depot and node-cycle assignment models respectively. The next two propositions show how to derive these inequalities from the disaggregated models and the linking constraints [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF] and [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], also indirectly proving their validity.

Proposition 20. Constraints (1c), (8a), (8e), (8h) and (10) imply constraints [START_REF] Glaab | A new variant of a vehicle routing problem: Lower and upper bounds[END_REF].

Proof. To derive constraints [START_REF] Glaab | A new variant of a vehicle routing problem: Lower and upper bounds[END_REF] from constraints of the DNDA models, we begin by observing that constraints (8h) and (8e) imply:

d∈V i∈V (d -1)s i d ≥ d∈V 3(d -1)y d
Rearranging the last term and using (1c), we get:

d∈V 3(d -1)y d = d∈V 3dy d - d∈V 3y d = d∈V 3dy d -3 d∈V y d = d∈V 3dy d -3p
This leads to:

d∈V i∈V (d -1)s i d ≥ d∈V 3dy d -3p
Observe also that constraints (8a) imply:

i∈V d∈V

s i d = d∈V i∈V s i d = n
The previous two inequalities lead to:

d∈V i∈V (d -1)s i d + d∈V i∈V s i d = i∈V d∈V ds i d ≥ d∈V 3dy d + n -3p
Finally, by using [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF], we obtain [START_REF] Glaab | A new variant of a vehicle routing problem: Lower and upper bounds[END_REF].

Proposition 21. Constraints (9a), (9d) and (11) imply constraints [START_REF] Glaab | The hamiltonian p-median problem[END_REF].

Proof. To derive constraints (13) from constraints of the DNCA models, we begin by observing constraints (9d) imply:

p k=1 i∈V (k -1)v k i ≥ p k=1 3(k -1) = p k=1 3k -3p
Observe also that constraints (9a) imply:

i∈V p k=1 v k i = p k=1 i∈V v k i = n
The previous two inequalities lead to:

p k=1 i∈V (k -1)v k i + p k=1 i∈V v k i = p k=1 i∈V kv k i ≥ p k=1 3k + n -3p
Finally, by using [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], we obtain [START_REF] Glaab | The hamiltonian p-median problem[END_REF].

However, while constraints (8h) and the enhanced right-hand side in constraints (9d) appear to have a positive impact on the computational times or on the LP relaxation bounds, some preliminary computational results suggest that constraints ( 12) and ( 13) result in minor improvements to the LP relaxation bounds ((13) results in slightly greater improvements than [START_REF] Glaab | A new variant of a vehicle routing problem: Lower and upper bounds[END_REF]) and no considerable improvements to the computational times when added to the aggregated node-depot and node-cycle assignment models.

On the other hand, these negative results do not prevent further study in the relation between aggregated and disaggregated models.

Preventing more than p cycles

In this section we address the HpMP and present constraints to prevent solutions with more than p cycles. It is far from clear how to write compact sets of such constraints without "directing" the graph, that is, for edge {i, j} we also distinguish whether the edge is used in the direction from i to j (arc (i, j) is used) or is used the direction from j to i (arc (j, i) is used). Thus, we introduce the binary variables x ij ∈ {0, 1}, ∀i, j ∈ V , which indicate whether arc (i, j) is in the solution, and add, to the BASE model, the sets of constraints:

u ij = x ij + x ji , ∀{i, j} ∈ E (14a) j∈V :j =i x ij = 1, ∀i ∈ V (14b) j∈V :j =i x ji = 1, ∀i ∈ V (14c) x ij ∈ {0, 1}, ∀i, j ∈ V : {i, j} ∈ E (14d)
Constraints (14b) and (14c) are usually included in MILP formulations for the ATSP. Observe that by using (14a) we can remove the u ij variables from the models and obtain "pure" directed models. In fact, in the computational testing for the HpMP, pure directed models where used. Observe also that the constraints x ij + x ji ≤ 1, ∀{i, j} ∈ E (that result from (14a)) were not added to the directed models since they are implied by the continuity constraints of the models.

We present, in the next two subsections, adaptations of two models known from the Asymmetric TSP to prevent solutions with more than p cycles. As far as we know, these adaptations are new and, as noted in the introduction, we describe two models with substantially different LP strengths in order to evaluate the contribution of this subsystem to the model for the HpMP.

The adapted Desrochers and Laporte formulation

Besides the y j and x ij variables already introduced for the directed BASE model, the formulation presented in this subsection also uses non-negative variables u i , ∀i ∈ V , that can be interpreted as indicating the position of a node in the cycle it is in (assuming an orientation in the cycle). Consider the following inequalities:

u j ≥ u i + 1 -M (1 -x ij + y j ) + (M -2)x ji , ∀i, j ∈ V, i = j (15a) (1 -y j ) ≤ u j ≤ (M -1)(1 -y j ), ∀j ∈ V (15b)
In inequalities (15a) and (15b), the value M equals n -3(p -1), which is the size of the largest cycle in any feasible solution. Inequalities (15b) indicate that a node i is a depot (y i = 1) if and only if its position, in the cycle, is zero (u i = 0) (observe that as a consequence, we have u j ≥ 1 for any node j that is not a depot). The requirement that y i = 1 if and only if u i = 0 is not needed for the validity of the model, but it helps the explanation of the inequalities (15a) in the case when y i = 1.

The inequalities (15a) are an adaptation of the well-known Desrochers and Laporte (DL) inequalities presented in [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF] and guarantee that the position of two consecutive nodes in a cycle must increase by at least 1, except if node j is a depot. The adaptation includes the y j variable in one of the terms of the inequality to guarantee this last condition. For a pair of nodes i, j such that x ij = 1, when y j = y i = 0, these inequalities have the standard interpretation, that is, the combined effect of the two inequalities for pairs i, j and j, i (with i and j not being depots) states that u j = u i + 1. If one of those nodes is a depot, however, the effect of the two constraints for this pair of nodes changes. Assuming again that x ij = 1, if y i = 1, the two inequalities for pairs i, j and j, i state that 1 ≤ u j ≤ M + 1, and if y j = 1, the two inequalities state that -1 ≤ u i ≤ M -1. In any case, we see that a cycle containing a depot is not forbidden.

These constraints, added to the BASE model, eliminate solutions with more than p cycles. To see this, consider a feasible solution with more than p cycles. If there are more than p cycles, at least one of these cycles does not have a depot. Using (15a) in a circular fashion along the arcs of the cycle we obtain a contradiction.

To conclude we make the following observations: i) Although we have "interpreted" the u i variables as indicating the position of a node in the cycle, they need not be defined as integer (and are not). Also, as the model is defined now, we do not guarantee that the first non-depot node in the cycle has position equal to 1. Thus, for a cycle (j, q, r) with j being the depot, we may have u j = 0, u q = 3 and u r = 4.

With respect to the HpMP, we maintain the same designation of the models that have been used for the HpMP ≥ , noting that the models for the HpMP are obtained from the models for the HpMP ≥ by replacing the edge variables u ij with the arc variables x ij (according to equality (14a)), by replacing constraints (1b) and (1d) with constraints (14b), (14c) and (14d), and by adding the u i variables and constraints (15a) and (15b).

The adapted multi-commodity flow formulation

Instead of the adapted DL constraints, we can use flows to prevent more than p cycles. There are several ways to derive flow based models for preventing more than p cycles. In this section, we present an adaptation of the well known multicommodity flow (MCF) model by Claus [START_REF] Claus | A new formulation for the travelling salesman problem[END_REF] for the ATSP. As noted in the introduction, for the ATSP, MCF based models are known to provide substantially better LP relaxation bounds than the ones provided by the DL model. Thus, although we do not expect the MCF based model to be competitive (in regards to computational times) with the DL based model, their inclusion in this study serves to have an idea of how the LP relaxation bounds for the whole model change when models for the subsystem preventing more than p cycles change.

Besides the variables x ij and y j , consider the variables g q ij , for all arcs (i, j) ∈ A and for all nodes q ∈ V : q = j, indicating whether arc (i, j) and node q are in one of the cycles and arc (i, j) is in the path from node q to the depot of the cycle. Clearly, the variables g q ij are not defined for q = j.

Consider the following constraints:

i∈V g q qi + y q = 1, ∀q ∈ V (16a) i∈V g q ji + y j ≤ 1, ∀q, j ∈ V : q = j (16b) i∈V g q ij -
i∈V g q ji ≤ y j , ∀q, j ∈ V : q = j (16c)

g q ij ≤ x ij , ∀(i, j) ∈ A, ∀q ∈ V : q = j (16d) g q ij ∈ {0, 1} ∀(i, j) ∈ A, ∀q ∈ V : q = j (16e)
Constraints (16a) state that for any node q, either the node is the depot of a cycle or it is the origin of flow in some cycle (that is, there is a path from node q to the depot of the same cycle). Constraints (16b) state that a given node j can only be the depot of a cycle if it is not in the path (more precisely, there is no arc leaving node j) from any other node q to the depot of the cycle it belongs to. Constraints (16c) state that if there exists an arc coming into node j in the path from a node q, then either node j is the depot of a cycle or there exists an arc leaving node j in the path from node q. Constraints (16d) are the standard constraints linking the g q ij variables with the x ij variables. Finally, constraints (16e) define the g q ij variables as binary.

To see that these constraints, when added to the BASE model, eliminate solutions with more than p cycles, observe again that as pointed out before, if there are more than p cycles, at least one of these cycles does not have a depot. Consider some node q in that cycle. The corresponding constraint (16a) guarantees that node q is the origin of a flow (since it is not a depot). By constraints (16c), that flow will traverse the entire cycle without being "absorbed" by any depot, and assuming node r is the node which precedes node q in that cycle, then, i∈V g q ir = 1 = i∈V g q ri holds. If node r precedes node q, then, x rq = 1, and, by constraints (14b), x ri = 0, ∀i ∈ V \ {q}, also holds. But constraints (16d) imply g q ri ≤ x ri = 0 ∀i ∈ V \{q}, and g q rq = 1 cannot hold since the variable g q rq is not defined. We therefore arrive at the conclusion that i∈V g q ri = 0 must hold, but this contradicts what was observed, since i∈V g q ri = 1 should hold. This contradiction is a consequence of considering a cycle with no depots, and therefore, these constraints eliminate solutions with more than p cycles.

Computational Results

In this section, the numerical results from our computational experiment are presented. Section 5.1 introduces the benchmark instances as well as the hardware and software configurations used in this computational experiment. Section 5.2 includes an analysis of the computational results obtained with all the models presented in this work for the HpMP ≥ , and has two subsections. In Section 5.2.1, we analyse the results obtained with the aggregated models for the HpMP ≥ , and in Section 5.2.2, we analyse the results obtained with the disaggregated models for the HpMP ≥ . Finally, Section 5.3 includes an overview of the results obtained for the HpMP. Regarding this last section, we observe that only a few models for the HpMP were tested, and these were selected based on the results reported in Section 5.2 for the HpMP ≥ .

Hardware / Software Configurations and Test Instances

For this computational experiment, we use a subset of the set of symmetric instances for the TSP from TSPLIB [START_REF] Reinelt | TSPLIB -a traveling salesman problem library[END_REF]. This subset includes 20 instances -gr24, fri26, bayg29, swiss42, att48, gr48, hk48, eil51, berlin52, brazil58, st70, eil76, pr76, rat99, kroA100, kroB100, kroC100, kroD100, kroE100 and rd100 -with the number of nodes ranging from 24 to 100. The instances whose names are bold have their edge weight functions modified to the euclidean distance (without rounding to integer values).

All tests were run on a computer with an Intel ® Core TM i7-4790 CPU, 8GB of DDR3-1600 RAM running Windows 10 Pro, version 21H2, within which CPLEX 20.1.0 Concert Technology for C++ was used. All CPLEX parameters are set to their default values, with the exception of the time limit, which is set to one hour -if an instance is not solved within the time limit, an optimal value is not obtained. Instead, an interval is obtained, to which the optimal value is guaranteed to belong.

Results for the HpMP ≥

In this subsection we analyse the results for the HpMP ≥ taken from the four classes of models. We compare the performance of the aggregated models in Section 5.2.1 and of the disaggregated models in Section 5.2.2. In this last subsection we also compare the performance of aggregated versus disaggregated models.

For each one of the instances described before, each model for the HpMP ≥ is tested for three different values of p, namely p *

1 = n 4 , p * 2 = p * 1 +p * 3 2
and p * 3 = n 3 . For this testing, p * 1 and p * 3 are considered as they are the two largest values of p considered in the benchmark instances for the HpMP, and p * 2 is an intermediate value. Additionally, p * 3 is also of interest, as it is the only value of p for which the constraints used to prevent more than p cycles are not needed. We observe that this section (and Section 5.3) only includes the average data -the detailed results can be found in the appendix.

Results for Aggregated Models

In this subsection we compare all the aggregated models defined in Section 3.1. We compare the models within each class (node-depot and node-cycle assignment models), and also between the different classes. We also evaluate the impact of augmenting these formulations with constraints (4e), (6e) and ( 2).

Table 1 presents the average LP and RN gaps as well as the CPU times to obtain the optimal solutions of the 10 models previously described, four node-depot models and six node-cycle models, and Figure 1 includes two charts which, for a certain time threshold t, indicate how many instances were solved using the corresponding model in up to t seconds. The LP gap can be calculated as follows: LP = 100 × z * -z l z * , where LP is the LP gap, z * is the optimal value of that instance and z l is the LP relaxation bound value. The RN (for "Root Node") gap can be calculated using a similar equation, replacing z l with z r , where z r is the lower bound obtained by limiting CPLEX to the root node of the branch-and-bound tree. In most cases, z r is considerably larger than z l since CPLEX may also add generic cuts which improve this lower bound and which are not considered when calculating z l . In Table 1, each row corresponds to one aggregated formulation presented in Section 3.1. For each row, the first three columns have the average LP gaps obtained with the corresponding model for all instances such that p = p * 1 , p = p * 2 and p = p * 3 , respectively. The next three columns have the average RN gaps obtained with the corresponding model for all instances, and the last three columns indicate the average CPU times obtained with the corresponding model for all instances. Next to each average CPU time entry we include the number of tests for which either the time limit was reached or CPLEX ran out of memory.

Average LP Gaps (%)

Average With respect to the aggregated node-depot assignment formulations, the computational results reported in Table 1 and Figure 1 indicate that the models ANDA + and SANDA + lead to lower computational times when compared with the corresponding non-enhanced models (although the difference is not very Clearly, the LP relaxation bounds of ANDA + and SANDA + are also always better than those of ANDA and SANDA respectively (this is obvious considering the definition of these models, although the difference is not very substantial).

A similar observation holds when we move from ANDA to SANDA and from ANDA + to SANDA + and the main conclusion is that the most competitive model from this class of models is the SANDA + model.

The addition of constraints (2) also seems to have a greater impact on the LP relaxation bounds and on the computational times than the addition of constraints (4e) (observe that the average gaps of ANDA + are slightly lower than those of SANDA, and the same applies to the average computational times). Finally, we also observe that, for p * 3 , all the average times are very similar. This is due to the fact that, for many instances with p = p * 3 , all models reach the time limit of one hour, which reduces the differences between the average times of all the models.

With respect to the aggregated node-cycle assignment formulations, these computational results indicate that although only small improvements are obtained in the LP relaxation bounds when adding constraints (2), (as was the case of the aggregated node-depot assignment models), the results also show that the addition of these constraints appears to be effective in the resolution of the instances (once again, similarly to what happened with the aggregated node-depot models). Model ANCA seems to benefit the most with the addition of constraints [START_REF] Bektaş | New path elimination constraints for multi-depot routing problems[END_REF].

The addition of the lifted constraints (6e) only leads to very small improvements of the LP relaxation bounds. However, similarly to what happened with the addition of constraints (2), the CPU times to obtain the optimal solution appear to improve with the inclusion of these lifted constraints. These constraints also seem to result in greater improvements to the computational times obtained with the aggregated node-cycle assignment models than constraints (4e) do for the aggregated node-depot assignment models.

As previously mentioned, the linking constraints (6c) and (5d), although not needed to define a valid model, help reduce the computational times. In addition to this, as mentioned earlier, constraints (6e) and ( 2) also further improve the computational times obtained with these models.

Despite resulting in considerably better computational times, SANCA + (the aggregated node-cycle assignment model with the lowest average computational times) is still not competitive with SANDA + (the aggregated node-depot assignment model with the lowest average computational times). We have already given one explanation for this behaviour in Section 3.1.3. One additional explanation for the differences in average CPU times between the aggregated node-cycle and node-depot assignment formulations might be related with the differences in average RN gaps between one class of formulations and the other. While the LP gaps between all the aggregated models are quite similar, the aggregated node-depot models seem to have better average RN gaps than their node-cycle counterparts. In other words, CPLEX appears to be able to generate stronger cuts from the node-depot formulations. This can also explain the lower computational times obtained with the aggregated node-depot formulations. Finally, we have also provided the bounds given by the BASE model and observe that the LP relaxation bounds obtained with BASE are not often considerably worse than those obtained with the proposed aggregated models. In fact, for all the instances in the benchmark set, the LP relaxation bounds given by ANCA are not better than those obtained with the BASE model. This was already pointed out in Section 3.1.3, where we have also stated that some testing with randomly generated instances shows that there are instances for which v(ANCA L ) > v(BASE L ). The small differences between the LP gaps obtained with BASE and with the other aggregated formulations suggests that further improvements for these aggregated models may be a topic of interest for future work.

Results for Disaggregated Models

In this subsection we start by comparing the disaggregated models, node-depot with node-cycle, and then compare the aggregated models with the disaggregated ones. Table 2 follows the same format as Table 1, and Figure 2 follows the same format as Figure 1 9) 2825 ( 15) 2885 ( 16) The results from our computational testing show that the LP relaxation bounds of the DNDA-model are only slightly worse than those of the DNDA model, but this difference is much larger when comparing DNDA-+ and DNDA + . This may explain why the CPU times obtained with DNDA-are considerably better than those obtained with DNDA and why this does not happen when comparing DNDA-+ with DNDA + . As a conclusion, the models DNDA-+ and DNDA + are the most competitive models, showing that the addition of constraints (8h) and (8i) is relevant, independently of using the complete or the reduced model. The results also indicate that for the node-cycle case, the reduced model is not as advantageous as in the node-depot case. Also, and although not considerably worse, the node-cycle models are not as good as the node-depot models. In terms of computational times, these results are interesting since the nodecycle models have fewer variables and constraints than the node-depot models. One explanation for this behaviour is that in the DNCA models, the v k i variables need to be defined as binary. This is in contrast with the DNDA models, where only the s d d variables need to be defined as binary. That is, the DNDA models are still valid if we relax the definition of the s i d , ∀i, d ∈ V : d ≤ i variables to be continuous.

Additionally, unlike all the models of the other classes and as mentioned in the beginning of Section 3.2.2, the DNCA models still allow many equivalent solutions, despite the fact that some of these equivalent solutions are eliminated by only defining the v k i variables for k ≤ i. These symmetries can have a negative impact on the computational times obtained with these models.

We also observe that when solving with the DNDA model, CPLEX has run out of memory for all instances with n ≥ 99 and p = p * 2 or p = p * 3 . A similar situation happens with the DNDA-model for all instances with n ≥ 99 and p = p * 3 . This is likely to be a consequence of the larger number of constraints and variables included in these disaggregated node-depot assignment models. However, the same does not happen when solving models DNDA-+ and DNDA + , as these rarely resulted in CPLEX running out of memory, most likely due to the improvements to the LP relaxation bounds. Unlike what happened with models DNDA-and DNDA, and similarly to models DNDA-+ and DNDA + , models DNCA-and DNCA also did not result in CPLEX running out of memory. This is likely to be a consequence of the lower number of variables and constraints these models have when compared with their node-depot counterparts.

We can also compare the results of the disaggregated models with the results from the aggregated models reported in the previous subsection. As noted in the introduction, and although providing better LP relaxation bounds (the model DNDA + provides substantially better LP bounds), for obtaining the optimal solutions the disaggregated models are not competitive and require considerably more CPU time when compared with the aggregated models. This difference between the CPU times obtained with the aggregated and disaggregated models might also be explained by the differences in RN gaps, as the average RN gaps obtained with the aggregated models are very similar to (if not lower than) the average RN gaps obtained with their disaggregated counterparts. The DNDA-+ and DNDA + formulations are the only exceptions, since they have the best average RN gaps out of all the models tested in this work. Finally, comparing with the LP gaps obtained with BASE, we observe that the LP bounds produced by DNDA + and DNDA-+ are considerably better, in contrast to what happens with the other models, either aggregated or disaggregated. This suggests that a more complete study on the relationship between disaggregated and aggregated models, similarly to what has been described in Section 3.3, may lead to improvements on the aggregated formulations.

Results for the HpMP

In this subsection we analyse the results obtained for the HpMP. Based on the results for the HpMP ≥ , we only compare the results from a few of the aggregated models, namely the models SANDA + and SANCA + .

For each one of the instances described in Section 5.1, we follow the literature and we test the models for the values p 1 = n 10 , p 2 = n 7 , p 3 = n 5 , p 4 = n 4 and p 5 = n 3 . Tables 3 (which includes average LP and RN gaps) and 4 (which includes average CPU times) follow the same format as Tables 1 and2 is the largest for p 5 (and even then, it is not substantial, at only 0.13% on average). Regarding average computational times, these follow a trend similar to the LP gap averages (that is, the differences between the CPU times grow with p). These also make it clear that SANDA + is the most competitive of these two models -even though there are some instances where the computational times obtained with SANCA + are lower than those obtained with SANDA + , those differences are typically very minor (for more detailed results, we forward the reader to the appendix). Finally, an interesting observation arises when comparing the models for the HpMP with the models for the HpMP ≥ -for all the instances with p = p 4 and p = p 5 , these models produce the same LP relaxation bounds and optimal values, with the single exception of instance fri26 with p = p 4 . In other words, for the vast majority of instances such that p ≥ p 4 , the adapted DL constraints (or any other constraints used to prevent more than p cycles) are not required, as the optimal solution will include p cycles regardless of whether such constraints are included. This behaviour was pointed out in the introduction of this paper and was used as a motivation to focus our study on the HpMP ≥ . Additionally, observe that the optimal values for all instances of the HpMP with p = p 5 must coincide with the optimal values of the corresponding instances of the HpMP ≥ with p = p * 3 , as in this case, there can never be solutions with more than p 5 cycles. We remark, however, that the RN gaps are better for the HpMP ≥ but these improvements are not propagated for the overall computational times.

Although we have made a study on compact formulations for the HpMP, we observe that the average gaps reported in Table 3, for the smaller values of p (p 1 to p 3 ), are comparable with the best LP relaxation bounds known from the literature (including the bounds associated to the models with exponentially sized set of inequalities reported by Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] and Barbato and Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF]). This also justifies the focus of this paper in studying these models for large values of p. The reported average gaps are, however, much worse for the biggest values of p, notably when compared with the bounds in the recent work by Barbato and Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF], where valid inequalities especially tailored for large values of p have been developed. As noted in the introduction, despite having a good performance to obtain the optimal integer solution when used within a ILP package, the adapted DL model is known to result in weak LP relaxation bound. Thus, we have also modelled the system preventing solutions with more than p cycles with the adapted MCF system (instead of the adapted DL model) described in Section 4.2. Since the HpMP model with the adapted MCF constraints contains many more variables than the HpMP model with the adapted DL constraints, we do not expect the former to be competitive with the latter to obtain the optimal solutions. However, we have compared the LP relaxation bounds with the adapted DL and with the adapted MCF systems using SANDA to prevent less than p cycles, in order to see whether some insight can be obtained with this experiment. A few computational results show that the differences are minimal, even for the smallest value of p for which the models for the HpMP were tested, that is p = n 10 . In order to focus our analysis, we have also extended this comparison for p = 2, p = 3 and p = 4, and we have only performed these tests for instances with n ≤ 76. These results are summarized in Table 5.

Looking at Table 5, we observe that the difference between the average gaps of both models is lower for larger values of p, and while it is substantial for p = 2, it is not so much so for p = 4. Also, while Average LP Gaps (%) p = 2 p = 3 p = 4 p = n 10 DL / SANDA 2.91% 2.10% 1.74% 1.47% MCF / SANDA 1.08% 1.45% 1.46% 1.07% Table 5: Average LP gaps obtained with models for the HpMP featuring the adapted DL and MCF constraints, for each value of p and for n ≤ 76 there is some difference between the gap averages for p = n 10 , we observe that these averages are skewed by the test results obtained for the instances with less than 30 nodes, since, for such instances, n 10 = 2. Indeed, if we calculate the average gaps for p = n 10 only for instances with more than 40 nodes (and up to 76 nodes), we see that the average LP gaps are now 1.36% and 1.29% for DL and MCF respectively. Considering what was said before, the results for p = n 10 also indicate that, for the values of p associated to the usual benchmark instances of the HpMP, the choice of the model for the system to prevent solutions with more than p cycles should be based more on computational issues. The results also indicate that for small values of p, we observe gains in the LP relaxation bounds when using MCF in place of DL. These results are consistent with what is known from the ATSP and indicate that if if the study is focused in small values of p : p > 1 then alternative models should be investigated.

Comparison with branch-and-cut methods from the literature

In this subsection, we compare the results obtained for the HpMP with the best formulation considered in this work (SANDA + ) with the results obtained with other (more sophisticated) branch-and-cut algorithms, namely, those presented by Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] and Barbato & Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF]. To allow for a more detailed comparison with the algorithm presented by Barbato & Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF], the formulation SANDA + for the HpMP was also tested on the instances considered in [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] with up to 400 nodes.

One remark is in order -the specifications of the computers in which these models were tested differ, and therefore, some of the differences in computational times may be explained by these specification differences. Regardless, we believe most of those differences to be explained by factors other than the specifications of the computers, since those specifications (in particular, regarding the CPU) are similar, with the exception of the CPU used to test Erdogan et al.'s [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF] algorithm, which is slower than the CPUs used to test the other algorithms.

The average data considered in this comparison can be found in Tables 6 and7. In Table 6, only instances with node counts between 42 and 100 were considered, since those are the instances for which every algorithm was tested. More detailed results can be found in the Appendix D. The results of Table 6 indicate that the compact formulation is not competitive for p 1 and p 5 . However, SANDA + is quite competitive for p 2 , p 3 and, to a lesser extent, p 4 . We also see that the algorithm proposed by Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] results in the lowest computational times for p 3 and p 4 , whereas the algorithm proposed by Barbato & Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] results in the lowest average computational times for every other value of p. 1403.12 (8) 2449.44 [START_REF] Marzouk | A branch-and-price algorithm for solving the hamiltonian p-median problem[END_REF] Table 7: Average CPU times obtained with algorithm F from [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] and SANDA + for all instances for which both algorithms were tested As for Table 7, it shows a clear trend: Barbato & Gouveia's [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] algorithm is substantially better if p = p 1 (that is, if it is close to 1) or if p = p 5 (that is, if it is equal or close to n 3 ), but that is not the case for the other three values of p -indeed, for p 3 and p 4 , SANDA + even results in lower CPU times (on average). We observe, however, that for the instances with more than 226 nodes, SANDA + is not competitive with Barbato & Gouveia's [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] algorithm.

In conclusion, the results obtained with the SANDA + formulation are comparable with the results obtained with the more sophisticated branch-and-cut algorithms presented by Erdogan et al. [START_REF] Erdogan | Exact and heuristic algorithms for the hamiltonian p-median problem[END_REF], Bektaş et al. [START_REF] Bektaş | Revisiting the hamiltonian p-median problem: A new formulation on directed graphs and a branch-and-cut algorithm[END_REF] and Barbato & Gouveia [START_REF] Barbato | The hamiltonian p-median problem: Polyhedral results and branch-and-cut algorithm[END_REF] for p 2 , p 3 and p 4 . Indeed, and as it can be seen in Appendix D, SANDA + was even able to solve some instances in record time -for instance, kroC100 and u159 with p = p 4 . It was also able to find new optimal solutions, for instance, for u159 with p = p 3 or d198 and pr226 with p = p 4 .

Conclusions

In this paper, we have introduced, studied and analysed several classes of compact formulations for the symmetric Hamiltonian p-Median Problem (HpMP). We have focused our study on compact formulations for eliminating solutions with less than p cycles since such formulations are less well known and studied than formulations which prevent solutions with more than p cycles. The advantage of such formulations is that they can be readily used in combination with off-the-shelf optimization software, unlike other types of formulations possibly involving the use of exponentially sized sets of variables or constraints which usually require the use of specialized methods, such as constraint separation, which may not always be easy to understand, implement or use. The proposed formulations are based on a common motivation, that is, they contain variables that assign labels to nodes, and prevent less than p cycles by stating that different depots must have different labels and that nodes in the same cycle must have the same label. We introduce and study aggregated formulations (which consider integer variables that represent the label of the node) and disaggregated formulations (which consider binary variables that assign each node to a given label). The aggregated models are new. The disaggregated formulations are not, although in all of them new enhancements have been included to make them more competitive with the aggregated models.

The main conclusion of this study is that, in the context of compact formulations, it is worth looking at the new and more compact models with the node variables. Despite the weaker LP relaxation bounds, the fewer variables and constraints lead to faster computational times, especially when solving instances with more than 50 nodes, and manage to be competitive with more sophisticated branch-and-cut algorithms from the literature for instances with up to around 200 nodes (unless p is very close to 1 or n 3 ). The statements made at the end of sections 5.2.1 and 5.2.2 suggest that further improvements for these aggregated models may be a topic of future research combined with a deeper understanding of the relationship between disaggregated and aggregated models.

A Test Results (HpMP ≥ )

This appendix includes the detailed test results obtained using the models proposed in this paper for the HpMP ≥ . A "N/A" entry in one of these tables indicates CPLEX ran out of memory while performing the corresponding test. We also recall that each test has a time limit of 1h, and that if that time limit is reached in a test, the corresponding entry in the table containing the CPU times is the gap, in %, between the obtained lower and upper bounds. Finally, for each row, the best result is in bold, unless all results are equal, in which case no result is in bold. This appendix includes the detailed test results obtained using the models proposed in this paper for the HpMP. A "N/A" entry in one of these tables indicates CPLEX ran out of memory while performing the corresponding test. We also recall that each test has a time limit of 1h, and that if that time limit is reached in a test, the corresponding entry in the table containing the CPU times is the gap, in %, between the obtained lower and upper bounds. Finally, for each row, the best result is in bold, unless all results are equal, in which case no result is in bold.

Proposition 17 .

 17 The formulation DNCA is valid for the HpMP ≥ .

Figure 1 :

 1 Figure 1: Number of solved instances of the benchmark set within fixed time thresholds for aggregated models for the HpMP ≥

Figure 2 :

 2 Figure 2: Number of solved instances of the benchmark set within fixed time thresholds for disaggregated models for the HpMP ≥

Figure 3 :

 3 Figure 3: Number of solved instances of the benchmark set within fixed time thresholds for aggregated models for the HpMP

Table 1 :

 1 Average LP and RN gaps and CPU times obtained with each aggregated model for the HpMP ≥ , for each value of p

	RN Gaps (%)	Average CPU Times (s)

Table 2 :

 2 Average LP and RN gaps and CPU times obtained with each disaggregated model for the HpMP

≥ , for each value of p

Table 3 :

 3 . Average LP and RN gaps obtained with SANDA + and SANCA + , for each value of p

			Average LP Gaps (%)		Average RN Gaps (%)
		p 1	p 2	p 3	p 4	p 5	p 1	p 2	p 3	p 4	p 5
	SANDA +	1.52% 1.17% 1.83% 3.11% 9.19% 0.94% 0.64% 1.21% 2.25% 6.75%
	SANCA +	1.52% 1.17% 1.84% 3.13% 9.32%	0.96% 0.63% 1.33% 2.43%	7.5%
					Average CPU Times (s)
				p 1	p 2	p 3	p 4	p 5
			SANDA +	417 (2) 22 (0) 70 (0) 247 (1) 2047 (10)
			SANCA +	461 (2)	33 (0) 303 (0) 462 (1)	2183 (11)

Table 4 :

 4 Average CPU times obtained with SANDA

+

and SANCA + , for each value of p

The main conclusion from these results is that the SANDA + formulation results in the best CPU times and has, on average, marginally better LP relaxation bounds than SANCA + . Regarding the LP relaxation bounds, there are no differences for the smallest values of p, and the difference between the average gaps

p

  Erdogan et al. Bektaş et al. Barbato & Gouveia SANDA

					+
	p 1	304.27 (0)	39.47 (0)	0.83 (0)	490.11 (2)
	p 2	320.06 (1)	31.18 (0)	24.09 (0)	26.06 (0)
	p 3	490.93 (2)	40.94 (0)	235.23 (1)	81.92 (0)
	p 4	273.24 (1)	60.29 (0)	513.47 (2)	290.33 (1)
	p 5	1655.14 (6)	2632.35 (9)	1266.72 (4)	2407.19 (10)

Table 6 :

 6 Average CPU times obtained with three algorithms from the literature and SANDA

	+ , for instances

Table 8 :

 8 LP gaps obtained with aggregated models for the HpMP ≥

	Model

Table 9 :

 9 RN gaps obtained with aggregated models for the HpMP ≥

	Model

Table 10 :

 10 CPU times (s) to reach optimality with aggregated models for the HpMP ≥

	Model

Table 11 :

 11 LP gaps obtained with disaggregated models for the HpMP ≥

	Model

Table 12 :

 12 RN gaps obtained with disaggregated models for the HpMP ≥

	Model

Table 13 :

 13 CPU times (s) to reach optimality with disaggregated models for the HpMP ≥

B Test Results (HpMP)

Table 14 :

 14 LP and RN Gaps (%) to reach optimality with aggregated models for the HpMP

			Model			Model
	Instance p SANDA					+
		p 1	0.5	0.3		p 1	58.3	30.1
		p 2	0.1	0.1		p 2	9.4	18.9
	gr24	p 3	0.1	0.1	st70	p 3	2.1	2.3
		p 4	1	1.7		p 4	6	14
		p 5	1.2	1.6		p 5	3.09%	4.38%
		p 1	0.8	0.4		p 1	7.3	57.9
		p 2	0.9	0.8		p 2	17.3	63.9
	fri26	p 3	0.3	0.4	eil76	p 3	370	685.3
		p 4	0.2	0.3		p 4	1.30%	1.02%
		p 5	0.1	0.1		p 5	19.01%	7.33%
		p 1	0.3	0.4		p 1	2.9	16.9
		p 2	0.2	0.2		p 2	10.5	26.6
	bayg29	p 3	0.3	0.3	pr76	p 3	379.4	1979
		p 4	5.7	6.8		p 4	55.2	88.2
		p 5	11.3	12.4		p 5	1372.7	1125.7
		p 1	2.3	1.2		p 1	5.5	3.6
		p 2	0.7	1		p 2	249.1	285.3
	swiss42	p 3	0.6	1.1	rat99	p 3	170.6	1642.3
		p 4	1.6	3.2		p 4	239.5	588.1
		p 5	5.9	8.8		p 5	2.14%	26.34%
		p 1	0.8	1.3		p 1	0.39%	0.37%
		p 2	0.5	0.8		p 2	5.5	25
	att48	p 3	5	7.3	kroA100	p 3	23.1	192.4
		p 4	12.6	16.2		p 4	471.5	1663.2
		p 5	2.13%	2.73%		p 5	29.41%	15.72%
		p 1	5.3	5.6		p 1	748.4	767.8
		p 2	0.5	0.9		p 2	95.2	104.4
	gr48	p 3	41.6	48.7	kroB100	p 3	1.9	9.6
		p 4	21.8	75.7		p 4	15.4	96.6
		p 5	261.3	341.7		p 5	6.89%	14.01%
		p 1	3.8	11.1		p 1	99	160.8
		p 2	0.3	0.5		p 2	8.1	25.9
	hk48	p 3	3.5	7.2	kroC100	p 3	116.9	675.3
		p 4	12.8	13.6		p 4	146.4	254.4
		p 5	33.8	124.6		p 5 50.20% 100.00%
		p 1	0.8	1.1		p 1	3.3	25.3
		p 2	2.6	14.1		p 2	6.7	5.9
	eil51	p 3	79.3	163.8	kroD100	p 3	124.7	231.4
		p 4	85.5	58.8		p 4	205	1495.3
		p 5	1612.6	1.50%		p 5 11.01%	28.24%
		p 1	2.1	2.2		p 1	10.7	114.3
		p 2	0.3	0.6		p 2	4.6	12.1
	berlin52	p 3	2.5	12.4	kroE100	p 3	15.2	178.2
		p 4	11.5	22.8		p 4	33.9	449.3
		p 5	1167.7	1981.5		p 5	3.14%	100.00%
		p 1	1.05%	0.92%		p 1	181.3	815.7
		p 2	20.4	22.6		p 2	11.2	51.9
	brazil58	p 3	1	2.2	rd100	p 3	55.1	213.8
		p 4	1.9	4.7		p 4	14.8	780.8
		p 5	468.3	458.7		p 5	5.85%	17.86%

+ SANCA + Instance p SANDA + SANCA

Table 15 :

 15 CPU times (s) to reach optimality with aggregated models for the HpMP p Erdogan et al. Bektaş et al. Barbato & Gouveia SANDA

	+

Table 18 :

 18 CPU times obtained with three algorithms from the literature and SANDA

+ (part 2)

Acknowledgements

L. Gouveia and F. Canas were supported by Portuguese National Funding from FCT -Foundation for Science and Technology, under projects UIDB/04561/2020 and UIDP/04561/2020.

C Test Results (DL and MCF comparison) This appendix includes the detailed test results obtained using the adapted DL and MCF formulations, with SANDA used to prevent less than p cycles. Only instances with n ≤ 76 were considered, and only the LP gaps were calculated. One last observation is in order: we were so far unable to determine the optimal values of instances brazil58 (with p = 4) and pr76 (with p = 2). However, we have obtained intervals to which the optimal values are guaranteed to belong. They are [21898.1533, 22211] and [104500, 106362.1561]. To calculate the gaps for these two instances, we have used the upper bounds of the corresponding intervals. Even if these upper bounds are not guaranteed to be the optimal values, the differences in the average gaps should not be substantial, as all other gaps were calculated using the corresponding optimal values.