Thomas Kleinert

Martine Labbé
email: martine.labbe@ulb.ac.be

Fränk Plein
email: frank.plein@ulb.ac.be

Martin Schmidt
email: martin.schmidt@uni-trier.de

There's No Free Lunch: On the Hardness of Choosing a Correct Big-M in Bilevel Optimization

Keywords: Mathematics Subject Classification. 90C46, 90C11, 91A20, 91A40 Bilevel optimization, Mathematical programs with complementarity constraints (MPCC), Bounding polyhedra, Big-M, Hardness

One of the most frequently used approaches to solve linear bilevel optimization problems consists in replacing the lower-level problem with its Karush-Kuhn-Tucker (KKT) conditions and by reformulating the KKT complementarity conditions using techniques from mixed-integer linear optimization. The latter step requires to determine some big-M constant in order to bound the lower level's dual feasible set such that no bilevel-optimal solution is cut off. In practice, heuristics are often used to find a big-M although it is known that these approaches may fail. In this paper, we consider the hardness of two proxies for the above mentioned concept of a bilevel-correct big-M . First, we prove that verifying that a given big-M does not cut off any feasible vertex of the lower level's dual polyhedron cannot be done in polynomial time unless P = NP. Second, we show that verifying that a given big-M does not cut off any optimal point of the lower level's dual problem (for any point in the projection of the high-point relaxation onto the leader's decision space) is as hard as solving the original bilevel problem.

Introduction

A bilevel optimization problem consists in a constrained optimization problem in which some constraints specify that a subset of variables constitutes an optimal solution of a second (auxiliary) optimization problem. Since the publication of the first and seminal paper [START_REF] Bracken | Mathematical programs with optimization problems in the constraints[END_REF], research on the subject has become increasingly important. Indeed, the bilevel structure allows the modeling of a large number of real-life problems involving two types of decision makers, a leader and a follower (or several followers) interacting hierarchically. Such optimization problems appear in many fields of application like energy markets [1, 10, 18-21, 23, 25], critical infrastructure defense [START_REF] Caprara | Bilevel Knapsack with Interdiction Constraints[END_REF][START_REF] Casorrán | A study of general and security Stackelberg game formulations[END_REF][START_REF] Denegre | Interdiction and discrete bilevel linear programming[END_REF][START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving Bayesian Stackelberg games[END_REF], or pricing [START_REF] Labbé | A bilevel model of taxation and its application to optimal highway pricing[END_REF][START_REF] Labbé | Bilevel programming and price setting problems[END_REF][START_REF] Rusmevichientong | A nonparametric approach to multiproduct pricing[END_REF].

Due to their ability of modeling hierarchical decision processes, bilevel optimization problems are inherently hard to solve. In [START_REF] Dempe | Bilevel Programming Problems[END_REF][START_REF] Hansen | New branch-and-bound rules for linear bilevel programming[END_REF] it is shown that even the easiest instantiation, i.e., bilevel problems with linear upper and lower level, is strongly NP-hard. Moreover, even checking local optimality for a given point is NP-hard as well [START_REF] Vicente | Descent approaches for quadratic bilevel programming[END_REF]. For other hardness results we refer to, e.g., [START_REF] Bard | Some properties of the bilevel programming problem[END_REF]. For general surveys of bilevel optimization see [START_REF] Bard | Practical bilevel optimization: algorithms and applications[END_REF][START_REF] Dempe | Bilevel Optimization: Theory, Algorithms, and Applications[END_REF][START_REF] Dempe | Foundations of bilevel programming[END_REF][START_REF] Dempe | Bilevel Programming Problems[END_REF] and [START_REF] Wen | Linear Bi-Level Programming Problems -A Review[END_REF] for a survey focusing on linear-linear (LP-LP) bilevel problems.

In this paper, we consider the LP-LP bilevel problem min x∈R n ,y∈R m c x + d y (1a)

s.t. Ax + By ≥ a, (1b)
y ∈ arg min

ȳ {f ȳ : Cx + Dȳ ≥ b} (1c) with c ∈ R n , d, f ∈ R m , A ∈ R k×n , B ∈ R k×m , a ∈ R k , C ∈ R ×n , D ∈ R ×m
, and b ∈ R . Problem (1a)-(1b) is the so-called upper level and (1c) is the lower level.

Here, we consider the optimistic version of bilevel optimization [START_REF] Dempe | Foundations of bilevel programming[END_REF]. This means that whenever the lower-level problem has multiple solutions y, the leader chooses the most favorable one in terms of the upper level's problem. LP-LP bilevel problems are often solved by a reformulation to an equivalent singlelevel problem. Usually, this is done by one of the following two approaches. One can replace the lower level with its primal and dual feasibility conditions as well as the strong-duality equation or one replaces the lower level with its Karush-Kuhn-Tucker (KKT) conditions. Both approaches have their drawbacks: The strong-duality based technique yields nonconvex bilinear terms whereas the KKT approach leads to the following mathematical program with complementarity constraints (MPCC):

min x,y,λ c x + d y (2a) s.t. Ax + By ≥ a, D λ = f, (2b)
0 ≤ Cx + Dy -b ⊥ λ ≥ 0. (2c)
Here, for two non-negative vectors 0 ≤ a, b ∈ R n , a ⊥ b abbreviates a i b i = 0 for all i = 1, . . . , n. The single-level reformulation (2) of the bilevel problem (1) explicitly models the optimistic version of bilevel optimization by construction. Consequently, there is no need anymore to distinguish between y and ȳ. In what follows, we thus call x the upper-level and y the lower-level (primal) variables. The lower-level dual variables are given by λ ∈ R . The hardness of Problem (2) stems from the complementarity conditions (2c). Since these constraints can be reformulated in a mixed-integer linear way, the KKT approach is often preferred in practice. Typically, one applies the reformulation introduced in [START_REF] Fortuny-Amat | A Representation and Economic Interpretation of a Two-Level Programming Problem[END_REF], which requires an additional binary variable z i ∈ {0, 1} for every i ∈ {1, . . . , } and the additional constraints

λ i ≤ M d z i , (Cx + Dy -b) i ≤ (1 -z i )M p , i ∈ {1, . . . , }, (3)
where M d and M p are sufficiently large constants, called big-M s. In this note, we focus on the big-M for bounding the lower-level dual variables, i.e., on M = M d . Since we show that finding M d is hard, this obviously implies that finding M d and M p is hard as well. Applying (3) requires to bound the lower level's dual polyhedron such that no point λ * that is part of an optimal solution (x * , y * , λ * ) of ( 2) is cut off. Stated differently, one needs to choose an M that preserves all bilevel-optimal points (x * , y * ). We call an M with this property a bilevel-correct big-M . When the dual of the lower level has a finite optimal value, there exists an optimal solution λ * that is a vertex of the associated feasible polyhedron. Hence, it is sufficient to obtain bounds on the dual variables that-independently of the upper-level decision-do not cut off (i) any feasible vertex of the lower level's dual polyhedron or (ii) any optimal vertex of the lower level's dual polyhedron. Note that these requirements do not take into account bilevel optimality but still preserve all optimal solutions (x * , y * , λ * ) of (2), i.e., (i) and (ii) still yield bilevel-correct big-M 's. This means, that (i) and (ii) are sufficient (but not necessary) conditions for a big-M to be bilevel-correct.

The choice of the big-M is often done heuristically, which may result in a severe issue: If the big-M is not chosen large enough, a "solution" of (2) with (2c) replaced by [START_REF] Bard | Practical bilevel optimization: algorithms and applications[END_REF] does not need to be a bilevel-optimal point. In fact, this point does not even need to be bilevel-feasible. See, e.g., [START_REF] Pineda | Solving Linear Bilevel Problems Using Big-Ms: Not All That Glitters Is Gold[END_REF], where a common heuristic for computing a big-M is shown to deliver wrong results.

The contribution of this note is twofold. First, in Section 2, we consider the hardness of verifying that a given big-M does not cut off any feasible vertex of the lower level's dual polyhedron. We show that there is no polynomial-time algorithm for this verification unless P = NP. Second, in Section 3, we show that validating that a given big-M does not cut off any optimal point of the lower level's dual problem (for any given feasible upper-level variable x) is as hard as solving the original bilevel problem. Both results together imply that there is no hope for an efficient, i.e., polynomial-time, general-purpose method for verifying or computing a correct big-M in bilevel optimization unless P = NP. Thus, our results strongly indicate that problem-specific bounds on the lower level's dual variables need to be investigated if the given bilevel problem is going to be solved using the KKT approach combined with the classical big-M linearization of KKT complementarity conditions.

Hardness of Bounding the Vertices of an Unbounded Polyhedron

Whenever the bilevel problem ( 1) is feasible, the lower-level primal and dual problem have a finite optimal solution. In particular, there is a vertex of the feasible region of the lower-level dual problem at which the optimal dual solution is attained. Thus, one way of preserving every bilevel-optimal solution in the KKT reformulation ( 2) is to choose a big-M such that no lower-level dual vertex is cut off. This bounding approach yields a bilevel-correct big-M . In this section, we show-even more generally-that bounding the vertices of an unbounded polyhedron is hard. This result is then applied to the lower level's dual polyhedron. Since the hardness result is mainly based on the unboundedness of this polyhedron, the question arises whether this situation frequently appears in practical LP-LP bilevel problems. It turns out that this is the case for almost all instances of bilevel test sets from the literature-which is also supported by the theoretical results in [START_REF] Williams | Boundedness relations for linear constraint sets[END_REF].

To obtain a hardness result in the Turing model of computation, we assume that all problem data are rational and thus are Turing representable. Let P (A, b) := {x ∈ Q n : Ax ≤ b} be an unbounded polyhedron defined by A ∈ Q k×n and b ∈ Q k . For M ∈ Q and j ∈ {1, . . . , n}, let Q j (A, b, M ) := {x ∈ Q n : Ax ≤ b, x j ≤ M } be the polyhedron obtained from adding the bound x j ≤ M to P (A, b). To validate a given big-M , we need to verify that for every j ∈ {1, . . . , n} the bound x j ≤ M is satisfied by all vertices of P (A, b). This results in the following decision problem.

Component-wise valid bound for the vertices of a polyhedron (CVBVP).

Input:

A ∈ Q k×n , b ∈ Q k , j ∈ {1, . . . , n}, M ∈ Q. Question: Does v ∈ Q j (A, b, M ) hold for every vertex v of P (A, b)?
We will see in the following that validating a big-M is related to the problem of finding an optimal vertex v in an unbounded polyhedron with respect to a linear objective function h v. If the polyhedron is bounded at least in the direction of optimization, then this problem is equivalent to linear optimization. However, in the general case of polyhedra that are unbounded in the direction of optimization, this is a difficult task. As shown in [START_REF] Fukuda | Analysis of Backtrack Algorithms for Listing All Vertices and All Faces of a Convex Polyhedron[END_REF], the decision problem that studies the existence of a vertex of a given polyhedron such that the corresponding objective function value is larger or equal to a certain threshold K is strongly NP-complete. The proof is based on a reduction from the Hamiltonian path problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]Problem GT39] and can easily be extended to the decision problem that decides whether a vertex with an objective function value strictly larger than a certain threshold exists:

