Hub Location with Protection Under Interhub Link Failures

Submitted to INFORMS Journal on Computing manuscript (Please, provide the manuscript number!) Authors are encouraged to submit new papers to INFORMS journals by means of a style file template, which includes the journal title. However, use of a template does not certify that the paper has been accepted for publication in the named journal. INFORMS journal templates are for the exclusive purpose of submitting to an INFORMS journal and should not be used to distribute the papers in print or online or to submit the papers to another publication. Hub Location with Protection under Inter-Hub Link Failures V´ıctor Blanco Institute of Mathematics (IMAG), Universidad de Granada, and Dpt. Quant. Methods for Economics & Business, Universidad de Granada [email protected] Elena Fern´andez Dpt. Statistics & Operations Research, Universidad de C´adiz [email protected] Yolanda Hinojosa Dpt. Applied Economics I, Universidad de Sevilla, and Institute of Mathematics (IMUS), Universidad de Sevilla yhino[email protected] This paper introduces the Hub Location Problem under Inter-Hub Link Failures, a hub location problem, in which activated inter-hub links may fail with a given probability. Two different optimization models are studied, which construct hub backbone networks protected under inter-hub link disruptions by imposing that for each commodity an additional routing path exists besides its original routing path. Both models consider the minimization of the fixed costs of the activated hubs and inter-hub links plus the expected value of the routing costs of the original and alternative paths. The first model builds explicitly the alternative routing paths, whereas the second model guarantees that for each commodity at least one alternative path exists using a large set of connectivity constraints, although the alternative paths are not built explicitly. The results of extensive computational testing allow to analyze the performance of the two proposed models and to evaluate the extra cost required to design a robust backbone network under inter-hub link failures. The obtained results support the validity of the proposal. Key words : Hub Location; Integer Programming; Robust Network Design; Disruptions 1. Introduction Hub location (HL) lies in the intersection of location analysis and network design, and produces challenging optimization problems with multiple applications, mostly in the fields of distribution and logistics (parcel delivery, air transportation, etc.) and telecommunications (Farahani et al. 2013). The increasing attention they have received in the last decades is thus not surprising (see e.g. Campbell and O’Kelly 2012, Erdo˘gan et al. 2022, Oliveira et al. 2022, Zheng et al. 2022). The 1 V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 2Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) interested reader is referred to the book chapter by Contreras and O’Kelly (2019) for a recent survey on HL, including terminology, modeling assumptions, etc. One of the current trends in HL is the search of models suitable for dealing with different sources of uncertainty (Alumur et al. 2012). Whereas some models in the literature consider uncertainty in demand (Contreras et al. 2011, Zetina et al. 2017), other models are concerned with the robustness of solution hub networks, by associating uncertainty with the possibility (probability) of disruption of the elements involved in the solution networks and looking for solutions that are robust under disruptions. Most hub networks are sensitive to failures in their links, being the impact on users particularly harmful in some cases. This is precisely the focus of this work, where we introduce the HL Problem under Inter-Hub Link Failures (HLPIHLF), a HL problem, in which activated inter-hub links may fail with a given probability. Our study can be very useful in typical HL applications, in which the total failure of a hub is highly unlikely, whereas partial failures occur only affecting some of the links incident with the hubs (certain air connections, train lines, etc.) We point out that by protecting inter-hub links under failure, we also partially protect hub nodes under failures. Examples of potential applications of the models that we study include the management of airlines and airport industries (Campbell et al. 2005), in which breakdowns in certain flight connections may occur, and passengers are directly affected, as well as rapid delivery packing systems (C¸etiner et al. 2010), where users pay for fast services and failures in the hub network may cause large delays. For dealing with the HLPIHLF we propose two alternative models, which guarantee that solution networks are protected under disruption of inter-hub links, in the sense that for each origin/destination demand pair (commodity) at least one alternative (backup) routing path exists. Both models can be seen as two-stage stochastic programming models, in which the a priori solution is determined by the strategic decisions associated with the selection of activated hub nodes and inter-hub links, together with an original plan given by a set of feasible routing paths, one for each commodity, containing exactly one inter-hub arc (directed link). The recourse action determines a backup plan, given by a set of alternative routing paths for the commodities, which can be used in case the inter-hub link of the original plan fails. The main difference between the two models is how backup paths are enforced. The first model imposes that the alternative routing path of each commodity contains exactly one inter-hub arc (as in the original plan) and builds it explicitly. The second model is more flexible, in the sense that it allows for arbitrarily large sequences of inter-hub arcs to be used in the alternative routing paths, although such paths are not built explicitly. This is achieved with a set of exponentially many (on the number of nodes of the network) constraints, by imposing that the cut-set of the V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 3 backbone network contains at least λlinks, for a given integer value of λ≥2. We study some properties of both models and propose a Mixed Integer Linear Programming (MILP) formulation for each case. For the second model, since it has exponentially many constraints, we also propose a branch-and-cut solution algorithm. In both cases we consider fixed costs for both activated hubs and activated inter-hub links. Thus, we do not fix the number of inter-hub links to activate, allowing for incomplete backbone networks. As usual in HL, the routing costs apply a discount factor αto inter-hub arcs. Extensive computational experiments have been carried out on a large set of benchmark instances based on the well-known CAB (O’Kelly 1987), AP (Ernst and Krishnamoorthy 1996), and TR (Tan and Kara 2007) datasets, for varying settings of the failure probabilities and other cost parameters. The obtained results are summarized and analyzed, comparing the computational performance of each of the models and the effect of the different parameters. Managerial insights are derived from the analysis of the characteristics of the solutions produced by each of the models and their comparison. In particular, we analyze the distribution of the costs among the different elements considered (routing costs and fixed costs for hubs and links), the number of activated hub and links, and the density of the obtained backbone networks. Finally, we have carried out an empirical analysis of the different models in terms of their efficiency and robustness. Specifically, for each of the proposed models, we analyze their a posteriori capability to re-route the commodities under multiple failure scenarios. The obtained results assess the validity of the proposal. The remainder of this paper is structured as follows. In Section 2 we review the related literature and in Section 3 we introduce the notation that we will use and formally define the HLPIHLF from a general perspective. In that section we also specify and motivate the modeling assumptions that we make. Most of them are usual in classical models in the HL literature. This allows us to better focus on the novel aspects of our proposal. Section 4 is devoted to the first HLPIHLF model that we analyze, in which it is assumed that the alternative paths for the commodities contain exactly one inter-hub arc. We study some of its properties and propose a MILP formulation for it. The model in which we impose that the backbone network is λ-connected is detailed in Section 5 where we also present a MILP formulation for it. Section 6 describes the computational experiments we have carried out and summarizes the obtained results. Some managerial insights from the analysis of the structure of the solution networks produced by each model are also derived in this section. Finally, Section 7 describes the empirical analysis that has been performed considering multiple failure scenarios. The paper closes in Section 8 with some conclusions. 2. Related literature Several works have studied models, in which it is assumed that activated hub nodes may fail (totally or partially) with a certain probability. An et al. (2015) propose a model, in which two backup hub V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 4Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) nodes are determined for each commodity. Rostami et al. (2019) assume that a finite set of hub breakdown scenarios is known and provide a two-stage formulation for the single allocation HL problem with possible hub breakdown. Cui et al. (2010) provide a stochastic model to determine a subset of the activated hub nodes of given size through which each commodity can be routed, in such a way that if the cheapest route fails, the commodity can be routed through the second cheapest, and so on, or through an emergency facility. The authors provide a MILP formulation for the problem as well as an approximate Lagrangean relaxation scheme for its solution based on the ideas in Snyder and Daskin (2005) for the p-median problem. The planar version of this model is also analyzed there. Kim and O’Kelly (2009) propose the reliable p-hub location problem (PHMR) and the p-hub mandatory dispersion (PHMD). In the PHMR the goal is to determine the location of pnodes based on the level of reliability to maximize the flows sent through the selected hubs. The PHMD imposes a certain minimum separation between the pselected hub nodes, and maximizes the reliability of the network. MILP formulations and heuristic approaches are provided for the problems, both in the single and the multiple-allocation framework. The reliability of hub backbone networks has been also studied in Zeng et al. (2010), Korani and Eydi (2021) and Li et al. (2022). Parvaresh et al. (2013) consider the multiple allocation p-hub median problem under intentional disruptions, in which the goal is to identify the optimal strategy for the location of phub nodes by minimizing the expected transportation cost, taking into account the worst-case situation when an interdictor attacks the backbone network. A bilevel mixed integer formulation is provided as well as a simulated annealing heuristic for its solution. The problem of designing robust networks under link failures has also been studied in the literature under different settings. Aneja et al. (2001) study the single-commodity maximum flow problem when link failures occur by means of the maximum residual flow problem, whose goal is to determine the maximal flow, in which the largest link flow is as small as possible. This problem is closely related to the network interdiction problem that consists of determining a certain number of links whose removal from the network minimizes the maximum amount of flow that one can send through the network (see e.g. Altner et al. 2010, Cormican et al. 1998, Royset and Wood 2007, Wood 1993). Ma et al. (2016) propose the Conditional Value-at-Risk Constrained Minimum Spanning k-Core Problem where the possibility of link disruptions is prevented when trying to construct minimum-cost subgraphs of a network with a minimum number of incident links at each node. Andreas and Smith (2008) study shortest path problems under link failures by imposing that the probability that all links that successfully operate in at least one path be greater than certain threshold value. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 5 We are only aware of one work (Mohammadi et al. 2019), dealing with potential failure of interhub links, which is the specific focus of our work. Nevertheless, in contrast to Mohammadi et al. (2019), who consider a p-hub single allocation pattern, and complete hub networks with no set-up costs for inter-hub links, we adopt a multiple allocation policy on a non-complete backbone network where the number of hubs is not a fixed parameter, and assume fixed costs for inter-hub links. Furthermore, Mohammadi et al. (2019) apply meta-heuristics to solve the models whereas we solve them to proven optimality. 3. Notation and definition of the problem Consider a graph N= (V,E), where the node set V={1,2,...,n}represents a given set of users and the edge set Ethe existing undirected connections or links between pairs of users. We assume that Nis a complete graph and that Econtains loops, i.e. for all i∈V, edge {i, i} ∈ E. We further assume that potential locations for hubs are placed at the nodes of the graph and the set of potential locations coincides with V. For each potential location k∈V, we denote by fkthe fixed cost for activating a hub at node k. Any pair of hub nodes can be connected by means of an inter-hub edge, provided that both end nodes kand lare activated as hub nodes as well. The set Ewill be referred to as the set of potential inter-hub edges or just as set of potential hub edges. Activated hub edges incur fixed costs; let hkl ≥0 be the fixed cost for activating hub edge {k, l} ∈ E. A set of activated hubs will be denoted by H⊆V, and a set of activated hub edges for Hby EH⊆E[H], where E[H] is the set of edges with both end nodes in H. Note that the assumption that Nis a complete graph implies no loss of generality, since (i) arbitrarily large fixed costs can be associated with nodes that are not potential hubs; and, (ii) arbitrarily large activation costs can be associated with non-existing hub edges. Activating a hub edge allows to send flows through it in either direction. Let A={(i,j)∪(j,i) : {i,j} ∈ E, i,j ∈V}be the (directed) arc set. Arcs in the form (i,i), i ∈V, are allowed and will also be called loops. We will use AH⊆Ato denote the set of inter-hub arcs induced by EH. Service demand is given by a set of commodities defined over pairs of users, indexed in a set R. Let D={(or,dr,wr) : r∈R}denote the set of commodities, where the triplet (or,dr,wr) indicates that an amount of flow wr≥0 must be routed from origin or∈Vto destination dr∈V. The origin/destination pair associated with a given commodity will also be referred to as its OD pair. Commodities must be routed via paths of the form π= (or,k1,...,ks,dr) with ki∈H, 1 ≤i≤s. Similarly to most HL problems, a routing path π= (or,k1,...,ks,dr) is feasible if (1) it includes at least two hub nodes, i.e. s≥2 (being possible that ki=ki+1 for some i= 1,···,s −1), and (2) the underlying edges of all traversed arcs other than the access and delivery arcs, (or,k1) and (ks,dr), respectively, are activated inter-hub edges, i.e., {ki, ki+1} ∈ EH, 1 ≤i≤s−1. Note that this V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 6Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) implies that in any feasible path all intermediate nodes are activated as hubs as well, i.e. ki∈H, 1≤i≤s. In the following, the set of feasible paths for a given commodity r∈Rwill be denoted by Π(H,EH)(r). Routing flows through the arcs of a hub-and-spoke network incurs different types of costs. These costs, which may depend on the type of arc, account for transportation costs as well as for some additional collection/handling/distribution costs at the end nodes of the arcs. As usual in the literature, we assume that transportation costs of flows routed through inter-hub arcs are subjected to a discount factor 0 ≤α≤1. In this work we will denote by cij ≥0 the unit routing cost through inter-hub arc (i,j)∈AH, which includes discounted transportation costs and handling costs, and we will denote by ¯cij ≥0 the unit routing cost for an access or a delivery arc (i,j)∈A\AH, which may also incorporate different discounted access or delivery costs. With the above notation, the routing cost of commodity r∈Rthrough a feasible path π= (or,k1,...,ks,dr)∈Π(H,EH)(r) is: Cr π=wr cork1+ s−1 X i=1 ckiki+1 +cksdr!,(1) where the first and last addends correspond to the access and delivery arcs, respectively, and the intermediate ones are the service costs through the backbone network (H,EH). Broadly speaking, under the above assumptions, the goal of a HL problem is to decide the location of the hub nodes Hand to select a suitable subset of hub edges EH, to optimally route the commodities through the backbone network (H,EH) induced by the activated hub nodes and hub edges, so as to minimize the sum of the overall fixed costs for activating hub nodes and hub edges, plus the commodities routing costs. With the above notation, this problem can be stated as: min H⊆V,EH⊆E[H]X k∈H fk+X e∈EH he+X r∈R min π∈Π(H,EH)(r)Cr π.(HLP) In the reminder of this paper we make the following assumptions: •All our models assume multiple allocation of OD nodes to open hubs. With the multiple allocation policy, in case two commodities share the same origin, each of them is allowed to route its demand using a different access arc. Similarly, two commodities with the same destination may be routed using different delivery arcs. This strategy provides greater flexibility to hub backbone networks than the single allocation pattern, in which OD nodes are assigned to a single hub node such that all the demand originated at (or arriving to) such a node is routed initially (or finally) through the same access (or delivery) arc. The multiple allocation policy is assumed in several so-called fundamental HL models (see, e.g., Contreras and O’Kelly 2019), and is commonly applied V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 7 by delivery companies that allocate origins to hub centers based on the destination of the deliveries, or by airlines that route passengers departing from the same city through different intermediate airports (hubs) based, for instance, on the continent of their final destination. •Loops (self-arcs) of the form (k,k) are considered as hub arcs. This is a slightly more general setting than in most HL problems studied in the literature, which only consider proper arcs (k,l) with k=las potential hub arcs. In this work, if a loop (k,k) is used in a routing path, it is required not only that kis activated as a hub node, but also that the loop {k,k}is activated as a hub edge as well. The rationale for this assumption is that loops allow one to model the different operations that are required when flows enter a hub network, as the preparation of the machinery, workforce, etc. As far as we know, loops have never been explicitly considered in the literature of HL, although loops are present in other types of network design problems (see, e.g., Wei et al. 2021, Huang and Tian 2017). •Activated hub edges in EHmay fail. We are given the failure probabilities for each potential inter-hub edge, pkl,{k, l} ∈ E. When k=lfailure of edge earises not only when the edge {k, l} can no longer be used, but also when, for any reason, the collection and distribution services at any end node of edge {k,l}cannot be carried out. Thus pkl represents the probability that any of these events happen. In case edge eis a loop, i.e., e={k,k}, with k∈H, then pkk represents the probability that the handling process carried out when kis used as the unique intermediate hub fails. •In case hub edge {k,l} ∈ EHfails, then the inter-hub arcs (k,l),(l,k)∈AHcan no longer be used for routing the commodities. In order to protect solution networks from failure we adopt a policy that does not alter the strategic decisions on the activated hubs and inter-hub edges and focuses solely on the operational decisions concerning the re-routing of affected commodities. Accordingly, we impose that, for each commodity, the backbone network (H,EH) contains a routing path and, in addition, a substitute path connecting orand dr, for each r∈R. Such a substitute path will be referred to as backup or alternative path whereas, the routing path used in case the edges do not fail will be called original or initial path. The allocation pattern may be different in the original and backup paths, so they may use different access/delivery arcs. •The original routing path of each commodity r= (or,dr,wr) contains exactly one inter-hub arc, which can be a loop. Most HL problems studied in the literature (see, e.g., Contreras and O’Kelly 2019), which do not consider loops explicitly, restrict the set of potential paths for routing the commodities to those using at most one hub arc. When loops are also considered as potential hub arcs as we do, potential routing paths are of the form (or,k,m, dr), where kand mare activated hubs, and (k,m) is an inter-hub arc, which reduces to a loop when k=m. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 8Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) Thus, for Hand EHgiven, the set of potential original paths for routing commodity r∈Ris given by Π(H,EH)(r) = {π= (or,k,l,dr) : k,l ∈H, {k, l} ∈ EH}. In such a case the routing cost of commodity rthrough path π= (or, k,l,dr) reduces to Cr kl := Cr π=wr(cork+ckl +cldr). •For each edge e={k,l} ∈ E, let Xkl denote the random variable modeling whether efails. We assume that the random variable Xkl follows a Bernoulli distribution with probability pkl, for each {k, l} ∈ E. Observe that when edges may fail with a given probability, the costs of feasible routing paths are also random variables, which will be denoted by Cr,r∈R. Furthermore, the probability distribution of Cr,r∈R, is dictated by the failure probability distribution of the involved inter-hub edges. In particular, when π0= (or,k, l, dr) is the original routing path of a given commodity r∈R, the expected routing cost of commodity rcan be calculated as: E[Cr|π0] = (1 −pkl)Cr kl +pklCr π0.(2) where π0∈Π(H,EH)(r)\{π0}is the backup path in case {k,l}fails. •The objective is to minimize the sum of the fixed costs of the activated hubs and inter-hub edges, plus the expected routing costs. In the following we deal with the problem of finding hub backbone networks protected against inter-hub edge failures under the above assumptions. That is, the HLPIHLF can be stated as: min H⊆V,EH⊆E[H]X k∈H fk+X e∈EH he+X r∈R min π0∈Π(H,EH)(r)E[Cr|π0].(HLPIHLF) Indeed, multiple choices fall within the above generic framework which differ from each other in how the alternative routing paths are obtained. In the following sections we propose two models for determining such backup paths, based on different assumptions, and provide mathematical programming formulations for each of them. 4. HLPIHLF with single inter-hub arc backup paths The HLPIHLF we address in this section enforces that the alternative paths for routing the commodities have the same structure as the original ones. That is, we assume that backup paths contain exactly one inter-hub arc (possibly a loop). This avoids having to use many transshipment points in case of failure. This model will be referred to as (HLPIHLF-1BP). Given a commodity r∈R, the backbone network (H, EH) and the original routing path π0= (or,k,l, dr), we assume that the backup path is in the form π0= (or,¯ k,¯ l,dr), with {k,l} ={¯ k,¯ l}. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 9 Thus, the expected routing cost of commodity ris: E[Cr|π0] = (1 −pkl)Cr kl +pklCr ¯ k¯ l.(3) Figure 1 illustrates the different situations that may arise in case a hub edge fails. Figure 1a shows a backbone network with four hub nodes (k,l,mand q) and five hub edges, two of them corresponding to loops, {m,m}and {l,l}, and the other ones corresponding to edges {k, q},{k,l} and {q,m}. The figure also depicts the origin (or) and destination (dr) of a given commodity r∈R and a possible path for this commodity through hub arc (k,l). We assume that this is the original path for the commodity r∈R. Access/distribution arcs are depicted as dashed lines and hub edges as solid lines. Figures 1b and 1c show different backup paths with a single inter-hub arc for commodity r∈R in case the inter-hub link {k,l}of the original path fails. In Figure 1b, the backup path uses the inter-hub arc (q,m) to re-route the commodity, whereas in Figure 1c the backup path uses the loop arc (l,l). Next we develop a mathematical programming formulation for the above problem, first introducing the decision variables associated with the design decisions on the elements of the network that are activated, hubs nodes and hub edges: zk=1 if a hub is opened at the potential hub node k, 0 otherwise for k∈V, ykl =1 if hub edge {k, l}is activated, 0 otherwise for {k, l} ∈ E. The formulation uses two additional sets of variables, which represent the original and alternative routing path for each commodity, respectively. In particular, for r∈Rand (k,l)∈A: xr kl =1 if the original routing path for commodity ris (or, k,l,dr), 0 otherwise, ¯xr kl =1 if the alternative path for commodity r∈Ris (or,k, l,dr), 0 otherwise. With these sets of decision variables the expected routing cost of commodity r∈Rcan be expressed as: X (k,l)∈A xr klhCr kl(1 −pkl) + pkl X (k′,l′)∈A\{(k,l)∪(l,k)} Cr k′l′¯xr k′l′i,(4) where the two addends in each term of the above expression correspond to the expected routing cost of the original and backup plan of commodity r, respectively. Observe that both terms only apply if, in the original plan, the commodity is routed through the inter-hub arc corresponding V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 16 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) Below we develop a MILP formulation for the HLPIHLF-λ, which incorporates λ-connectivity by means of the family of Constraints (14) introduced above. The formulation uses the same z,y, and xvariables as before. Still, since backup paths are no longer made explicit, variables ¯xused in formulation HLPIHLF-1BP of the previous section are no longer needed. As explained, the y variables will be used to impose the λ-connectivity condition, which will be stated by means of an exponential set of constraints. Given that backup routes are no longer made explicit, we no longer have closed expressions for their expected routing costs and we must estimate their values. Let us denote by ¯ Cr kl an estimation of the backup routing cost of commodity rwhen the hub arc (k,l) of its original routing path fails. The resulting formulation for the HLPIHLF-λis: min X k∈V fkzk+X {k,l}∈E hklykl+X r∈RX (k,l)∈A(1 −pkl)Cr kl +pkl ¯ Cr klxr kl s.t. X (k,l)∈A xr kl = 1 ∀r∈R(2.1) xr kl +xr lk ≤ykl ∀r∈R,{k,l} ∈ E, k =l(2.2) xr kk ≤ykk ∀r∈R,k ∈V(2.3) ykl ≤zk∀{k,l} ∈ E(2.4) ykl ≤zl∀{k,l} ∈ E,k =l(2.5) y(δ(S)) + ykk ≥λ(zk+zl−1) ∀S⊂V, |S| ≥ 2, k ∈S, l /∈S(2.6) y(δ(k)) + ykk ≥λzk∀k∈V(2.7) xr kl ∈ {0,1} ∀r∈R,(k,l)∈A(2.8) zk∈ {0,1} ∀k∈V(2.9) ykl ∈ {0,1} ∀{k,l} ∈ E, (2.10) where Constraints (2.1)-(2.5) are similar to (1.1)-(1.6) but referring to the original path only, and (2.6) and (2.7) are the λ-connectivity constraints described above. Note that once the hub backbone network (H, EH) is obtained by solving the above problem, one can explicitly compute a backup path for a commodity r∈R, whose original path is (or, kr, lr, dr), by solving (in polynomial time) a shortest path problem from source orto destination dron the graph Gr= (Vr,Er) with nodes Vr={or,dr} ∪ Hand edges Er=n{or,h}:h∈Ho∪Eh∪n(h,dr) : h∈Ho\{kr,lr}. We next show that, in some cases, Constraints (2.6) can be reinforced. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 17 Proposition 5.1 Let S⊂Vbe a nonempty subset of nodes with |S| ≤ λ−1. Then, the associated Constraints (2.6) are dominated by the following set of inequalities, which are also valid for (HLPIHLF−λ): y(δ(S)) + ykk ≥λzk,∀k∈S. (2.11) Proof. We first see that (2.11) are valid for (HLPIHLF−λ). Taking into account that λconnectivity with λ≥2 implies that any feasible solution has at least λopen hubs and that |S| ≤ λ−1, when k∈Sis activated as a hub node (i.e. zk= 1), there will be at least one more open hub ¯ l /∈S(i.e. z¯ l= 1). That is, when zk= 1 there will be at least one active constraint in the set (2.6) with right-hand-side value λ(zk+z¯ l−1) = λzk. When kis not activated as a hub node (i.e. zk= 0) none of Constraints (2.6), nor (2.11) will be active. Furthermore, (2.11) dominate (2.6), since zl≤1 implies that λzk≥λ(zk+zl−1). Therefore, the result follows. □ Note that when Sis a singleton, i.e., S={k}, the set of Constraints (2.11) reduces precisely to (2.7). 5.1. Incorporation of λ-cutset constraints: a branch-and-cut approach As already mentioned, in (HLPIHLF−λ), the size of the family of Constraints (2.6) is exponential in the number of potential hub nodes, n. It is thus not possible to solve the formulation directly with some off-the-shelf solver, even for medium size instances. In this section we present an exact branch-and-cut algorithm for this formulation in which, as usual, the family of constraints of exponential size (2.6) is initially relaxed. The strategy that we describe below is embedded within an enumeration tree and it is applied not only at the root node but also at all explored nodes. Our separation procedure is an adaptation of the separation procedure for classical connectivity constraints (Padberg and Gr¨otschel 1985), and follows the same vein of those applied to more general connectivity inequalities in node and arc routing problems (see, e.g., Belenguer and Benavent 1998, Ar´aoz et al. 2009, Rodr´ıguez-Pereira et al. 2019, for further details). The initial formulation includes all Constraints (2.1)-(2.5), and (2.7). Furthermore, all integrality conditions are relaxed. Let (x,y,z) be the solution to the current linear programming relaxation and let G(y) = (V(y),E(y)) denote its associated support graph where E(y) consists of all the edges of Esuch that ¯ykl >0 and V(y) the set of end nodes of the edges of E(y). Each edge (k,l)∈E(y) is associated with a capacity ¯ykl. The separation for inequalities (2.6) is to find S⊂V, and k∈S,l∈ S, with y(δ(S)) < λ(zk+zl−1)−ykk or to prove that no such inequality exists. Note that, when they exist, violated λ-connectivity Constraints (2.6) can be identified from a tree of min-cuts associated with G(y) relative to the capacities vector y,T(y). V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 18 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) Therefore, to solve the above separation problem we proceed as follows. For each min-cut δ(S) of T(y) of value y(δ(S)), we identify (k,l)∈arg max{λ(zk+zl−1) −ykk :k∈S,l ∈ S}. Then, if y(δ(S)) < λ (zk+zl−1) −ykk, the inequality (2.6) associated with Sand (k,l) is violated by y. In case |S| ≤ λ−1, according to Proposition 5.1, we incorporate cuts in the shape of (2.11) instead of (2.6). In our computational experiments, we use the procedure proposed by Gusfield (1993) to identify T(y). Such an algorithm computes V(y) max-flows in G(y), so its overall complexity is O(|V(y)| × |V(y)|3). 6. Computational Experience In this section we report the results of an extensive battery of computational tests, which have been carried out to analyze the performance of the two modeling approaches for obtaining robust hub networks protected under inter-hub failures, discussed in the previous sections. For the experiments, we have used a large set of benchmark instances based on the wellknown CAB (O’Kelly (1987)), AP (Ernst and Krishnamoorthy (1996)) and TR (Tan and Kara (2007)) datasets (taken from ORLIB http://people.brunel.ac.uk/~mastjjb/jeb/orlib/ and https://ie.bilkent.edu.tr/~bkara/dataset.php), for varying settings of the failure probabilities and other parameters as described below. All instances were solved with the Gurobi 9.1.1 optimizer, under a Windows 10 environment on an Intel(R) Core(TM) i7-6700K CPU @ 4.00 GHz 4.01 GHz processor and 32 GB of RAM. Default values were used for all parameters of Gurobi solver and a computing time limit of 7200 seconds was set. 6.1. Instance generation We have generated several instances based on the entire CAB, AP and TR datasets with a number of nodes (n) initially ranging in {10,15,20,25}for the instances based on the CAB and TR datasets and in {10,20,25}for the instances based on the AP dataset. Let c′ kl be the standard unit transportation costs provided in ORLIB for CAB and AP instances or the travel distances provided for the TR instances and let wbe the normalized vector of service demands, taken from the provided datasets. Table 1 summarizes the main characteristics of the testing instances and the selected parameters. A detailed description of how the different parameters have been generated can be found in Section Instance Generation of the Online Supplement. The files of the randomly generated probabilities, as well as the Python codes used for such a generation, are available in the Github repository https://github.com/vblancoOR/HLPLF. For each combination of parameters n, α, ρ, and each dataset (CAB, AP with type T and L hub nodes activation costs, and TR) five different instances have been generated for scenario of failure V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 19 Instances n λ α β CAB {10,15,20,25} {2,3,4} {0.2, 0.5, 0.8} {0.5,1} AP {10,20,25,40,50} TR {10,15,20,25,40,50} Costs Instances ckl fkhkl ¯ Cr kl CAB α(ak+c′ kl +dl) 100          100ckl/wkl maxw if k=l, 100 ckl/¯ w maxw if k=l. (1 + β)Cr kl AP Tand Ldata files TR Data file Data file ak=dk= min{minj=kc′ kj,minj=kc′ jk};¯ w:= Mean of w;maxw= max{cij wij :i,j ∈V, wij >0} Failure probabilities Random probabilities (RP) pkl ∼U[0,ρ], ρ ∈ {0.1,0.3} Clustered probabilities (CP) pkl ∈ {0.1,0.2,0.3} Same probability (SP) pkl =ρ, ρ ∈ {0.1,0.3} Table 1 Summary of instances and parameters. probabilities RP and one instance has been considered for scenario SP. Five instances have been also generated for scenario CP and each combination of n, α, and each dataset. Thus, (HLPIHLF-1BP), hereafter called M1, has been solved on a total of 714 instances. Concerning formulation (HLPIHLF-λ), we considered three different values for the parameter λ, namely λ= 2, λ= 3 and λ= 4 (we call the corresponding models M2 2, M2 3 and M2 4, respectively) and two values for β. Thus, (HLPIHLF-λ) has been solved on a total of 4284 instances. Additionally, for comparative purposes, we have solved 42 instances of the uncapacitated HL problem, in which no protection under failures is considered. This model will be referred to as M0. Finally, to test the scalability of our formulations, a second experiment was carried out on a set of larger instances (n∈ {40,50}) based on the AP and TR datasets considering only (HLPIHLF-λ), V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 20 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) CPUTime MIPGAP %Solved RP CP SP RP CP SP RP CP SP 0,1 0,3 0,1 0,3 0,1 0,3 0,1 0,3 0,1 0,3 0,1 0,3 nαData Av Max Av Max Av Max Av Max Av Max 10 0,2 APT4 7 15 1 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 APL7 93 13 2 3 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 CAB 143 5198 TL 1 1 0,0 0,0 0,5 1,1 1,1 1,9 0,0 0,0 0,0 0,0 100 40 0 100 100 TR 4 9 44 0 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 0,5 APT6 21 14 1 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 APL10 112 11 2 3 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 CAB 62 5391 1864 1 1 0,0 0,0 1,0 3,5 0,0 0,0 0,0 0,0 0,0 0,0 100 40 100 100 100 TR 5 11 63 0 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 0,8 APT7 13 15 1 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 APL9 99 12 1 3 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 CAB 25 TL 727 1 1 0,0 0,0 0,6 1,2 0,0 0,0 0,0 0,0 0,0 0,0 100 0 100 100 100 TR 5 19 53 0 1 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 15 0,2 CAB TL TL TL 2 7 0,6 1,0 8,3 9,6 9,4 10,8 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 51 71 287 4 6 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 0,5 CAB 4577 TL TL 5 8 0,4 0,7 9,8 12,5 7,6 8,4 0,0 0,0 0,0 0,0 40 0 0 100 100 TR 52 89 619 3 5 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 0,8 CAB 511 TL TL 4 6 0,0 0,0 8,8 13,3 3,7 5,5 0,0 0,0 0,0 0,0 100 0 0 100 100 TR 55 138 412 3 4 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 20 0,2 APT1812 TL 6336 43 30 0,0 0,0 5,9 15,8 1,4 2,0 0,0 0,0 0,0 0,0 100 0 20 100 100 APLTL TL TL 3443 2080 14,8 15,4 20,5 21,8 17,8 19,1 0,0 0,0 0,0 0,0 0 0 0 100 100 CAB TL TL TL 13 121 3,1 3,7 14,4 16,5 15,4 17,9 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 360 1402 4131 32 32 0,0 0,0 0,0 0,0 1,1 5,3 0,0 0,0 0,0 0,0 100 100 80 100 100 0,5 APT1520 TL 6611 53 34 0,0 0,0 2,9 7,4 1,3 3,4 0,0 0,0 0,0 0,0 100 0 20 100 100 APLTL TL TL 4308 1257 13,1 14,5 20,3 21,9 17,2 19,2 0,0 0,0 0,0 0,0 0 0 0 100 100 CAB TL TL TL 20 51 3,0 3,8 15,7 19,5 15,5 19,8 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 363 1850 5497 21 36 0,0 0,0 0,0 0,0 3,6 12,0 0,0 0,0 0,0 0,0 100 100 60 100 100 0,8 APT1693 7137 5780 40 33 0,0 0,0 3,8 12,1 1,6 4,3 0,0 0,0 0,0 0,0 100 20 40 100 100 APLTL TL TL 2343 1809 12,6 13,6 19,0 20,6 15,6 16,5 0,0 0,0 0,0 0,0 0 0 0 100 100 CAB TL TL TL 12 27 2,8 3,9 16,4 19,3 16,2 18,3 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 422 1864 5625 19 33 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 0,0 100 100 100 100 100 25 0,2 APT4774 6835 4731 128 129 2,4 11,8 12,5 21,6 0,0 0,0 0,0 0,0 0,0 0,0 80 20 100 100 100 APLTL TL TL TL TL 16,6 17,9 22,8 24,1 19,8 20,4 11,5 11,5 14,2 14,2 0 0 0 0 0 CAB OoM TL TL 62 3208 4,3 5,4 19,5 21,7 19,2 23,2 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 2332 TL TL 119 169 0,0 0,0 12,0 13,2 17,3 22,4 0,0 0,0 0,0 0,0 100 0 0 100 100 0,5 APT6389 TL 4789 129 179 5,0 9,5 12,6 19,3 0,0 0,0 0,0 0,0 0,0 0,0 40 0 80 100 100 APLTL TL TL TL TL 14,7 16,3 20,8 24,4 18,1 20,2 8,6 8,6 13,4 13,4 0 0 0 0 0 CAB OoM OoM TL 98 2403 4,7 6,7 23,2 24,2 20,4 23,3 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 1935 TL TL 82 128 0,0 0,0 11,5 13,3 13,3 14,3 0,0 0,0 0,0 0,0 100 0 0 100 100 0,8 APT6465 6789 3490 125 316 5,1 8,7 10,3 20,6 0,0 0,0 0,0 0,0 0,0 0,0 40 20 100 100 100 APLTL TL TL 2927 TL 13,8 14,8 20,5 21,8 17,5 22,5 0,0 0,0 11,9 11,9 0 0 0 100 0 CAB OoM OoM TL 71 2740 7,1 9,2 24,0 25,1 18,9 20,6 0,0 0,0 0,0 0,0 0 0 0 100 100 TR 2483 TL TL 53 80 0,0 0,0 12,7 14,6 14,2 15,1 0,0 0,0 0,0 0,0 100 0 0 100 100 Table 2 Results for (HLPIHLF-1BP). which, as we will see, is the most promising formulation, for λ∈ {2,4}and β= 1. We have solved a total of 612 instances in this second study. Overall, 5652 instances have been solved. 6.2. Numerical results with (HLPIHLF-1BP) and (HLPIHLF-λ) In the following sections we present and compare the results obtained in our computational study with (HLPIHLF-1BP) and (HLPIHLF-λ). 6.2.1. Results with (HLPIHLF-1BP) The results produced by (HLPIHLF-1BP) are summarized in Table 2. In that table, “RP”, “CP” and “SP” stand for the scenarios with random failure probabilities (with ρ= 0.1 and ρ= 0.3), clustered failure probabilities and same failure probability (with ρ= 0.1 and ρ= 0.3), respectively. The values of n,αand “Data” indicate the number of V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 21 nodes in the network, the value for the discount factor applied to the routing cost through inter-hub arcs, and the dataset that has been used to obtain the costs and the flows, respectively. “APT” and “APL” refer to AP dataset using type T and type L fixed costs for the hub nodes, respectively. For scenarios RP and CP, the information contained in each row refers to average values over the five instances with the corresponding combination of parameters, whereas for scenario SP the values of the entries correspond to the unique instance with this combination of parameters. The maximum MIPGap values over the five instances are also shown under “Max” columns. The numerical results are grouped in three blocks of columns. Block “CPUTime” gives the computing times, in seconds, required to solve the instances, block “MIPGap” the percentage MIP gaps returned by Gurobi at termination, and block “%Solved” the percentage of instances solved to proven optimality within the time limit. An entry “TL” in the CPUTime block means that the time limit of 7200 seconds was reached in all five instances of the group. The “OoM” entry indicates that the flag “Out of memory” was the output of the solver in at least one of the instances in the row, and then, the remaining information of the row refers to average values (or maximum values when corresponding) over the solved instances only (even if none of these instances could be solved to proven optimality). We observe that there is no significant difference between the maximum and the average MIP gap values. Thus, in what follows we only report average values in our tables and figures. Table 2 shows a different performance of (HLPIHLF-1BP) among the instances corresponding to the different scenarios. Based on the computing times, MIPGaps, and percentage of solved instances, scenario SP produces the easiest instances, for all configurations of parameters, as expected. One can observe that, for instances with the same probability, all instances generated from the CAB, APT, and TR datasets with up to n= 25, as well as the instances generated from APLwith up to n= 20 have been optimally solved within the time limit. On the other hand, note that instances based on CAB dataset are more difficult to solve than instances based on TR and AP datasets. TR based instances are the easiest to solve: all instances with up to n= 15, 95% for n= 20 and 35% for n= 25 have been optimally solved within the time limit. Regarding AP based instances, all instances with n= 10 have been optimally solved within the time limit, although APTinstances consumed, in general, less computing time. The difference between APTand APLinstances becomes more evident for n > 10, since approximately 50% of the APTinstances were optimally solved whereas no APLinstance with random or clustered probability (scenarios RP and CP) was solved to proven optimality within the time limit. As mentioned before, CAB instances are the most difficult ones. For n= 10, 30% of these instances could not be optimally solved solved within the time limit. This percentage increases up to 90% for n= 20 and up to 100% for n= 25. Additionally, for n= 25, the execution was stopped due to an Out of Memory flag for 30% of the CAB instances under a random probability (RP) scenario. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 22 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 6.2.2. Results with (HLPIHLF-λ)Using formulation (HLPIHLF-λ), all instances could be solved to proven optimality within the maximum computing time for all three considered values of λ∈ {2,3,4}. This means that (HLPIHLF-λ) is notably easier to solve than (HLPIHLF-1BP). This can be explained by its smaller number of decision variables. The difficulty of (HLPIHLF-λ) increases with the value of λ, as reflected by a decrease in its performance for higher values of this parameter. This could be expected, as (HLPIHLF-λ) becomes more restrictive as the value of λincreases. When λ= 2, the average computing time over all the instances is approximately 64 seconds, being two seconds for the CAB instances, 70 seconds for the APTinstances, 89 seconds for the APLinstances, and 102 for the TR instances. Unlike (HLPIHLF-1BP), CAB instances are computationally less demanding than AP and TR instances. This behavior was also observed for λ= 3 and λ= 4. For λ= 4 the average computing time over all the instances is approximately 218 seconds, being 30 seconds for the CAB instances, 168 seconds for the APTinstances, 225 seconds for the APLinstances, and 438 for the TR instances. We have observed that the value of αalso affects the performance of (HLPIHLF-λ), instances being more difficult for smaller αvalues, specially for the AP instances. We also noted that, unlike (HLPIHLF-1BP), with (HLPIHLF-λ) there seem to be no noticeably differences among scenarios. Detailed results on the computing times required in the implementation of (HLPIHLF-λ) with the tested instances can be found in the Online Supplement. Given the good performance of (HLPIHLF-λ) with instances up to 25 nodes, we carried out a second set of computational experiments, on larger instances (n∈ {40,50}) based on the AP and TR datasets. For this second set of experiments we only considered (HLPIHLF-λ) for λ∈ {2,4} and β= 1, but did not consider (HLPIHLF-1BP), since it most instances could not be optimally solved with it within the time limit already for n= 25. Table 3 summarizes the obtained results. We can observe that, for λ= 2, all the TR instances, 99% of the APTinstances, and 80% of the APLinstances have been solved to proven optimality within the time limit, whereas for λ= 4 the percentage of solved instances decreased to 70% for APT, 40% for APL, and 6% for TR. This shows that (HLPIHLF-λ) is able to solve larger instances with up to n= 50, even if instances become more challenging as the value of the parameter λincreases. 6.3. Managerial Insight In this section we derive some managerial insight from the results obtained in our first set of experiments, i.e., n≤25 when the instances were solved with both formulations, as well from the solutions of these instances for M0 (the uncapacitated HLP with no protection under failures). V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 23 CPUTime MIPGap %Solved RP SP RP SP RP SP nαData 0.1 0.3 CP 0.1 0.3 0.1 0.3 CP 0.1 0.3 0.1 0.3 CP 0.1 0.3 λ= 2 40 0.2 APT371 363 382 361 310 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL1608 1689 5418 1181 1059 0.00 0.00 4.43 0.00 0.00 100 100 60 100 100 TR 870 775 773 666 537 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 0.5 APT344 479 313 322 432 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL3579 3155 2487 2602 1351 0.00 3.16 0.00 0.00 0.00 100 80 100 100 100 TR 434 600 601 337 414 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 0.8 APT331 376 325 303 337 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL2067 1663 2126 1604 692 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 TR 499 489 430 512 390 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 50 0.2 APT3986 3470 5128 4147 3651 0.00 0.00 2.46 0.00 0.00 100 100 80 100 100 APL6427 5362 6472 3273 TL 13.79 6.56 12.40 0.00 15.64 20 60 40 100 0 TR 2877 2520 3385 2180 3817 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 0.5 APT2635 3789 4319 4823 1590 0.00 0.00 2.02 0.00 0.00 100 100 80 100 100 APL3819 3406 TL 2761 2131 0.00 3.45 17.09 0.00 0.00 100 80 0 100 100 TR 1296 3037 3264 1201 1349 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 0.8 APT2862 3732 2836 3736 2474 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL5872 5151 7183 6279 1172 9.93 2.99 14.02 0.00 0.00 40 80 20 100 100 TR 1010 1495 1356 852 1050 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 λ= 4 40 0.2 APT2919 2907 2300 2582 6060 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APLTL TL 5643 5427 5674 49.57 50.23 19.51 0.00 0.00 0 0 60 100 100 TR 6958 7044 6002 TL TL 13.64 9.16 5.89 7.14 10.61 20 20 60 0 0 0.5 APT2807 2335 2293 2186 2939 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL5589 4865 3875 5535 2960 10.03 10.03 0.00 0.00 0.00 80 80 100 100 100 TR TL 7030 7024 TL TL 20.15 14.00 12.62 22.24 13.99 0 20 20 0 0 0.8 APT2315 2400 2272 1640 1750 0.00 0.00 0.00 0.00 0.00 100 100 100 100 100 APL4959 5691 3512 6253 2767 19.83 9.12 0.00 0.00 0.00 60 80 100 100 100 TR TL 6991 7139 TL TL 17.73 13.56 12.25 17.51 14.20 0 20 20 0 0 50 0.2 APT7089 7198 6966 7134 6004 30.03 40.96 27.71 0.00 0.00 40 20 40 100 100 APLTL TL TL TL TL 49.99 49.23 50.31 48.24 46.73 0 0 0 0 0 TR TL TL TL TL TL 17.12 13.76 16.23 17.77 16.54 0 0 0 0 0 0.5 APT6482 6995 6711 5310 TL 20.67 39.79 25.85 0.00 49.76 60 20 40 100 0 APLTL TL TL TL TL 49.68 48.38 50.03 48.53 48.19 0 0 0 0 0 TR TL TL TL TL TL 24.59 16.95 17.05 22.29 20.88 0 0 0 0 0 0.8 APT6933 7081 6259 7202 7201 28.66 39.67 0.00 29.26 46.94 40 20 100 0 0 APLTL TL TL TL TL 48.91 47.81 47.93 48.24 46.94 0 0 0 0 0 TR TL TL TL TL TL 25.94 21.03 22.78 27.43 23.51 0 0 0 0 0 Table 3 Average results for (HLPIHLF-λ) for n≥40. We first analyze the contribution to the objective function value of the different types of costs: routing costs, fixed costs for activating hubs (Hubs Costs) and fixed costs for activating inter-hub edges (Links Costs). Since results are similar for APLand APTdata sets, for this analysis we differentiate, for each formulation, between datasets CAB, AP and TR, as well as among the three values of the αparameter. The obtained results indicate that the percent contribution of the fixed costs for activating interhub edges varies from 0.1% for CAB instances to 13% for TR instances in M0. As expected, the percent contribution of hub fixed costs depends on the value of the parameter α, but mainly on the dataset and on the model. For the instances based on the CAB dataset, the percent contribution of hub fixed costs varies from 20% with M0, M1, and M2 2 for α= 0.8, to 40% with M2 4 for α= 0.2. For the instances based on the AP dataset, the percent contribution of hub fixed costs varies from 45% with M0 and α∈ {0.5,0.8}to 75% with M2 4. Regarding the instances based on the TR dataset, the percent contribution of hub fixed costs varies from 21% with M0 and α∈ {0.5,0.8}to V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 24 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 69% with M2 4 for α= 0.2. The percent contribution of routing costs also depends on the value of the parameter α, on the dataset, and on the model. For the instances based on the CAB dataset, this percentage varies from 55% to 80%, with the highest values for M0, M1, and M2 2 and α= 0.8. For the instances based on the AP dataset, the percent contribution of routing costs varies from 15% to 55%, corresponding the highest values to M0 with α= 0.8. Finally, for the instances based on the TR dataset, the percent contribution of routing costs varies from 22% to 68%, corresponding again the highest values to M0 with α= 0.8. In the Online Supplement we provide the bar diagrams of the percent contribution of each type of cost for the different models. For analyzing the structure of the backbone networks produced by the different tested models we have studied the number of open hubs and activated inter-hub edges, together with two density indices, denoted by I1and I2.I1indicates the density of the backbone network including loops and is computed as the ratio between the number of activated inter-hub edges and the total number of edges in a complete graph with #Hnodes if loops were included: I1=2#Lk #H(#H+ 1). I2indicates the density of the backbone network when loops are not considered, and is computed as the ratio between the number of non-loop inter-hub edges activated and the number of edges in a complete graph with #Hnodes if loops were excluded: I2=2(#Lk −#Lp) #H(#H−1) . Note that values of the two indices range in [0,1], i.e., 0 ≤I1, I2≤1. The obtained results show that the number of open hubs mostly depends on the model and also on the dataset. For all the AP instances the number of open hubs is always one for M0, two for M1 and M2 2, three for M2 3, and four for M2 4. For the CAB instances this number ranges, in average, in [2, 5] for M0, in [2.18, 4.29] for M1, in [2, 4.76] for M2 2, in [3, 4.76] for M2 3, and in [4, 4.74] for M2 4. For the TR instances, the number of open hubs is always one for M0, two for M1, three for M2 2, four for M2 3, and five for M2 4. Concerning the number of activated inter-hub edges, our results indicate that this number is smaller for M0 than for M1 and that, for M1, this number is similar to that for M2 2 but smaller than that for M2 λfor λ > 2, since, as expected, this number increases with the value of λ. Additionally, with M0 most of the activated inter-hub edges are loops, specially with the AP and the TR instances. In fact, we have observed that M0 gives a density index I1close to 1 and a density index I2close to 0, for both the AP instances and TR instances. This indicates that most V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 25 of the activated inter-hub edges are loops. Instead, for the TR instances the optimal backbone networks produced by M1 have both I1=I2= 1. We finally point out that, as expected, the density of the backbone network for M2 λincreases with the value of the parameter λ. In the Online Supplement we provide with several tables and figures that support the above analysis. In particular, we report the average values for the number of hubs (#H), inter-hub edges including loops (#Lk), and loops (#Lp) activated in the optimal backbone networks for the different models; and depict bar diagrams of the values for the density indices I1and I2, for the different types of instances. 7. The price of robustness For assessing the robustness and reliability of the hub network models proposed in this paper, we evaluate the so-called price of robustness (see e.g. Bertsimas and Sim 2004), defined as the extra cost incurred to design a robust network. In our case, robustness translates into protecting the backbone network under inter-hub edge failures, which, essentially, is attained by incorporating additional inter-hub edges to the backbone network. We thus start our analysis by comparing the overall fixed cost for each of the proposed models M1, M2 2, M2 3, and M2 4, with that of the unprotected network obtained with M0. The interested reader is referred to the Online Supplement for detailed information about these costs. As expected, constructing networks that are robust under inter-hub edge failures has a significant impact in the design costs (fixed costs for activating hub nodes and inter-hub edges) of the network. We have observed that, in this type of instances, the value of αscarcely affects the design costs in any of the models. The percent deviation of the activated hubs and inter-hub edges fixed costs with respect those of M0 ranges, in average, from around 65% for M1 to around 335% for M2 4, increasing for the λ-connected models with the value of λ. This can be easily explained as larger backbone networks are required as λincreases. Nevertheless, as shown by the results of the experiment that we report next, in case of failure, this increase in the design costs strengthens the possibility of being able to re-route all the commodities of a posteriori solutions (after the occurrence of a failure in the inter-hub edges). For this experiment, we have used all n= 20 instances of the TR dataset, whose optimality is guaranteed for all models. For each of the 255 instances generated for the TR dataset with n= 20, we have simulated the following scenarios for potential failures of the inter-hub edges of the backbone networks produced by the different models: •Failure scenario 1 (FS1): Only activated inter-hub edges may fail. Each activated hub edge {k, l}in a backbone network fails (and removed from the backbone network) according to a Bernoulli distribution with probability pkl. V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures 32 Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) de Camargo RS, Miranda Jr G, Ferreira RPM, Luna H (2009) Multiple allocation hub-and-spoke network design under hub congestion. Computers & Operations Research 36(12):3097–3106. Erdo˘gan G, Battarra M, Rodr´ıguez-Ch´ıa AM (2022) The hub location and pricing problem. European Journal of Operational Research 301(3):1035–1047. Ernst AT, Krishnamoorthy M (1996) Efficient algorithms for the uncapacitated single allocation p-hub median problem. Location Science 4(3):139 – 154. Farahani RZ, Hekmatfar M, Arabani AB, Nikbakhsh E (2013) Hub location problems: A review of models, classification, solution techniques, and applications. Computers & Industrial Engineering 64(4):1096– 1109. Gusfield D (1993) Very simple methods for all pairs network flow analysis. SIAM Journal on Applied Mathematics 19(1):143–555. Huang X, Tian YP (2017) Localization in sensor networks with communication delays and package losses. 2017 IEEE 56th Annual Conference on Decision and Control (CDC), 3974–3979 (IEEE). Kim H, O’Kelly ME (2009) Reliable p-hub location problems in telecommunication networks. Geographical Analysis 41(3):283–306. Korani E, Eydi A (2021) Bi-level programming model and KKT penalty function solution approach for reliable hub location problem. Expert Systems with Applications 184:115505. Li Y, Li X, Shu J, Song M, Zhang K (2022) A general model and efficient algorithms for reliable facility location problem under uncertain disruptions. INFORMS Journal on Computing 34(1):407–426. Ma J, Pajouh FM, Balasundaram B, Boginski V (2016) The minimum spanning k-core problem with bounded cvar under probabilistic edge failures. INFORMS Journal on Computing 28(2):295–307. Mar´ın A, C´anovas L, Landete M (2006) New formulations for the uncapacitated multiple allocation hub location problem. European Journal of Operational Research 172(1):274–292. Mohammadi M, Jula P, Tavakkoli-Moghaddam R (2019) Reliable single-allocation hub location problem with disruptions. Transportation Research Part E: Logistics and Transportation Review 123:90–120. O’Kelly ME (1987) A quadratic integer program for the location of interacting hub facilities. European Journal of Operational Research 32(3):393–404. O’Kelly ME (1992) Hub facility location with fixed costs. Papers in Regional Science 71(3):293–306, ISSN 1435-5957, URL http://dx.doi.org/10.1007/bf01434269. Oliveira FA, de S´a EM, de Souza SR (2022) Benders decomposition applied to profit maximizing hub location problem with incomplete hub network. Computers & Operations Research 142:105715. Padberg M, Gr¨otschel M (1985) Polyhedral computations. The Traveling Saliesman Prolbem: A Guided Tour of Combinatorial Optimization, 307–360 (Wiley). V. Blanco, E. Fern´andez and Y. Hinojosa: Hub Location under Link Failures Article submitted to INFORMS Journal on Computing; manuscript no. (Please, provide the manuscript number!) 33 Parvaresh F, Hashemi Golpayegany SA, Moattar Husseini SM, Karimi B (2013) Solving the p-hub median problem under intentional disruptions using simulated annealing. Networks and Spatial Economics 13(4):445–470, URL http://dx.doi.org/10.1007/s11067-013-9189-3. Rodr´ıguez-Pereira J, Fern´andez E, Laporte G, Benavent E, Mart´ınez-Sykora A (2019) The Steiner traveling salesman problem and its extensions. European Journal of Operational Research 278:615–628. Rostami B, K¨ammerling N, Buchheim C, Clausen U (2019) Reliable single allocation hub location problem under hub breakdowns. Computers & Operations Research 96:15–29. Royset JO, Wood RK (2007) Solving the bi-objective maximum-flow network-interdiction problem. INFORMS Journal on Computing 19(2):175–184. Snyder LV, Daskin MS (2005) Reliability models for facility location: the expected failure cost case. Transportation Science 39(3):400–416. Tan PZ, Kara BY (2007) A hub covering model for cargo delivery systems. Networks: An International Journal 49(1):28–39. Wei N, Walteros JL, Worden MR, Ortiz-Pe˜na HJ (2021) A resiliency analysis of information distribution policies over mobile ad hoc networks. Optimization Letters 15(4):1081–1103. Wood RK (1993) Deterministic network interdiction. Mathematical and Computer Modelling 17(2):1–18. Zeng B, An Y, Zhang Y, Kim H (2010) A reliable hub-spoke model in transportation systems. Proceedings of the 4th international symposium on transportation network reliability, Minneapolis, Minnesota, USA, 22–23. Zetina C, Contreras I, Cordeau JF, Nikbakhsh E (2017) Robust uncapacitated hub location. Transportation Research Part B: Methodological 106:393–410. Zheng Y, Ji Y, Shen Y, Liu B, Du Y (2022) Hub location problem considering spoke links with incentivedependent capacities. Computers & Operations Research 148:105959.