Optimal vertex of a polyhedron (OVP). Input:

A ∈ Q k×n , b ∈ Q k , h ∈ Q n , K ∈ Q. Question: Is there a vertex v of P (A, b) with h v > K?
As pointed out above, w.r.t. the linearization (3), we are interested in the special case h = e j . The related decision problem is the following:

Component-wise optimal vertex of a polyhedron (COVP). Input:

A ∈ Q k×n , b ∈ Q k , j ∈ {1, . . . , n}, K ∈ Q. Question: Is there a vertex v of P (A, b) with v j > K?
We now show that even for this subclass of instances, the decision problem is strongly NP-complete. In what follows, vert(P (A, b)) denotes the set of vertices of the polyhedron P (A, b).

Theorem 1. COVP is strongly NP-complete.

Proof. We prove the result for j = 1. For any other j ∈ {2, . . . , n}, e 1 can be replaced with e j in the proof.

It is clear that the problem is in NP. We prove its hardness by reduction from

OVP. Let A ∈ Q k×n , b ∈ Q k , h ∈ Q n , K ∈ Q be a
given OVP instance and assume h = 0. Otherwise, the corresponding instance is trivial. Now take j ∈ {1, . . . , n} with h j = 0. We construct a basis of Q n by replacing e j with h. If we put h as the first basis vector, the corresponding linear transformation is given by the inverse of matrix B = [h, e 1 , . . . , e j-1 , e j+1 , . . . , e n ] and can be computed in polynomial time. Using this basis change, we now linearly transform the hyperplanes defining the polyhedron P (A, b) and the objective function vector h of the given OVP instance.

We construct an instance of COVP by à = AB -, b = b, K = K. Note that B -1 h = e 1 holds. It remains to show that there exists a vertex v of P (A, b) with h v > K if and only if there exists a vertex ṽ of P ( Ã, b) with ṽ1 > K. Let v ∈ vert(P (A, b)) such that h v > K and define ṽ := B v. Then,

Ãṽ = AB -B v = Av ≤ b = b, ṽ1 = h B -B v = h v > K = K.
Thus, ṽ ∈ P ( Ã, b) and it is clear that ṽ is also a vertex of P (A, b). Conversely, let ṽ ∈ vert(P ( Ã, b)) with ṽ1 > K and define v := B -ṽ. Then,

Av = AB -ṽ = Ãṽ ≤ b = b, h v = h B -ṽ = ṽ1 > K = K.
Thus, v is a vertex of P (A, b).

Using problem COVP, we can deduce the complexity of CVBVP.

Theorem 2. CVBVP is strongly coNP-complete.

Proof. The decision problem CVBVP, i.e., the complement of CVBVP, is to find a vertex v ∈ vert(P (A, b)) such that v j > M holds. This is equivalent to COVP with K = M .

Finally, we can state the main result of this section. Proof. Assume a polynomial-time algorithm exists. Then, for every j ∈ {1, . . . , n} we can efficiently decide whether vert(P (A, b)) ⊆ Q j (A, b, M ) holds. This implies that we can decide CVBVP in polynomial time, and thus P = coNP must also hold. Since P is closed under taking the complement, it follows that P = NP.

Corollary 1. Let A ∈ Q k×n , b ∈ Q k ,
As a final remark, note that to compute the tightest possible big-M such that no vertex of P (A, b) is cut off, we can set

M := max j∈{1,...,n} max x∈vert(P (A,b)) x j . (4)
It is equivalent to solving COVP for j ∈ {1, . . . , n} and taking the maximum value. Thus, ( 4) cannot be computed in polynomial time, unless P = NP.

Valid Bounds for Bilevel-Feasible Solutions

Recall that a big-M is bilevel-correct, if it does not cut off any bilevel-optimal solution. For this, it is sufficient to find a big-M that maintains at least one optimal lower-level dual vertex for every feasible upper-level decision. This means that, in contrast to the big-M of Section 2, we now allow to cut off lower-level dual vertices that do not correspond to an optimal solution.

Here and in what follows we denote the high-point relaxation of the bilevel problem [START_REF] Ambrosius | Endogenous Price Zones and Investment Incentives in Electricity Markets: An Application of Multilevel Optimization with Graph Partitioning[END_REF] as

H := {(x, y) ∈ R n × R m : Ax + By ≥ a, Cx + Dy ≥ b}
and the corresponding projection onto the space of x-variables is defined as

H x := {x ∈ R n : ∃y with (x, y) ∈ H}.
For the sake of simplicity, we make the following assumption. Assumption 1. For every upper-level decision x ∈ H x , the lower-level problem (1c) admits a unique solution y and satisfies the linear independence constraint qualifications (LICQ) at y.

The assumption of a unique lower-level solution is justified in this scope given that the bilevel problem becomes even more difficult to analyze otherwise; see, e.g., Chapter 7 in [START_REF] Dempe | Foundations of bilevel programming[END_REF]. Moreover, the LICQ guarantees the uniqueness of the lower-level dual optimal solution for every upper-level decision x ∈ H x ; see, e.g., Chapter 12 of [START_REF] Nocedal | Numerical Optimization[END_REF]. Let us note that the assumption that LICQ holds at the unique primal solution of the lower level can be very strong. However, we will show the hardness of choosing a big-M that does not cut off any optimal vertex of the lower level's dual polyhedron under the simplifying Assumption 1. Thus, the problem of choosing such a big-M is hard also for the situation in which Assumption 1 is dropped.

We start by introducing a validity criterion for the big-M proxy discussed in this section. To this end, define the lower-level optimal value function ϕ(x) for any upper-level decision x ∈ H x by means of its dual as

ϕ(x) := max λ {(b -Cx) λ : D λ = f, λ ≥ 0}. (5)
Further, for any upper-level decision x ∈ H x , M ∈ R, and i ∈ {1, . . . , }, let ϕ i (x, M ) be the optimal value function of the lower level's dual problem with the additional bound λ i ≤ M , i.e.,

ϕ i (x, M ) := max λ {(b -Cx) λ : D λ = f, λ ≥ 0, λ i ≤ M }, (6)
where we formally set 6) is infeasible. Under Assumption 1, all bilevel-feasible solutions remain the same after adding the big-M bounds to the lower level's dual problem if and only if for every upper-level decision x ∈ H x and for every i ∈ {1, . . . , }, the lower-level optimal value stays unchanged, i.e., if ϕ(x) = ϕ i (x, M ) holds. We now collect some simple observations on these two optimal value functions that are used afterward.

ϕ i (x, M ) = -∞ if Problem (
Observation 1. Given an upper-level decision x ∈ H x and i ∈ {1, . . . , }, the following properties hold:

(a) ϕ(x) ≥ ϕ i (x, M ) for every M ∈ R.

(b) ϕ i (x, •) is monotonically increasing. (c) Suppose there exists an M ∈ R with ϕ(x) = ϕ i (x, M ). Then, ϕ(x) = ϕ i (x, M ) holds for every M ≥ M .

Lemma 1. Suppose that Assumption 1 holds and let an upper-level decision x ∈ H x and M = M (x) ∈ R be given. Then, for every i ∈ {1, . . . , }, ϕ(x) = ϕ i (x, M (x)) holds if and only if M (x) ≥ max{λ * i (x) : i ∈ {1, . . . , }}, where λ * (x) is the unique optimal solution of the lower level's dual problem (5) corresponding to x.

Proof. If M (x) < max{λ * i (x) : i ∈ {1, . . . , }}, then there is an i ∈ {1, . . . , } such that the optimal solution of the lower level's dual problem is cut off by the bound λ i ≤ M (x), which again is equivalent to ϕ(x) > ϕ i (x, M (x)).

In particular, this implies that for every fixed upper-level decision, we can validate a given big-M by computing the corresponding unique optimal solution of the lower level's dual problem and by verifying that it satisfies the bounds λ i ≤ M for all i ∈ {1, . . . , }.

For the case that all upper-level decisions are taken into account, the next result gives a necessary and sufficient condition for the property that a big-M does not cut off any bilevel-feasible point.

Proof. Observe that the first constraint of (2b) defines the domain of the upperlevel decisions x, whereas the second constraint together with (2c) determine the lower-level primal-dual optimal solution (y, λ) corresponding to x. The final result then follows by Lemma 1 and Property (c) in Observation 1.

Theorem 3 implies that validating a big-M requires optimizing different objective functions over a set of constraints that are equivalent to feasibility of the original bilevel problem. In [START_REF] Karp | On the Computational Complexity of Combinatorial Problems[END_REF], linear 0-1-feasibility has been shown to be NP-complete. It is thus possible to adapt the techniques from [START_REF] Audet | Links Between Linear Bilevel and Mixed 0-1 Programming Problems[END_REF] to show the NP-completeness of LP-LP bilevel feasibility by reduction from linear 0-1-feasibility. Similarly to Corollary 1, we can thus state that there is no polynomial-time validation of a given big-M w.r.t. [START_REF] Bracken | Mathematical programs with optimization problems in the constraints[END_REF] unless P = NP. On the other hand, computing the tightest big-M w.r.t. the proxy considered in this section requires solving a maximization problem over all bilevel-feasible solutions for every i ∈ {1, . . . , } and taking the maximum objective value. Computing this big-M is therefore as hard as solving the initial problem and there is little hope of doing it efficiently, unless the original bilevel problem (1) can be solved in polynomial time.

Conclusion

Many applications of LP-LP bilevel optimization make use of the KKT reformulation of the lower-level problem together with a big-M linearization of the KKT complementarity constraints. This results in a single-level mixed-integer linear problem that can, in principle, be solved with state-of-the-art solvers. However, to guarantee bilevel feasibility of a solution obtained by this approach, one needs to validate the bilevel-correctness of the big-M that is used to bound the lower level's dual variables-a necessary task that is not always carried out in practice.

In general, such a big-M is bilevel-correct if it does not cut off any bilevel-optimal point. In this note we considered two proxies for this type of correctness and proved that even validating that a given big-M does not cut off any feasible or optimal vertex of the lower level's dual polyhedron cannot be done in polynomial time unless P = NP. Both proxies abstract from upper-level optimality and, thus, are only sufficient but not necessary conditions for a big-M to be bilevel-correct. Hence, validating that a given big-M preserves all bilevel-optimal points can be expected to be at least as hard since one needs to take into account another-i.e., the upper level's-optimization problem on top of what needs to be considered for the two proxies.

Our results strongly suggest that the popular big-M approach needs to be applied very carefully. If the bilevel-correctness of the chosen big-M is not guaranteed by problem-specific insights, it cannot be formally guaranteed that the obtained "solutions" are indeed bilevel-optimal. In such cases, we suggest to better resort to exact approaches that do not rely on big-M 's like, e.g., the kth best algorithm ( [START_REF] Bialas | Two-level linear programming[END_REF] or [START_REF] Bard | Practical bilevel optimization: algorithms and applications[END_REF]Chapter 5.3.1]) or branch-and-bound methods ( [START_REF] Bard | A branch and bound algorithm for the bilevel programming problem[END_REF][START_REF] Hansen | New branch-and-bound rules for linear bilevel programming[END_REF], or [START_REF] Bard | Practical bilevel optimization: algorithms and applications[END_REF]Chapter 5.3.2]). Moreover, identifying reasonably generic sub-classes of bilevel optimization problems for which it is easy to determine a bilevel-correct big-M is subject to future research.

  and M ∈ Q. Then, there exists no polynomialtime algorithm for checking whether vert(P (A, b)) ⊆ n j=1 Q j (A, b, M ), unless P = NP.

Theorem 3 .

 3 Let M ∈ R be given and suppose that Assumption 1 holds. Then, for every upper-level decision x ∈ H x and for every i ∈ {1, . . . , }, ϕ(x) = ϕ i (x, M ) holds if and only if M ≥ max i∈{1,..., } max x,y,λ {λ i : (2b), (2c)} .

Acknowledgments

This research has been performed as part of the Energie Campus Nürnberg and is supported by funding of the Bavarian State Government. The authors thank the DFG for their support within project A05, B08, and Z01 in CRC TRR 154. Fränk Plein thanks the "Fonds de la Recherche Scientifique" (F.R.S.-FNRS) for financial support. Martine Labbé has been partially supported by the Fonds de la Recherche Scientifique -FNRS under Grant(s) no PDR T0098.18.


Víctor Bucarey
email: vbucarey@vub.be

Bernard Fortz
email: bernard.fortz@ulb.ac.be

Natividad González-Blanco

Martine Labbé
email: mlabbe@ulb.ac.be

Juan A Mesa
email: jmesa@us.es

Benders decomposition for Network Design Covering Problems

Keywords: Benders decomposition, Network Design, Rapid Transit Network

We consider two covering variants of the network design problem. We are given a set of origin/destination (O/D) pairs and each such O/D pair is covered if there exists a path in the network from the origin to the destination whose length is not larger than a given threshold. In the first problem, called the maximal covering network design problem, one must determine a network that maximizes the total demand of the covered O/D pairs subject to a budget constraint on the design costs of the network. In the second problem, called the partial covering network design problem, the design cost is minimized while a lower bound is set on the total demand covered.

After presenting formulations, we develop a Benders decomposition approach to solve the problems.

Further, we consider two different stabilization methods to determine the Benders cuts as well as the addition of cut-set inequalities to the master problem. Computational experiments show the efficiency of these different aspects.

Introduction

Infrastructure network design constitutes a major step in the planning of a transportation network since the performance, efficiency, robustness, and other features strongly depend on the selected points and the way of connecting them, see [START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF]. For instance, the main purpose of a rapid transit network is to improve the mobility of the inhabitants of a city or metropolitan area. This improvement could lead to lower journey times, less pollution and/or less energy consumption which drives the communities to a more sustainable mobility.

Since it is generally too expensive to connect all the existing facilities, one must determine a subnetwork that serves at best the traffic demand. Depending on the application, different optimality measures are considered. In particular, in the field of transportation, and especially in the area of passengers transportation, the aim is to get the infrastructure close to potential customers. In this framework, [START_REF] Schmidt | Location of speed-up subnetworks[END_REF] propose to minimize the maximum routing cost for an origindestination pair when using the new network. Alternatively, the traffic between an origin and a destination may be considered as captured if the cost or travel time when using the network is not larger than the cost or travel time of the best alternative solution (not using the new network). In this case, [START_REF] Perea | Transportation infrastructure network design in the presence of modal competition: computational complexity classification and a genetic algorithm[END_REF] and [START_REF] García-Archilla | Grasp algorithms for the robust railway network design problem[END_REF] propose to select a sub(network) from an underlying network with the aim of capturing or covering as much traffic for a reasonable construction cost. This paper is devoted to this problem, called the Maximum Covering Network Design Problem (MC) as well as to the closely related problem called, Partial Covering Network Design Problem (PC), in which one minimizes the network design cost for building the network under the constraint that a minimum percentage of the total traffic demand is covered.

When designing an infrastructure network, the demand is given by pairs of origin-destination points, called O/D pairs, and each such pair has an associated weight indicating the traffic between the origin and the destination. Usually this demand is encoded using an origin-destination matrix.

When planning a new network, often there exists a network already functioning and offering its service to the same set of origin-destination pairs. For example, in order to improve the mobility of a big city or metropolitan area, a new rapid transit systems is planned. The current transit system could be more dense than the rapid transit planned but slower since uses the same right-of-way than the private traffic. Thus, in some way both systems compete with each other and both compete with the private mode of transportation. A similar effect occurs with mobile telecommunication operators.

Therefore, the traffic between an origin and a destination is distributed among the several systems that provide the service.

There are mainly two ways of allocating the share of each mode. The first one is the binary all-ornothing way, where the demand is covered by each mode only if the mode covers the demand point within a range of quality service, as in [START_REF] Church | The maximal covering location problem[END_REF]. The second one is some continuous function, for example multi-logit probability distributions, as in [START_REF] Cascetta | Transportation systems analysis: models and applications volume 29[END_REF]. Both mode-share are based on the comparison of distances, times, costs, generalized costs or utilities. In this paper, we consider a binary one, where each O/D pair is covered only if the time spent into the network is below a threshold. This threshold represents the comparison between the time spent in the proposed network and a private mode assigning the full share to the most beneficial one.

Since most network design problems are NP-hard (see e.g. [START_REF] Perea | Transportation infrastructure network design in the presence of modal competition: computational complexity classification and a genetic algorithm[END_REF]), recent research efforts have been oriented to apply metaheuristic algorithms to obtain good solutions in a reasonable computational time. Thus, in the field of transportation network design, Genetic Algorithms [START_REF] Król | The design of a metro network using a genetic algorithm[END_REF]), Greedy Randomized Adaptive Search Procedures [START_REF] García-Archilla | Grasp algorithms for the robust railway network design problem[END_REF]), Adaptive Large Neighborhood Searchs [START_REF] Canca | An adaptive neighborhood search metaheuristic for the integrated railway rapid transit network design and line planning problem[END_REF]) and Matheuristics [START_REF] Canca | Integrated railway rapid transit network design and line planning problem with maximum profit[END_REF] have been applied to rapid transit network design problems with applications to medium-size instances.

In this article, after presenting models for problems (MC) and (PC), we propose exact methods based on Benders decomposition. This type of decomposition has been applied to many problems in different fields, see [START_REF] Rahmaniani | The benders decomposition algorithm: A literature review[END_REF] for a recent literature review on the use of Benders decomposition in combinatorial optimization. One remarkable recent contribution applied to set covering and maximal covering location problems appears in [START_REF] Cordeau | Benders decomposition for very large scale partial set covering and maximal covering location problems[END_REF].

Benders decomposition for network design problems have been studied since the 80s. In [START_REF] Magnanti | Tailoring benders decomposition for uncapacitated network design[END_REF], authors minimize the total building cost of an uncapacitated network subject to the constraints that all O/D pairs must be covered. Given the structure of the problem, the Benders reformulation is stated with one sub-problem for each O/D pair. A Benders decomposition for a multi-layer network design problem is presented in [START_REF] Fortz | An improved benders decomposition applied to a multilayer network design problem[END_REF]. Benders decomposition were also applied in [START_REF] Botton | Benders decomposition for the hop-constrained survivable network design problem[END_REF] in the context of designing survivable networks. Authors in [START_REF] Costa | Benders, metric and cutset inequalities for multicommodity capacitated network design[END_REF] have studied multi-commodity capacitated network design and the strength of different Benders cuts; in [START_REF] Marín | Urban rapid transit network design: accelerated benders decomposition[END_REF] a multi-objective approach is solved through Benders decomposition where coverage is maximized, but at the same time total cost design is minimized. To best of our knowledge, this is the first time that a branch-and-Benders-cut approach is applied to network design coverage problems. We also study for the first time in this context the inclusion of facet-defining cuts [START_REF] Conforti | Facet" separation with one linear program[END_REF].

This paper provides several contributions. First we present new mathematical integer formulations for the network design problems (MC) and (PC). The formulation for (MC) is stronger than a previously proposed one, see e.g. [START_REF] Marín | Urban rapid transit network design: accelerated benders decomposition[END_REF] and [START_REF] García-Archilla | Grasp algorithms for the robust railway network design problem[END_REF] (although the proposed formulation is not the main purpose of these papers), while (PC) was never studied to the best of our knowledge. Our second contribution consists of polyhedral properties that are useful from the algorithmic point of view. A third contribution is the study of exact algorithms for the network design based on different Benders implementations. We propose a normalization technique and we study the facet-define cuts. Our computational experiments show that our Benders implementations are competitive with exact and non-exact methods in the literature.

All our computational experiments were performed on a computer equipped with a Intel Core i5-7300 CPU processor, with 2.50 gigahertz 4-core, and 16 gigabytes of RAM memory. The operating system is 64-bit Windows 10. Codes were implemented in Python 3.8. These experiments have been carried out through CPLEX 12.10 solver, named CPLEX, using its Python interface. CPLEX parameters were set to their default values and the models were optimized in a single threaded mode. We used their LazyConstraintCallback function to separate integer solution and the UserCutCallback to separate fractional ones.

In tables reporting these results, t denotes the average values for solution times given in seconds, gap denotes the relative optimality gap in percent, LP gap denotes the LP gap in percent and cuts is the number of cuts generated.

The structure of the paper is as follows. In Section 2, we state the mixed integer formulation for (MC) and (PC). We also study some polyhedral properties of the formulation and propose a simple algorithm to find an initial feasible solution for both problems. In Section 3, we study different Benders implementations and some algorithmic enhancements. Also we discuss some improvements based on cut-set inequalities. A computational study is detailed in Section 4. Finally, our conclusions are presented in Section 5.

Problem formulations and some properties

In this section we present mixed integer formulations for the Maximal Covering Network Design Problem (MC) and the Partial Set Covering Network Design Problem (PC). We also describe some pre-processing methods and finish with some polyhedral properties. We first introduce some notation.

We consider an undirected graph denoted by N = (N, E), where N and E are the sets of potential nodes and edges that can be built. Each element e ∈ E is denoted by {i, j}, with i, j ∈ N . We use the notation i ∈ e if node i is a terminal node of e.

The mobility patterns in the metropolitan area are represented by a set W ⊂ N × N of origin/destination pairs named as O/D pairs. Each w = (w s , w t ) ∈ W is defined by an origin node w s ∈ N , a destination node w t ∈ N , an associated demand g w > 0 and a private utility u w > 0. This parameter translates the fact that there already exists a different mode of transportation, referenced as private mode, competing with the network to be built in an all or nothing way. In other words, an O/D pair (w s , w t ) will travel on the newly built network if it contains a path between w s and w t of length shorter than the private utility u w . We then say that the O/D pair is covered. Costs for building nodes, i ∈ N , and edges, e ∈ E, are denoted by b i and c e , respectively. The total building cost cannot exceed the budget C max .

For each e = {i, j} ∈ E, we define two arcs: a = (i, j) and â = (j, i). The resulting set of arcs is denoted by A. The length of arc a ∈ A is denoted by d a . For each O/D pair w ∈ W we define a subgraph N w = (N w , E w ) containing all feasible nodes and edges for w, i.e. that belong to a path in N whose total length is lower than or equal to u w . We also denote A w as the set of feasible arcs. In Section 2.2, we describe how to construct these subgraphs. We use notation δ + w (i) (δ - w (i) respectively) to denote the set of arcs going out (in respectively) of node i ∈ N w . In particular, δ - w (w s ) = ∅ and δ + w (w t ) = ∅. We also denote by δ w (i) the set of edges incident to node i in graph N w .

Mixed Integer Formulations

We first present a formulation of the Maximal Covering Network Design Problem (MC), whose aim is to design an infrastructure network maximizing the total demand covered subject to a budget constraint:

(MC) max

a∈δ + w (i) f w a - a∈δ - w (i) f w a =            z w , if i = w s , -z w , if i = w t , 0, otherwise, w ∈ W, i ∈ N w , (2.4) f w a + f w â ≤ x e , w ∈ W, e = {i, j} ∈ E w : a = (i, j), â = (j, i), (2.5) a∈A w d a f w a ≤ u w z w , w ∈ W, (2.6) y i , x e , z w ∈ {0, 1}, i ∈ N, e ∈ E w , w ∈ W, (2.7)
f w a ∈ {0, 1}, a ∈ A w , w ∈ W, (2.8)
where y i and x e represent the binary design decision of building node i and edge e, respectively. Mode choice variables z w take value 1 if the O/D pair w is covered and 0 otherwise. Variables f w a are used to model a path between w s and w t , if possible. Variable f w a takes value 1 if arc a belongs to the path from w s to w t , and 0 otherwise. For each f w a such that a / ∈ A w , this variable is set to zero.

The objective function (2.1) maximizes the demand covered. Constraint (2.2) limits the total building cost. For each pair w, expressions (2.4), (2.5) and (2.6) guarantee demand conservation and link flow variables f w a with decision variables z w and design variables x e . Constraints (2.6) upper bound the length of the path for each pair w. Variable z w will take value 1 only if there exists a path between w s and w t with length shorter than u w . This path is represented by variables f w a . Finally, constraints (2.7) and (2.8) state that variables are binary.

The Partial Covering Network Design Problem (PC), which minimizes the total building cost of the network subject to a minimum coverage level of the total demand, can be formulated as follows:

(PC) min where β ∈ (0, 1] and Z total = w∈W g w . Here, the objective function (2.9) minimizes the design cost.

Constraint (2.10) imposes that a proportion β of the total demand is covered.

In previous works by [START_REF] Marín | Urban rapid transit network design: accelerated benders decomposition[END_REF] and [START_REF] García-Archilla | Grasp algorithms for the robust railway network design problem[END_REF] , constraints

(2.5) and (2.6) are written in a different way. For example, in García-Archilla et al. ( 2013), these constraints were written as

f w a + z w -1 ≤ x a , w ∈ W, e = {i, j} ∈ E w : a = (i, j), (2.11) a∈A w d a f w a + M (z w -1) ≤ u w z w , w ∈ W, (2.12)
where the design variable x a is defined per arc. Given that z w -1 ≤ 0, expressions (2.5) and (2.6) are stronger than (2.11) and (2.12), respectively.

In addition, constraint (2.12) involves a "big-M" constant. Our proposed formulation does not need it, which avoids the numerical instability generated by this constant. In Table 1, we compare both formulations for (MC). We use CPLEX to solve some instances described in Subsection 4.1. Average solution times in seconds and percent LP gaps are shown for instances with 10 and 20 nodes. We also tested instances with 40 nodes but most of them were not solved to optimality within one hour. In that case, we provide the optimality gap instead of the solution time. We consider 5 instances of each size. Note that constraints (2.12) are equivalent to constraints (2.6) by setting M = 0. We tested several positive values for M . We observed that our proposed formulation is not only stronger than the one proposed in García-Archilla et al. ( 2013), but it is also more computationally efficient.

Preprocessing methods

In this section we describe some methods to clean up instances before solving. First, we describe how to build each subgraph N w = (N w , E w ). Then for each problem, (MC) and (PC), we sketch a method to eliminate O/D pairs which will never be covered.

To create N w we only consider useful nodes and edges from N . For each O/D pair w, we eliminate all the nodes i ∈ N w that do not belong to any path from w s to w t shorter than u w . Then, we define E w as the set of edges in E incident to the non eliminated nodes. Finally, the set A w is defined with the duplicated directed version of edges in E w with the exception of arcs in the form (i, w s ) and (w t , i). We describe this procedure in Algorithm 1.

We assume that there are no nodes or edges which cannot be built because their construction cost is higher than the budget.

Algorithm 1 Preprocessing I for w ∈ W do

N w = N for i ∈ N do compute the shortest path for the O/D pairs (w s , i) and (i, w t )
if the sum of the lengths of these paths is greater than u w then

N w = N w \ {i} E w = E w \ δ(i) end if end for A w = {(i, j) ∈ A : {i, j} ∈ E w , j = w s , i = w t } end for return {N w = (N w , E w ), A w } w∈W
Next, we focus on (MC). We can eliminate an O/D pair w that is too expensive to be covered.

That means, the O/D pair w is deleted from W if there is no path between w s and w t satisfying: i. its building cost is less than C max ; or ii. its length is less than u w .

This can be checked by solving a shortest path problem with resource constraints that can be done in a pseudo-polynomial time. [START_REF] Desrochers | An algorithm for the shortest path problem with resource constraints[END_REF] shows how to adapt Bellman-Ford algorithm to solve it. However, for the size of graphs that we are considering, we solve it as a feasibility problem.

For each w, we consider the feasibility problem associated to constraints (2.2) (2.3), (2.4), (2.5), (2.6) and (2.7), with z w fixed to 1. If this problem is infeasible, then the O/D pair w is deleted from W . Otherwise, there exists a feasible path denoted by Path w . We denote by ( E w , N w ) the subgraph induced by Path w .

Polyhedral properties

Both formulations (MC) and (PC) involve flow variables f w a whose number can be huge when the number of O/D pairs is large. To circumvent this drawback we use a Benders decomposition approach for solving (MC) and (PC).

In this subsection, we present properties of the two formulations that allow us to apply such a decomposition in an efficient way. The first proposition shows that we can relax the integrality constraints on the flow variables f w a . Let (MC-R) and (PC-R) denote the formulations (MC) and (PC) in which constraints (2.8) are replaced by nonegativity constraints, i.e.

f w a ≥ 0, w ∈ W, a ∈ A. (2.13)
We denote the set of feasible points to a formulation F by F(F ). Further, let Q be a set of points

(x, z) ∈ R q × R p .
Then the projection of Q onto the x-space, denoted P roj x Q, is the set of points given by P roj

x Q = {x ∈ R q : (x, z) ∈ Q for some z ∈ R p }. Proposition 1. P roj x,y,z F(M C) = P roj x,y,z F(M C -R) and P roj x,y,z F(P C) = P roj x,y,z F(P C - R).
Proof. We provide the proof for (MC), the other one being identical.

First, F(M C) ⊆ F(M C -R) implies P roj x,y,z F(M C) ⊆ P roj x,y,z F(M C -R).
Second, let (x, y, z) be a point belonging to P roj x,y,z F(M C -R). For every O/D pair w ∈ W such that z w = 0, f w = 0. In the case where z w = 1, there exists a flow f w a ≥ 0 satisfying (2.4) and (2.5) that can be decomposed into a convex combination of flows on paths from w s to w t and cycles.

Given that the flow f w a also satisfies (2.6), then a flow of value 1 on one of the paths in the convex combination must satisfy this constraint. Hence by taking f w a equal to 1 for the arcs belonging to this path and to 0 otherwise, we show that (x, y, z) also belongs to P roj x,y,z F(M C).

Based on Proposition 1, we propose a Benders decomposition where variables f w a are projected out from the model and replaced by Benders feasibility cuts. As we will see in Section 3.3, we also consider the Benders facet-defining cuts proposed in [START_REF] Conforti | Facet" separation with one linear program[END_REF]. To apply this technique it is necessary to get an interior point of the convex hull of P roj x,y,z F(M C -R) (resp. P roj x,y,z F(P C -R)). The following property give us an algorithmic tool to apply this technique to (MC).

Proposition 2. After pre-processing, the convex hull of P roj x,y,z

F(M C -R) is full-dimensional.
Proof. To prove the result, we exhibit |N | + |E| + |W | + 1 affinely independent feasible points:

• The 0 vector is feasible.

• For each i ∈ N , the points:

y i = 1, y i = 0, i ∈ N \ {i}, x e = 0, e ∈ E, z w = 0, w ∈ W
• For each e = {i, j} ∈ E, the points:

y k = 1, k ∈ e, y k = 0, k ∈ N \ {i, j}, x e = 1, x e = 0, e ∈ E \ {e}, z w = 0, w ∈ W
• For each w ∈ W , the points:

y i = 1, i ∈ N w , y i = 0, i ∈ N \ N w , x e = 1, e ∈ E w , x e = 0, e ∈ E \ E w , z w = 1, z w = 0, w ∈ W \ {w}
Clearly these points are feasible and affinely independent. Thus the polytope is full-dimensional.

The proof of Proposition 2 gives us a way to compute an interior point of the convex hull of

P roj x,y,z (F(M C -R)). The average of these |N | + |E| + |W | + 1 points is indeed such an interior point.
This is not the case for (PC) as we show in Example 1.

Example 1. Consider the instance of (PC) given by the data presented in Table 2 and Figure 1. We take the case where half of the population at least must be covered, that is β = 0.5. In order to satisfy the trip coverage constraint (2.10), the O/D pair W = (1, 4) must be covered. That is z (1,4) = 1 is an implicit equality. Furthermore, the only path with a length is less than or equal to u (1,4) = 15 is composed of edges {1,2} and {2,4}. Hence, x {1,2} , x {2,4} , y 1 , y 2 and y 4 must take value 1. In consequence, the polytope associated to (PC) is not full-dimensional.

Origin Destination u w g w We can compute the dimension of the convex hull of P roj x,y,z F(P C -R) in an algorithmic fashion.

We find feasible affinely independent points and at the same time we detect O/D pairs which must be covered in any feasible solution. Due to the latter, there are a subset of nodes and a subset of edges that have to be built in any feasible solution. This means that there is a subset of design variables y i , i ∈ N , x e , e ∈ E and mode choice variables z w , w ∈ W that must take value 1. At the opposite to (MC), a solution to (PC) with all variables set to 0 is not feasible. However, the solution obtained by serving all O/D pairs and building all nodes and edges is feasible. Therefore, we start with a solution with all variables in x, y, z set to 1 and we check, one by one, if it is feasible to set them to 0. By setting one variable x e or y i to 0, it may become impossible to cover some O/D pair w. In this case, we say that edge e and node i is essential for w. To simplify the notation, we introduce the binary parameters θ w e and θ w i taking value 1 if edge e (respectively node i) is essential for w. These new points are stored in a set L. Each time the algorithm finds a variable that cannot be set to 0, we store it in sets N , Ē, W , respectively. At the end of the algorithm, the dimension of the convex hull of P roj x,y,z

F(P C -R) is dim(P x,y,z ) = |N | + |E| + |W | -(| N | + | Ē| + | W |).
This procedure is depicted in Algorithm 2.

Algorithm 2 allows : i) to set some binary variables equal to 1, decreasing the problem size; and ii) to compute a relative interior point of the convex hull of P roj x,y,z F(P C -R), necessary for the facetdefining cuts, as explained below in Section 3.3. The relative interior point is given by the average of the points in set L.

Example 1 cont. Regarding the previous example and following Algorithm 2, the O/D pair (1, 4) must be covered, z (1,4) = 1. Due to that, as its shortest path in the networks (N (1,4) , E (1,4) \ {{1, 2}}) and (N (1,4) , E (1,4) \ {{2, 4}}) is greater than u (1,4) = 15, variables x {1,2} , x {2,4} , y 1 , y 2 , y 4 are set to 1.

Finally, the dimension of this polyhedron is dim(P x,y,z ) = 4 + 4 + 3 -(3 + 2 + 1) = 5.

The relative interior point computed is:

x {1,2} = 1, x {2,4} = 1, x {1,3} = 5 6 , x {3,4} = 2 3 , y 1 = 1, y 2 = 1, y 3 = 5 6 , y 4 = 5 6 , z (1,4) = 1, z (2,4) = 5 6 , z (3,4) = 1 2 .

Setting an initial solution

We determine an initial feasible solution for (MC) and (PC) with a simple greedy heuristic in which we sequentially select O/D pairs with best ratio demand over building cost. More precisely,

given the potential network N = (N, E), we compute for each O/D pair w the ratio r w =

g w C(Pathw)
, where C(Path w ) is the cost of a feasible path for w. We order these ratios decreasingly. We use this initial order in the heuristic both for (MC) and (PC). For (MC) the method proceeds as follows. It starts with an empty list of nodes and edges built, an empty list of O/D pairs covered, and a total cost set to 0. For each O/D pair w, in decreasing order of r w , the heuristic tries to build Path w considering edges and nodes that are already built. If the additional cost plus the current cost is less than the budget C max , nodes and edges in Path w are built and the O/D pair w is covered (i.e. z w = 1). The Algorithm 2 Computing the dimension of the polytope of (PC)

Initialization: Set N = ∅, Ē = ∅, W = ∅ and L = ∅ Add to set L: (y i = 1, i ∈ N, x e = 1, e ∈ E, z w = 1, w ∈ W ) . for w ∈ W do if w∈W \{w } g w ≥ β Z total then Add to set L: y i = 1, i ∈ N, x e = 1, e ∈ E, z w = 0, z w = 1, w ∈ W \ {w } . else W = W ∪ {w } for e = {i, j} ∈ E do
Compute shortest path between w s and w t in the graph (N w , E w \ {e}).

if the length of the shortest path is greater than u w or There is no path between w s and

w t then Ē = Ē ∪ {e} and N = N ∪ {i, j} end if end for end if end for for e ∈ E \ Ē do Add to set L: (y i = 1, i ∈ N, x e = 1, e ∈ E \ {e }, x e = 0, z w = 1 -θ we , w ∈ W ) .
end for

for i ∈ N \ N do
Add to set L:

y i = 0, y i = 1, i ∈ N \ {i }, x e = 0, i ∈ e, x e = 1, i / ∈ e, z w = 1 -θ wi , w ∈ W .
end for

dim(P x,y,z ) = |N | + |E| + |W | -(| N | + | Ē| + | W |).
return N , Ē, W , L and dim(conv(P x,y,z )).

total cost, the lists of nodes and edges built are updated. Otherwise we proceed with the next O/D pair. At the end of the algorithm we have an initial feasible solution.

To get an initial solution for (PC) we start with a list of all the O/D pairs covered and the amount of population covered equal to Z total . For each O/D pair w, in decreasing order of r w , the algorithm checks if by deleting the O/D pair w from the list, the coverage constraint (2.10) is satisfied. If so, the O/D pair w is deleted from the list and the amount of population covered is updated. Finally, the algorithm builds the union of the subgraphs ( N w , E w ) induced by Path w for all the O/D pairs covered.

Pseudo-codes for both routines are provided in Appendix A. In Section 4, we will show the efficiency of adding this initial solution at the beginning of the branch-and-benders-cut procedure.

Benders Implementations

In the following, we describe different Benders implementations for (MC) and (PC) obtained by projecting out variables f w a . These implementations are used as sub-routines in a branch-and-Benderscut scheme. This scheme cuts infeasible solutions along the branch-and-bound tree. Depending on the implementation, solutions can be cut whenever an integer solution is found or at any node in the tree. In the case of (MC), the master problem that we solve is:

(M -MC) max The master problem for (PC), named (M -PC), is stated analogously.

In Section 3.1, we discuss the standard Benders cuts obtained by dualizing the respective feasibility subproblem. Then, in Section 3.2 we propose a way of generating normalized subproblems, we named them normalized Benders cuts. In Section 3.3, we apply facet-defining cuts in order to get stronger cuts, as it is proposed in [START_REF] Conforti | Facet" separation with one linear program[END_REF]. Finally, we discuss an implementation where at the beginning cut-set inequalities are added to enhance the link between z and x, and then Benders cuts are added.

LP feasibility cuts

Since the structure of the model allows it, we consider a feasibility subproblem made of constraints (2.4), (2.5), (2.6) and (2.13) for each commodity w ∈ W , denoted by (SP) w . As it is clear from the context, we remove the index w from the notation. The dual of each feasibility subproblem can be expressed as:

(DSP)

w max α,σ,υ z α w s - e∈E x e σ e -u z υ (3.2) s.t. α i -α j -σ e -d a υ ≤ 0, a = (i, j) ∈ A : e = {i, j} (3.3) σ e , υ ≥ 0, e ∈ E (3.4)
where vector α is related to constraints (2.4), σ is the dual variable vector corresponding to the set of constraints (2.5) and υ is the dual variable of constraint (2.6). Since one of constraints in (2.4) is linearly dependent, we set α w t = 0. Given a solution of the master problem (x, y, z), there are two possible outcomes for (SP) w :

1. (SP) w is infeasible and (DSP) w is unbounded. Then, there exists an increasing direction (α, σ, υ)

with positive cost. In this case, the solution (x, y, z) is cut by:

(α w s -u υ) z - e∈E σ e x e ≤ 0 (3.5)
2. (SP) w is feasible and consequently, (DSP) w has an optimal objective value equal to zero. In this case, no cut is added.

Normalized Benders cuts

The overall branch-and-Benders-cut performance heavily relies on how the cuts are implemented.

It is known that feasibility cuts may have poor performance due to the lack of ability of selecting a good extreme ray (see for example [START_REF] Fischetti | A note on the selection of benders' cuts[END_REF]; [START_REF] Ljubić | Exact approaches to the single-source network loading problem[END_REF]). However, normalization techniques are known to be efficient to overcome this drawback [START_REF] Magnanti | Accelerating benders decomposition: Algorithmic enhancement and model selection criteria[END_REF]; [START_REF] Balas | Lift-and-project for mixed 0-1 programming: recent progress[END_REF][START_REF] Balas | A precise correspondence between lift-and-project cuts, simple disjunctive cuts, and mixed integer gomory cuts for 0-1 programming[END_REF]. The main idea is to transform extreme rays in extreme points of a suitable polytope. In this section we study three ways to normalize the dual subproblem described above.

First, we note that the feasibility subproblem can be reformulated as a min cost flow problem in N w with capacities x and arc costs d a .

(NSP) w min

x,y,z,f a∈A

d a f a (3.6) s.t (2.4), (2.5), (2.13)
The associated dual subproblem is:

(DNSP) w max α,σ z α w s - e∈E σ e x e (3.7) s.t α i -α j -σ e ≤ d a , a = (i, j) ∈ A : e = {i, j} (3.8) σ e ≥ 0, e ∈ E (3.9)
Whenever z > 0, the primal subproblem (NSP) w may be infeasible. Subproblems (NSP) w are no longer feasibility problems, although some of their respective dual forms can be unbounded. As splitting demand constraint has to be satisfied there are two kind of cuts to add:

1. (NSP) w is infeasible and (DNSP) w is unbounded. In this case, the solution (x, y, z) is cut by the constraint:

α w s z - e∈E σ e x e ≤ 0 (3.10)
2. (NSP) w is feasible and (DNSP) w has optimal solution. Consequently, if their solutions (α, σ)

and (x, y, z) satisfy that α w s z -e∈E σ e x e > u z then, the following cut is added

(α w s -u) z - e∈E σ e x e ≤ 0. (3.11)
We refer to this implementation BD Norm1.

In this situation, there still exists dual subproblems (DNSP) w with extreme rays. We refer to BD Norm2 as second dual normalization obtained by adding the dual constraint α w s = u + 1. In this case, every extreme ray of (SP) w correspond to one of the extreme points of (NSP) w . A cut is added whenever the optimal dual objective function is positive. This cut is in the form:

z - e∈E σ e x e ≤ 0 (3.12)
We finally tested a third dual normalization, BD Norm3, by adding constraints We tested the three dual normalizations described above for (MC) using randomly generated networks with 10, 20 and 40 nodes, as described in Subsection 4.1. Table 3 shows average values obtained for solution time in seconds and number of cuts needed for this experiment. The only one that seems competitive is BD Norm1. We observed that cut coefficients generated with BD Norm1 are mainly 0's or 1's. In the case of BD Norm2 and BD Norm3 we observe that coefficients generated are larger than the ones generated by BD Norm1, so they may induce numerical instability. This situation is similar for the case of (PC).

σ e ≤

Facet-defining Benders cuts

Here we describe how to generate Benders cuts for (MC) based on the ideas exposed in [START_REF] Conforti | Facet" separation with one linear program[END_REF]. The procedure for (PC) is the same. Given an interior point or core point, named as (x in , y in , z in ), of the convex hull of feasible solutions and a solution of the LP relaxation of the current restricted master problem, an exterior point, named (x out , y out , z out ), a cut that defines a facet or an improper face of the polyhedron defined by the LP relaxation of P roj x,y,z F(M C) is generated.

We denote the difference x outx in by ∆x. Analogously we define ∆y and ∆z. The idea is to find the furthest point from the core point, feasible to the LP-relaxation of P roj x,y,z F(M C) and lying on the segment line between the core point and the exterior point. This point is of the form (x sep , y sep , z sep ) = (x out , y out , z out ) -λ(∆x, ∆y, ∆z). The problem of generating cuts like this reads as follows:

(SP -CW)

w min f,λ λ (3.14) s.t. a∈δ + w (i) f a - a∈δ - w (i) f a =      z out -λ ∆z, if i = w s , 0, otherwise,
(3.15)

f a + f â ≤ x out e
λ ∆x e , e = {i, j} ∈ E : a = (i, j), â = (j, i), (3.16)

a∈A d a f a ≤ u z out -u ∆z λ, (3.17) 0 ≤ λ ≤ 1, (3.18) f a ≥ 0, a ∈ A (3.19)
In order to get the Benders feasibility cut we solve its associated dual:

(DSP -CW) w max α,σ,υ z out α w s - e∈E x out e σ e -u z out υ (3.20) s.t. ∆z α w s - e∈E ∆x e σ e -u ∆z υ ≤ 1, (3.21) α i -α j -σ e -d a υ ≤ 0, a = (i, j) ∈ A : e = {i, j}, σ e , υ ≥ 0, e ∈ E
Given that (SP-CW) w is always feasible (λ = 1 is feasible) and that its optimal value is lower bounded by 0, then, both (SP -CW) w and (DSP-CW) w have always finite optimal solutions. Whenever the optimal value of λ is 0, (x out , y out , z out ) is feasible. Cuts are added if the optimal value of (DSP-CW) w is strictly greater than 0. The new cut has the same form as in (3.5). Note that this problem can be seen as a dual normalized version of (SP) w with the dual constraint (3.21).

Core points for both formulations can be obtained by computing the average of the points described in the proof of Proposition 2 for (MC) and the average of the points in list L obtained by applying Algorithm 2.

Cut-set inequalities

By projecting out variable vector f, information regarding the link between vectors x and z is lost.

Cut-set inequalities represent the information lost regarding the connectivity for the O/D pair w in the solution given by design variable vector x. Let (S, S C ) a (w s , w t )-partition of N w for a fixed O/D pair w. That is (S, S C ) satisfies: i. w s ∈ S; ii. w t ∈ S C , with S C = N \ S its complement. A cut-set inequalities is defined as

z w ≤ {i,j}∈E w : i∈S, j∈S C
x {i,j} , p ∈ P, (S, S C ) a (w s , w t )-partition of N w (3.22)

This type of constraints have been studied in several articles, for instance [START_REF] Barahona | Network design using cut inequalities[END_REF]; [START_REF] Koster | Extended cutset inequalities for the network power consumption problem[END_REF]; [START_REF] Costa | Benders, metric and cutset inequalities for multicommodity capacitated network design[END_REF]. Note that is easy to see that cut-set inequalities belong to the LP-based Benders family. Let (S, S C ) be a (w s , w t )-partition in the graph N w for w ∈ W . We take the following dual solution:

• α i = 1 if i ∈ S; α i = 0 if i ∈ S C .
• σ e = 1 if e = {i, j} ∈ E w , i ∈ S, j ∈ S C ; σ e = 0, otherwise.

• υ = 0.

Note that this solution is feasible to (DSP) w and it induces a cut as in (3.22). In order to improve computational performance, we test two approaches to include these inequalities:

1. We implement a modification of the Benders callback algorithm with the following idea. First, for each w ∈ W , using the solution vector (x, y) from the master, the algorithm generates a network (N w , E w ) with capacity 1 for each edge built. Then, a Depth-First Search (DFS) algorithm is applied to obtain the connected component containing w s . If the connected component does not contain w t , a cut in the form (3.22) is added. Otherwise, we generate a Benders cut as before.

This routine is depicted in Algorithm 3.

We tested this implementation with subproblems in the form (DSP-CW) w . We observe that by using Algorithm 3 with CW the convergence is slower and we generate more cuts. This might Algorithm 3 Callback implementation with cut-set inequalities.

Require: (x e , e ∈ E, z w , w ∈ W ) from the master vector solution (x, y, z).

for w ∈ W do Build graph (N w (x), E w (x)) induced by the solution vector x from the master.

Compute the connected component S in (N w (x), E w (x)) containing w s .

if w t is included in S then
Add the cut z w ≤ {i,j}∈E w : i∈S, j∈S C

x {i,j} else Solve the corresponding subproblem ((DSP) w , (DNSP) w , (DSP-CW) w ) and add cut if it is necessary.

end if end for return Cut.

be due to the fact that these cuts do not include information about the length of the path in the graph, but only information regarding the existence of the path. These preliminary results 2. We add to the Master Problem the cut-set inequalities at the origin and at the destination of each O/D pair w ∈ W at the beginning of the algorithm. In particular, this valid inequalities has the form:

are shown in

         z w ≤ e∈δ(w s ) x e , z w ≤ e∈δ(w t ) x e , (3.23)
This means that for each O/D pair to be covered, there should exist at least one edge incident to its origin and one edge incident to its destination, i.e. each O/D pair should have at least one arc going out of its origin and another one coming in its destination.

Computational Results

In this section, we compare the performance of the different families of Benders cuts presented in Section 3 using the branch-and-Benders-cut algorithm (noted B&BC).

Data sets: benchmark networks and random instances

We divide the tested instances into two groups: benchmarks instances and random instances. Our benchmarks instances are composed by the Sevilla García-Archilla et al. ( 2013) and Sioux networks "Hellman (Accessed June 16th, 2020).

The Sevilla instance is composed partially by the real data given by the authors of [START_REF] García-Archilla | Grasp algorithms for the robust railway network design problem[END_REF]. From this data, we have used the topology of the underlying network, cost and distance vector for each arc and the demand matrix. This network is composed of 49 nodes and 119 edges.

Originally, the set of O/D pairs W was formed by all possible ones (49 • 48 = 2352). However, some entries in the demand matrix of this instance are 0 and we exclude them from the analysis. We consider as private utility u twice the shortest path in the underlying network. Each node cost is generated according to an uniform distribution U(2000, 4000). The available budget has been fixed as 30% of the cost of building the whole underlying network and the minimum demand to be covered as

β = 0.5.
For the Sioux instance, the topology of the network is described by 24 nodes and 38 edges. Set W is also formed by all possible O/D pairs (38 • 37 = 1406). With respect to the parameters, they have been chosen in the same manner as for random instances.

We generate our random instances as follows. We consider planar networks with a set of n nodes, with n ∈ {10, 20, 40, 60}. Nodes are placed in a grid of n square cells, each one of 10 units side. For each cell, a point is randomly generated close to the center of the cell. For each setting of nodes we consider a planar graph with its maximum number of edges, deleting each edge with probability 0.3.

We replicated this procedure 10 times for each n, so that the number of nodes is the same while that the number of edges may vary. Therefore, there are 40 different underlying networks. We name these instances as N 10, N 20, N 40 and N 60. We provide the average cycle availability, connectivity and density for random instances networks in Table 5. A couple of them are depicted in Figure 2. Construction costs b i , i ∈ N , are randomly generated according to an uniform distribution U(7, 13).

Then, each node costs 10 monetary units on average. Construction cost of each edge e ∈ E, c e , is set to its Euclidean length. These two parameters are rounded to integer numbers. It means that building the links cost 1 monetary unit per length unit. We set C max equal to 50% of the cost of building the whole underlying network considered. We denote this total cost as T C, so C max = 0.5 T C.

To build set W , we randomly pick each possible O/D pair of nodes with probability 0.5. In consequence, this set has n(n-1) 2 pairs on average. Parameter u w is set to 2 times the length of the shortest path between w s and w t . We name this shortest distance SP w . Finally, the demand g w for each O/D pair w is randomly generated according to the uniform distribution U(10, 300).

Algorithmic setting

Our preliminary experiments have shown that including cuts only in integer nodes is more efficient than including them in nodes with fractional solutions. Thus, in our experiments we only separate integer solutions unless we specify the opposite. We study the different implementations of B&BC proposed in Sections Sections 3.1, 3.2 and 3.3. We name them with the following nomenclature:

• BD Trd: B&BC algorithm using the feasibility subproblems structure (DSP) w , and its corresponding feasibility cuts (3.5).

BD Norm and BD CW are able to solve all the instances in N40 within the time limit for both problems.

For problems in set N60, we also have better gaps after one hour in comparison with AUTO BD. This difference is significant in (MC) with a reduction of more than 5%, but it is less significant for (PC).

For N60 we compare the performance by setting an initial feasible solution IS and adding cuts at the root node RNC. We perform this experiment by computing the optimality gap after one hour limit.

First we note that we obtain worse solutions by adding RNC in both problems with all the algorithms tested. For (MC) we observe that adding an initial solution is only profitable for BD CW+CS, obtaining in average a 3.5% better optimality gap. In the following, we analyze the performance of algorithms BD Norm+CS BD CW+CS when changing parameters C max , β and u in the corresponding models. In Tables 14 and15, we report average solution times and cuts needed to obtain optimal solutions for N40 for different values of these parameters.

The instances are grouped by the three different increasing values of the available budget C max (Table 14.a) or β ( For (PC), Table 15.a shows that for β = 0.7 both algorithms take less time to solve at optimality in comparison with β = 0.3 and β = 0.5. BD CW+CS is 5 minutes faster in average than BD Norm+CS with β = 0.5. For β = 0.3 the result is the opposite, BD Norm+CS is 100 seconds faster in average than BD CW+CS. By varying u, we observe that the less the difference between public and private mode distances in the underlying network, the longer it will take to get optimality.

Performance of algorithms in benchmark instances

We start analyzing the Sevilla instance. Tables 17 and18 show some results for this instance solved with BD CW+CS. Based on this case, figures in Tables 17 and18 show the solution graphs for different parameter values. Points not connected in these graphs refer to those nodes that have not been built. The O/D pairs involving some of these nodes are thus not covered. They have been drawn

to represent these not covered areas. Data corresponding to each case are collected at the bottom of its figure, in which v(ILP) refer to objective value. For model (MC), parameter cost represents the cost of the network built, and, for (PC), Z makes reference to the demand covered. For (MC) we see that the solution times vary from 21 seconds to 2243 seconds depending on the parameters C max and u. For (PC) these times are in the range of 353 seconds to 1358 seconds. We observe that smaller values of C max carry bigger solution times as in random instances. We also observe that, as opposite to random instances, higher values of β are translated in larger solution times. Besides, in this instance, for both models, the shorter is the parameter u the larger are the solution times.

Furthermore, we compare the performance of the GRASP algorithm from García-Archilla et al.

(2013) and our implementation BD CW+CS. We implemented the GRASP algorithm to run 5 times and return the best solution. Table 16 shows solution times, best value for GRASP (Best Value), the optimality gap, and the optimal value computed with BD CW+CS. On the one hand, we observed that the more time BD CW+CS takes to compute the optimal solution the larger is the gap of the solution returned by GRASP. This happens for smaller values of the budget C max and utility u. On the other hand, for problems where GRASP obtains small optimality gap, BD CW+CS is more efficient to compute the optimal solution. In other words, since GRASP is a constructive algorithm, it is not competitive for instances whose optimal solution captures most of the demand. Our exact method is able to get the best quality solution, with a certificate of optimality in reasonable times. Given that network design problems are strategic decisions, having the best quality decision is often more important than the computational times. However, having efficient exact methods as the proposed in this article, allows decision makers to perform sensitivity analysis with optimality guarantees in reasonable times.

Conclusions

In a lower bound in the total population covered. We state integer programming formulations that are stronger than existing ones for both problems. We provide some polyhedral properties of these formulations useful from the algorithmic point of view. We develop exact methods based on Benders decomposition. We also discuss some preprocessing routines to scale-up the instances solved. This preprocessing techniques play a key role in order to get information about the instances and to get better algorithmic performance. Our computational results show that the techniques developed in this article allow to obtain better solutions in less time than the techniques in the existing literature.

Figure 1 :

 1 Figure 1: Graph of Example 1.

Figure 2 :

 2 Figure 2: Example of underlying networks with |N |=20 and |N |=40.

Finally

  we discuss the results for the Sioux instance. They have also been obtained by using BD CW+CS. We observe for (MC), as in the Sevilla network, that smaller values of C max and u the larger solution times. The same is true varying β in (PC), but not for u. It takes less time if the difference between both modes of transport is smaller or larger than 2 SP . These results are summarized inTables B.19 and B.20 in Appendix B.

C

  max = 0.3 T C, u = 2 SP T C = 295297826, Z total = 64011 t = 313.07 , cuts = 7149 cost = 88939.76, v(ILP) = 62294 C max = 0.2 T C, u = 2 SP C max = 0.4 T C, u = 2 SP t = 1036.11, cuts = 8428 t = 21.36, cuts = 2259 cost = 53362.53, v(ILP) = 52274 cost = 124353.84, v(ILP) = 64011 C max = 0.3 T C, u = 1.5 SP C max = 0.3 T C, u = 3 SP t= 2243.83, cuts = 8537 t = 113.88, cuts = 6400 cost = 88800.91, v(ILP) = 59958 cost = 88879.76, v(ILP) = 62670

Table 1 :

 1 Comparing the performance of the two different types of mode choice and capacity constraints for (MC) within a time limit of 1 hour. The majority of N40 instances were not solved to optimality, then the average gap is shown.

		t	LP gap	t	LP gap
	N10	0.17	43.21	0.26	96.43
	N20	5.78	56.33	228.22	106.71
		gap	LP gap	gap	LP gap
	N40	11.74	68.15	54.85	137.13

Network

Formulation using (2.5)-(2.6) Formulation using (2.11)-(2.12)

Table 3 :

 3 Comparing the performance of the three dual normalization within a time limit of 1 hour for (MC). N10,

	1,	e ∈ E,	(3.13)

Table 4 :

 4 Table 4, which shows average values obtained for solution times in seconds and the number of cuts needed. Comparing the performance of the Algorithm 3. N10, N20 and N40 refer to networks with 10, 20 and 40 nodes respectively.

	Network	BD CW	Algorithm 3+BD CW
		t	cuts	t	cuts
	N10	0.23	48	0.15	46
	N20	2.47	411	2.53	500
	N40	619.31 3486 722.02	3554

Table 5 :

 5 Cycle availability, connectivity and density parameters for the underlying networks in random instances.

  The impact of adding an initial solution for (PC) is significant for BD Trd+CS, BD Norm+CS and BD CW+CS obtaining in average solutions with a gap around 4% smaller. This improvement is not significant for BD Auto. In summary, in the set of instances N60 we have that the best algorithms are BD CW+CS+IS for (MC) and BD CW+CS+IS and BD Norm+CS+IS for (PC).

	Network	CPLEX	Auto BD	BD Trd	BD Norm	BD CW
		t		t	cuts	t		cuts	t	cuts	t	cuts
	N10	0.18		0.43	27	0.25	92	0.24	91	0.19	94
	N20	6.77		4.51	273	3.89	620	3.18	590	3.34	641
	N40	1646.93 617.85 1967 1095.25 3990 541.03 3677 457.81 4137
	Network	Auto BD+CS		BD Trd+CS	BD Norm+CS	BD CW +CS
			t	cuts	t	cuts	t	cuts	t	cuts
	N10	0.32	12	0.21	49	0.28	52	0.23	54
	N20	3.94	178	2.29	382	2.50	383	1.85	416
	N40	484.95 1248 637.49 2378 575.87 2530 272.39 3186

Table 6 :

 6 Comparing the performance of the three algorithms for (MC).

		CPLEX Auto BD BD Trd BD Norm BD CW
	without CS	3	10	9	8	8
	+CS	-	10	10	10	10

Table 7 :

 7 Instances N40 solved for (MC) within a time limit of 1 hour.

Table 8 :

 8 Computing gaps to solve N60 (MC) instances comparing the performance of three families of Benders cuts.

	AUTO BD+CS+IS	BD Trd+CS+IS	BD Norm+CS+IS	BD CW+CS+IS
	gap	cuts	gap	cuts	gap	cuts	gap	cuts
	without RNC 32.90	4987	27.23	9038	26.94	9469	22.27	11151
	+RNC	-	37.88	8054	37.92	8230	33.58	10834

Table 9 :

 9 Computing gaps to solve N60 (MC) instances comparing the performance of three families of Benders cuts.

	Network	CPLEX	Auto BD	BD Trd	BD Norm	BD CW
		t		t	cuts	t	cuts	t	cuts	t	cuts
	N10	0.1765	0.29	16	0.24	92	0.28	89	0.20	91
	N20	6.7281	4.87	305	3.55	607	4.68	681	2.15	606
	N40	2153.1458 504.06 1752 657.59 4470 514.42 4246 837.41 4412
	Network	Auto BD+CS	BD Trd+CS	BD Norm+CS	BD CW+CS
			t	cuts	t	cuts	t	cuts	t	cuts
		N10	0.28	11	0.16	56	0.20	57	0.145	54
		N20	4.12	213	3.11	497	3.43	495	2.070	461
		N40	439.23 1527 261.74 3528 323.21 3583 197.55 3949

Table 10 :

 10 Comparing the performance of the three algorithms for (PC).

		CPLEX Auto BD BD Trd BD Norm BD CW
	without CS	3	9	8	8	8
	+CS	-	10	10	10	10

Table 11 :

 11 Instances N40 solved for (MC) within a time limit of 1 hour.

		Auto BD	BD Trd	BD Norm	BD CW
		gap	cuts	gap	cuts	gap	cuts	gap	cuts
	without CS 20.49 7009 20.40 14784 21.41 15501 19.93 15116
	+CS	15.92 5109 14.89 12354 14.09 11687 14.50 11744

Table 12 :

 12 Computing gaps to solve N60 (PC) instances comparing the performance of three families of Benders cuts

	AUTO BD+CS+IS	BD Trd+CS+IS	BD Norm+CS+IS	BD CW+CS+IS
	gap	cuts	gap	cuts	gap	cuts	gap	cuts
	without RNC 15.86	4372	11.06	8961	10.47	8490	10.44	9683
	+RNC	-	20.93	10971	21.28	11449	19.94	11053

Table 13 :

 13 Computing gaps to solve N60 (PC) instances comparing the performance of three families of Benders cuts.bigger the values of C max the shorter is the average solution time. Table14.b. shows that the larger is the parameter u the shorter is the solution time for BD Norm+CS. This behavior seems to be different if we are using BD CW+CS, which takes less time if the difference between both types of transport is smaller or larger than 2 SP .

	C max	BD Norm+CS t cuts	BD CW+CS t cuts	u	BD Norm+CS t cuts	BD CW+CS t cuts
	0.3 T C 1053.56 1580 873.58 2017	1.5 SP 802.05 2792 495.84 3041
	0.5 T C	622.45 2634 375.30 3358	2 SP	622.46 2634 375.30 3358
	0.7 T C	151.24 3970 177.90 5035	3 SP	591.02 2674 490.28 3173
		a.			b.

Table 14 :

 14 Sensitivity analysis for (MC) with |N | = 40.

Table 15 :

 15 Sensitivity analysis for (PC) with |N | = 40.

Table 16 :

 16 this article, we have studied two variants of the Network Design Problem: Maximal Covering Network Design Problem where we maximize the population covered under a budget constraint; and Partial Set Covering Network Design Problem where the total building cost is minimized satisfying Sensitivity analysis for GRASP algorithm García-Archilla et al. (2013) with Sevilla instance.

	C max	u		GRASP		BD CW+CS
			t	Best Value gap	t	v(ILP)
	0.2 T C		110.829	48629	6.97 1036.11	52274
	0.3 T C	2 SP	260.220	59828	3.96 313.07	62294
	0.4 T C		396.226	63546	0.72	21.36	64011
	0.3 T C	1.5 SP 267.275 3 SP 225.312	55778 62049	6.97 2243.83 0.99 113.88	59958 62670

Table 17 :

 17 Sensitivity analysis Sevilla Network with (MC).

Acknowledgments

Víctor Bucarey and Martine Labbé have been partially supported by the Fonds de la Recherche Scientifique -FNRS under Grant(s) no PDR T0098.18. Natividad González-Blanco and Juan A. Mesa are partially supported by Ministerio de Economía y Competitividad (Spain)/FEDER under grant MTM2015-67706-P.

• BD Norm: B&BC algorithm using the normalized subproblems structure (DNSP) w , and its corresponding cuts (3.10) and (3.11).

• BD CW: B&BC algorithm using the subproblems structure (DSP -CW) w , and feasibility cuts (3.5).

We compare our algorithms with the direct use of CPLEX, and the automatic benders proposed by CPLEX, noted by AUTO BD. CPLEX provides different implementations of its automatic Benders depending on the information that the user provides to the solver: i. CPLEX attempts to decompose the model strictly according to the decomposition provided by the user; ii. CPLEX decomposes the model by using this information as a hint and then refines the decomposition whenever possible; iii.

CPLEX automatically decomposes the model, ignoring any information supplied by the user. We tested these three possible settings, and the only competitive was the first one.

Furthermore we tested the following features:

• CS: If we include cut-set inequalities at each origin and destination as in (3.23).

• IS: If we provide an initial solution to the solver.

• RNC: If we add Benders cuts at the root node.

Performance of algorithms on random instances

All the experiments were performed with one hour time limit. Tables in this section show average values obtained for solution times in seconds, percent relative gaps, and number of cuts needed.

We consider in the average only the instances solved at optimality for all the algorithms. First, we compare the performance of CPLEX for formulations (MC) and (PC) and the different B&BC implementations. All the algorithms are able to solve at optimality instances N10 and N20 in less than 7 seconds for (MC) and (PC). Instances in set N40 were not all solved at optimality neither for (MC) nor (PC) (see first block rows in Table 7 andTable 11). For (MC) the fastest algorithm was BD CW in sets N10, N20 and N40 for the instances solved at optimality. This does not happen for (PC), since we can observe that AUTO BD is slightly faster. AUTO BD solved more instances for both, (MC) and (PC). This trend is confirmed in (MC) for instances in set N60 were the optimality gap obtained after one hour is smaller in AUTO BD as it is shown in Table 8. However, for (PC) the gap after one hour for BD CW is slightly better than the other methods in this family (see Table 12). Now we refer to the effect in the performance by including CS. In general, solution times decrease as the amount of cuts required as well when CS is considered. The only exception occurs with BD Norm for (MC) where the computing time is slightly larger. Despite the fact that BD CW generates a larger amount of cuts, this is the most efficient method in terms of computational times for (MC) and (PC) in families N10, N20 and N40. In contrast to the case without CS, by adding these initial cuts BD Trd,

Appendix A. Pseudo-code for initial feasible solutions

In this section we provide the pseudo-codes to get an initial feasible solution for (MC) and (PC) described in Section 2.4. We denote as N s , E s and W s the set of indices of design and mode choice variables set to 1 at the end of each algorithm.

Algorithm 4 Initial Feasible Solution for (MC)

end if end for

x e = 1 for e ∈ E s , 0 otherwise. :

end if end for

x e = 1 if e ∈ w∈Ws E w , 0 otherwise.