# Protecting and Expanding Railway Networks to Hedge against Disruptions: A Model with Travel Segments

Authors: Stefano Starita; Chandra Ade Irawan; Pavitra Jindahra
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

Railway systems are instrumental to countries’ economies. They are also vulnerable to man-made and natural disruptions and therefore need to be protected. This paper introduces an optimization model to identify the set of fortifications and network expansions to mitigate the impact of worst-case disruptions. The model incorporates multiple travel segments allowing to study the impact of different route choice behaviours. A Benders decomposition approach is developed to solve the bilevel problem alongside a more efficient nested GRASP heuristic. The algorithms are tested and compared on a set of artificial grids. A case study is developed using the Central London Tube network to provide managerial insights.


---

## Body

## 1. Introduction

Transporting people and goods on rail is crucial for the well-being of countries. Currently, more than 4 billion passengers and 11 billion tonnes per kilometre are moved on rail worldwide [1]. The impact on economic development is further highlighted by the massive investments in railway projects envisioned by the belt and road initiative [2]. Aside from long-distance transportation, the crucial role of rail in moving people around urban areas is especially evident in cities with capillary metropolitan systems, such as New York, London, or Tokyo, with a daily ridership of millions. Developing countries, particularly in South East Asia, are significantly expanding their urban railway networks to address their chronic traffic congestion issues. Bangkok, for instance, is planning on increasing its number of mass rapid transport lines from 5 to 15 and boosting the track length from 159 km to more than 400 km [3]. Similar projects, albeit on a smaller scale, are planned in Hanoi [4]. Building and expanding urban railway systems have also proved to have a beneficial impact on the reduction of air pollution [5]. The size of railway systems makes them vulnerable to disruptions. This is particularly true when the vast portion of the infrastructure is exposed and easy to access. The list of past disruptions is sadly long and includes man-made attacks [6–8], natural events BBC [9], and unintentional disruptions ABC [10]. It is economically unfeasible to protect a high standard of an entire railway system, yet there are mitigating strategies that can be employed. These strategies are threat-dependent. Flood control measures can be implemented for systems vulnerable to flooding, while surveillance and access control measures are suitable to deal with human threats. We group these strategies together and refer to them as fortifications. They can be defined as targeted investments to reduce the vulnerability of an asset (stations, tracks, trains, etc.). Another approach to hedge against disruptions is to increase a railway’s network reliability by expanding its design and adding new connections and paths between stations. The concept of network expansion here is considered generic on purpose. In fact, depending on the stations being connected, an expansion could be a dedicated bus, a walking tunnel, or a line extension. In this work, both fortification and network expansions are considered as mitigation strategies. The disruption of critical infrastructure, like urban railways, can have a dramatic long-lasting socio-economic impact. The low probability of high-consequence disruptive events prompts decision-makers to be risk-averse while devising mitigation strategies [11]. Because of this issue, we focus on worst-case disruptions. This is also further motivated by these systems often targeted by terrorism and other intentional disruptions.One important aspect that needs to be accounted for when focusing on urban railways is the presence of different travel behaviours among a group of customers. Often, users have several choices they can make to reach their destinations. An example is choosing the transportation mode. Even within the same transportation mode, different paths might be available. For example, focusing on railway transportation, the quickest trip and the trip with minimum line changes do not necessarily coincide. The purpose of the trip will affect the choice a traveler will make, and the disutility incurred if the trip has to be cancelled. To the best of our knowledge, no attempt has been made to incorporate customer segments into mathematical models to find strategies to mitigate disruptions. This is crucial as an inaccurate estimation of pre- and postdisruption travel demand can lead to suboptimal solutions, ultimately leaving the railway infrastructure less reliable.To sum up, the contributions of this paper are as follows:(i)
A novel mathematical model to identify railway network fortification and expansion decisions accounting for different customer segments(ii)
Exact and approximate procedures to efficiently solve the aforementioned model(iii)
Managerial insights obtained on the Central London Tube network case studyThe remainder of the paper is organised as follows: In Section2, a review of the relevant literature is provided. Section 3 introduces and describes the mathematical models. Exact and approximate solution approaches are given in Section 4. Computational analysis and managerial insights are discussed in Section 5. Finally, conclusions are drawn along with suggested future research directions in Section 7.


## 2. Literature Review

Over the past two decades, the interest in Operations Research in disruption management has seen a fast growth. Several research branches have appeared, focusing on different critical infrastructures and types of disruptions. Many papers, for instance, assume disruptions to be probabilistic and study the problem of protecting a system with stochastic models. Several examples can be found in supply chain systems [12, 13] and transportation systems [14–16]. In this work, we focus on worst-case disruptions. Two types of game theoretical frameworks are popular within this context: interdiction and fortification [17]. The first one models the problem from an intelligent attacker’s point of view, whose objective is to disrupt a number of system’s assets to inflict the highest damage. This attacker, typically referred to as an interdictor, is used as a proxy to identify the most vulnerable assets of a system. Interdiction models often are modelled as bilevel optimization programs, with one level representing the interdictor another, called the user, computing the system’s performance. The fortification framework adds another level, mimicking a benevolent actor aiming to distribute protective measures to minimize the impact of the interdictor’s actions.Interdiction and fortification problems are described with the leader-follower game theoretical framework. As a result, they are modelled as bilevel or trilevel programming problems. More details are provided in Section3 and Appendix A.These problems have been studied within different contexts (e.g., logistics and supply systems, and electrical grids) and under different assumptions (e.g., complete vs partial protection).In the following, we review the interdiction and fortification research streams that are most relevant to this work.Research on worst-case disruptions in supply chains gained momentum after Church et al. The authors of [18] studied the interdiction problem on p-median and maximal covering problems. Resulting models were customized afterwards, for instance by introducing uncertainty in the intensity of disruptions [19], considering partial disruptions [20], and focusing on hub-and-spoke systems [21]. Most of the efforts are being directed toward fortification problems. Scaparra and Church [22] study the problem of allocating fortification resources on a supply chain facing worst-case disruptions. This work has been extended in several ways. For instance, Liberatore et al. [23] consider the number of disruptions to be uncertain.Shifting the focus to network problems, several seminal works appeared in the past, studying interdiction problems on the flow [24–26] and shortest path networks [27]. These works were more recently extended and specialized. Focusing on flow networks, Wood [28] provides a mathematical programming formulation for the interdiction problem where the attacker aims to maximize the flow on a network. This work is extended by Cormican et al. [29], by assuming uncertainty in the interdiction outcomes and arc capacities. Bell [30] proposes a bilevel formulation to estimate the reliability of a single origin-destination network. In their model, the interdictor does not directly disrupt a link, but it selects the route utilities by choosing from a set of scenarios. Lim and Smith [31] focus on the flow interdiction problem with multiple commodities studying both total and partial disruptions. A dynamic version of the network flow interdiction problem is studied by Rad and Kakhki [32] by adding traversal times to arcs. Schäfer et al. [33] introduce a biobjective formulation to study interdiction on a network with two flows where arcs have two separate capacities. Myung and Kim [34], Murray et al. [35], and Matisziw and Murray [36] study similar flow interdiction problems assuming that origin-destination paths are precomputed and measuring the amount of postdisruption unserved flow. Scaparra et al. [37], and Starita and Scaparra [38] extend these models by adding fortification decisions and assuming that the protection budget is spread on a planning horizon, respectively. Hien et al. [39] study a flow network fortification problem assuming demand uncertainty and considering both arc capacity expansion and arc fortifications to thwart attacks from the interdictor.Focusing on shortest path networks, Israeli and Wood [40] provide a mixed-integer formulation for the interdiction problem where an attacker aims to disrupt arcs in order to increase the distance between two nodes. Bayrak and Bailey [41] propose an interdiction model where the information about arc lengths between the interdictor and the user traveling the network is asymmetric. Sefair and Smith [42] provide a dynamic formulation allowing the interdictor to disrupt arcs once the user travels to any point in the network. Wei et al. [43] introduce a model to minimize the resources of an interdictor aiming to limit the network capacities to a given threshold. The fortification problem on shortest path networks is first introduced by Cappanera and Scaparra [44]. This work is extended by Sadeghi et al. [45] to account for partial disruptions and protections.More closely relevant to this paper is the line of research that studied interdiction and fortification in demand-based networks, while customizing more generic models by incorporating features typical of transportation systems. Alderson et al. [46] and Starita and Scaparra [47] assess vulnerabilities in road networks by studying interdiction problems that account for congestion on the roads, using the popular Bureau of Public Records (BPR) formula. Küçükaydın and Aras [48] study an interdiction model for railways, where both stations and tracks are capacitated and subjected to disruption. Sarhadi et al. [49, 50] propose trilevel fortification models to allocate protective resources on a rail-truck intermodal terminal network. The objective is to identify the set of terminals to fortify so that the postdisruption shipping cost of the container flow is minimized. Focusing on urban railway systems, Jin et al. [51] study a fortification problem where stations are subject to several disruption intensities affecting their capacities. A comparison between vulnerability metrics typically used for railway network risk assessments and interdiction models is provided by Starita et al. [52]. They highlight that interdiction models, while requiring higher computational costs, yield better results. Starita et al. [52] introduce a fortification model that measures the postdisruption railway network’s performance by combining connectivity, demand loss, and travel time.Another related research stream is the one investigating how to design critical infrastructure robust to disruptions. This problem has been widely studied for supply chains, with several papers proposing modelling solutions to design supply networks where facilities can fail with some probability. Within this stream, seminal papers by Daskin [53] and Drezner [54] were followed by many works investigating the problem from different angles. For a more comprehensive discussion on the topic, the reader is referred to Scaparra and Church [55]. While several works can be found on the design of transportation networks robust to uncertain demand [56, 57], only a small number of them explicitly address disruptions of the network’s assets. Among them, Laporte et al. [58] propose an integer program to design a robust railway that provides different routes to address arc failures and limit congestion. Cadarso and Marín [59] introduce a network design model to trade-off between robustness and recovery of a transportation system subject to a set of disruption scenarios. Luo and Xu [60] propose a stochastic model with the aim of minimizing the postdisruption unsatisfied demand of a rail-bus system, by integrating the original network with new additional bus connections. Laporte et al. [61] formulate an interdiction problem to design a railway transit network. They assume that only one link at a time can fail, and the objective is to find the network that maximizes the trip coverage.Within the supply systems context, several authors have proposed models to combine fortification with design measures to mitigate disruptions, both focusing on probabilistic [62, 63] and worst-case disruptions [64–67]. On the other hand, in the distance and flow-based networks, this issue is somewhat neglected. This is particularly true when focusing on worst-case disruptions. Lou and Zhang [68] study both road design and fortification strategies to address random and strategic disruptions, but they do so by introducing two separate optimization models. Fotuhi and Huynh [69] introduce an intermodal freight network expansion design. Their objective is to locate new terminals and devise which of the existing ones should be enhanced or fortified. They assume uncertain demand and probabilistic disruptions. Perea and Puerto [70] propose a model to jointly design and protect a railway network. They use trip coverage as a performance metric while assuming single-arc worst-case disruptions.One further limitation arising from the review of the railway disruption management literature is the typical assumptions on travelers’ behaviour. For example, no paper considers different travel segments. Moreover, route choice models adopted are generally simplistic, assuming users to pick the shortest/cheapest route or entire demand to be served as long as a connection is available. One diversion from this assumption is found in Starita and Scaparra’s study [71]. They propose a fortification model where the amount of travel demand on a given path changes with the length of the shortest available path. However, their model is limited in that the proportion of demand is precomputed as a parameter of a given route and does not change with the available route menu.More complex route choice models can be found in the railway network design literature (e.g., Wang and Meng) [72]. However, to the best of our knowledge, in this stream, no paper addresses disruption issues.Summing up, this work aims to address a few gaps in the disruption management literature for railway systems. It does so by accounting for different traveler segments, providing a methodology that allows for complex route choice behaviour to be incorporated, and using both fortification and design decisions to hedge against disruptions.


## 3. Mathematical Models

This section introduces the interdiction and fortification nonlinear models and their linearised counterparts.The assumptions common to all formulations are listed as follows:(i)
The railway system is modelled as an undirected graph where nodes and edges represent stations and track connections, respectively(ii)
Every origin-destination pair is connected by one or many precomputed paths(iii)
Passengers are grouped into travel segments, defined by different route choice behaviours(iv)
The demand for each segment is known(v)
The problem is studied under worst-case disruption scenarios(vi)
The performance of the railway system is measured in terms of total travel cost/timeIn the following formulations, the railway system is modelled as an undirected graph, with known origin-destination demand, where nodes and edges represent stations and track connections, respectively.


### 3.1. Railway Interdiction with Customer Segmentation (RICS)

The interdiction problem defined below models an intelligent actor whose aim is to disrupt edges so as to maximize the postdisruption economic loss.The model is built with the following assumptions:(i)
A disrupted edge becomes inoperable, rendering unavailable any paths using that edge(ii)
The interdictor can disrupt only a limited predefined number of edgesWe consider a setL of customer segments l, with a demand of dωl for any pair ω. Customers in the same group l share the same likelihood of selecting longer paths and incur the same cost if they cannot make the trip. Each segment is defined by a choice model computing the share of demand traveling on any given path. Specifically, the function Prly∈0,1 represents the share of the customer from the segment l choosing path r, given the route options available (represented by the binary vector y). A segment l is also defined by its penalty unit cost ϕl, which is paid on the share of demand that is lost, identified by P¯ωl. The rest of the notation is summarised as follows:Set and Indices:(i)
E set of edges, indexed by e(ii)
Ω set of origin-destination pairs, indexed by ω(iii)
Rω set of paths connecting pair ω, indexed by r(iv)
Er subset of edges forming the path r(v)
L set of customer segments, indexed by lParameters:(i)
dωl is the total demand for customer segment l and pair ω(ii)
cr is the cost of the path r(iii)
ϕl is the unit cost incurred when a customer belonging to the segment l abandons the trip(iv)
u is the maximum number of edges that can be disrupted simultaneouslyDecision Variables:(i)
xe is equal to 1 if the edge e is disrupted; 0 otherwise(ii)
yr is equal to 1 if the path r is disrupted; 0 otherwiseAn integer programming formulation for the RICS model is given as follows:(1)RICSmaxx,y∑l∈L∑ω∈Ω∑r∈RωdωlcrPrly+ϕldωlP¯ωly,(2)s.t.∑e∈Exe≤u,(3)yr≤∑e∈Erxe∀ω∈Ω,r∈Rω,(4)xe∈0,1∀e∈E,(5)pωi∈0,1∀ω∈Ω,i∈Mω.The objective of the interdictor is to maximize the total expected cost computed as the sum of travel and penalty costs. Constraint (2) states that up to u edges can be removed from the network. Inequalities (3) enforce that path r can be disrupted only if at least one of its edges is removed. Constraints (4) and (5) define variables xe and yr to be binary.


### 3.2. Linearisation via Interdiction Plans Enumeration

With a few exceptions, choice models are typically nonlinear. This assumption renders RICS’ objective function nonlinear. In order to linearise the problem, we propose an approach that relies on enumerating all the possible interdiction plans. Specifically, for each pairω, we enumerate all the path status combinations where each path can be either nondisrupted or disrupted. Such combinations are indexed by i and are referred to as interdiction plans. They are defined via parameter air which takes 1 when path r, connecting pair ω, is disrupted under combination i, and 0 otherwise. Table 1 shows an example of such plans assuming two paths are connecting a pair.Table 1
Example of disruption plans when two paths are available.
ir1r21Nondisrupteda1r1=0Nondisrupteda1r2=02Disrupteda1r1=1Nondisrupteda1r2=03Nondisrupteda1r1=0Disrupteda1r2=14Disrupteda1r1=1Disrupteda1r2=1Each interdiction plani identifies the menu of routes that will be available for the user to choose from. For example, when i=1, both paths r1 and r2 can be used. However, if we focus on i=3, the only path r2 is available. For any of these plans, we can compute the choice probabilities. For example, the share of segment l choosing the path r2, under the disruption plan i=3, is defined by Pr2lyr1=a1r1=0,yr2=a1r2=1.Formally, letMω hold all the interdiction plans for any given ω (i.e., Mω=0,1Rω). For any pair ω and plan i, the precomputed choice probabilities Pirl are defined as Pirl=Prly=ai.Finally, letpωi be the binary variables set to 1 when the plan i is selected for the pair ω. Having precomputed the choice probabilities, variables yr are no longer necessary. A linearisation of RICS is given as follows:(6)RICS−Lmaxx,p∑ω∈Ω∑i∈MωPωi∑l∈L∑r∈RωdωlcrPirl+dωlϕlP¯il,s.t.2,4,(7)∑i∈Mωairpωi≤∑e∈Erxe∀ω∈Ω,r∈Rω,(8)∑i∈Mωpωi=1∀ω∈Ω,(9)pωi∈0,1∀ω∈Ω,i∈Mω.The objective is now to select the interdiction plans maximising the system’s cost. Constraints (7) enforce that a plan i where the path r is disrupted can be selected only if one of the edges along the path r is removed. Constraints (8) state that only one plan must be selected for any pair. Selection variables pωi are defined as binary by (9).


### 3.3. Hedging against Disruptions

In this section, we model a second benevolent actor, referred to as a network planner. Its aim is to select the right strategies to mitigate worst-case disruptions. Mathematically, this results in a bilevel model where the two levels represent the actions of the network planner and the interdictor. These actors have opposite aims, leading to a bilevel min-max model. The network planner can combine two mitigation strategies: fortifying edges and adding new edges.The assumptions of the model are given as follows:(i)
A fortified edge is immune to disruptions(ii)
The network planner can fortify a limited predefined number of edges(iii)
The network planner can add to the network a limited predefined number of new edgesThe model uses the following additional notation:Parameters.(i)
s is the maximum number of edges that can be fortified(ii)
t is the maximum number of new edges that can be built(iii)
E′ is the set of candidate new edges, indexed by e(iv)
Rω′ is the set of additional paths resulting by adding every edge in E′, connecting pair ω, indexed by r(v)
Eω′ is the set containing only the new edges being part of the path r(vi)
Mω′=0,1Rω∪Rω′Decision Variables.(i)
ze is equal to 1 if the edge e is fortified; 0 otherwise(ii)
we is equal to 1 if the edge e is added; 0 otherwiseThe bilevel model to hedge against interdictions (HRICS) is described as follows:(10)HRICSminz,wfx,y,z,w,(11)s.t.∑e∈E∪E′ze≤s,(12)∑e∈E′we≤t,(13)ze∈0,1∀e∈E∪E′,(14)we∈0,1∀e∈E′,(15)fx,y,z,w=maxx,y∑l∈L∑ω∈Ω∑r∈Rω∪Rω′dωlcrPrly+ϕldωlP¯ωly,(16)s.t.∑e∈E∪E′xe≤u,(17)yr≤∑e∈Erxe∀ω∈Ω,r∈Rω∪Rω′,(18)yr≥1−∑e∈Er′weEr′∀ω∈Ω,r∈Rω′,(19)xe≤1−ze∀e∈E∪E′,(20)xe∈0,1∀e∈E∪E′,(21)yr∈0,1∀ω∈Ω,r∈Rω∪Rω′.The objective (10) is to minimize the total expected system’s cost. Up to s edges can be fortified, and t new edges can be added, as respectively stated by constraints (11) and (12). Variables ze and we must be binary (13) and (14). The inner level is the interdiction problem as modelled by RICS, with an additional set of constraints to model the impact of the network planner actions on the interdictor. Constraints (16) and (17) share the same interpretation with (2) and (3), respectively, while accounting for new sets of edges and paths resulting from the network planner’s decisions. Constraints (18) state that a new path is unavailable unless every edge composing it is added. Constraints (19) enforce that fortified edges cannot be disrupted. Finally, as in RICS formulation, xe and yr are binary variables (20) and (21). A linear version of the model can be obtained by adopting the same strategy devised for RICS-L. We refer to this model as HRICS-L, while omitting the formulation for brevity.


## 3.1. Railway Interdiction with Customer Segmentation (RICS)

The interdiction problem defined below models an intelligent actor whose aim is to disrupt edges so as to maximize the postdisruption economic loss.The model is built with the following assumptions:(i)
A disrupted edge becomes inoperable, rendering unavailable any paths using that edge(ii)
The interdictor can disrupt only a limited predefined number of edgesWe consider a setL of customer segments l, with a demand of dωl for any pair ω. Customers in the same group l share the same likelihood of selecting longer paths and incur the same cost if they cannot make the trip. Each segment is defined by a choice model computing the share of demand traveling on any given path. Specifically, the function Prly∈0,1 represents the share of the customer from the segment l choosing path r, given the route options available (represented by the binary vector y). A segment l is also defined by its penalty unit cost ϕl, which is paid on the share of demand that is lost, identified by P¯ωl. The rest of the notation is summarised as follows:Set and Indices:(i)
E set of edges, indexed by e(ii)
Ω set of origin-destination pairs, indexed by ω(iii)
Rω set of paths connecting pair ω, indexed by r(iv)
Er subset of edges forming the path r(v)
L set of customer segments, indexed by lParameters:(i)
dωl is the total demand for customer segment l and pair ω(ii)
cr is the cost of the path r(iii)
ϕl is the unit cost incurred when a customer belonging to the segment l abandons the trip(iv)
u is the maximum number of edges that can be disrupted simultaneouslyDecision Variables:(i)
xe is equal to 1 if the edge e is disrupted; 0 otherwise(ii)
yr is equal to 1 if the path r is disrupted; 0 otherwiseAn integer programming formulation for the RICS model is given as follows:(1)RICSmaxx,y∑l∈L∑ω∈Ω∑r∈RωdωlcrPrly+ϕldωlP¯ωly,(2)s.t.∑e∈Exe≤u,(3)yr≤∑e∈Erxe∀ω∈Ω,r∈Rω,(4)xe∈0,1∀e∈E,(5)pωi∈0,1∀ω∈Ω,i∈Mω.The objective of the interdictor is to maximize the total expected cost computed as the sum of travel and penalty costs. Constraint (2) states that up to u edges can be removed from the network. Inequalities (3) enforce that path r can be disrupted only if at least one of its edges is removed. Constraints (4) and (5) define variables xe and yr to be binary.


## 3.2. Linearisation via Interdiction Plans Enumeration

With a few exceptions, choice models are typically nonlinear. This assumption renders RICS’ objective function nonlinear. In order to linearise the problem, we propose an approach that relies on enumerating all the possible interdiction plans. Specifically, for each pairω, we enumerate all the path status combinations where each path can be either nondisrupted or disrupted. Such combinations are indexed by i and are referred to as interdiction plans. They are defined via parameter air which takes 1 when path r, connecting pair ω, is disrupted under combination i, and 0 otherwise. Table 1 shows an example of such plans assuming two paths are connecting a pair.Table 1
Example of disruption plans when two paths are available.
ir1r21Nondisrupteda1r1=0Nondisrupteda1r2=02Disrupteda1r1=1Nondisrupteda1r2=03Nondisrupteda1r1=0Disrupteda1r2=14Disrupteda1r1=1Disrupteda1r2=1Each interdiction plani identifies the menu of routes that will be available for the user to choose from. For example, when i=1, both paths r1 and r2 can be used. However, if we focus on i=3, the only path r2 is available. For any of these plans, we can compute the choice probabilities. For example, the share of segment l choosing the path r2, under the disruption plan i=3, is defined by Pr2lyr1=a1r1=0,yr2=a1r2=1.Formally, letMω hold all the interdiction plans for any given ω (i.e., Mω=0,1Rω). For any pair ω and plan i, the precomputed choice probabilities Pirl are defined as Pirl=Prly=ai.Finally, letpωi be the binary variables set to 1 when the plan i is selected for the pair ω. Having precomputed the choice probabilities, variables yr are no longer necessary. A linearisation of RICS is given as follows:(6)RICS−Lmaxx,p∑ω∈Ω∑i∈MωPωi∑l∈L∑r∈RωdωlcrPirl+dωlϕlP¯il,s.t.2,4,(7)∑i∈Mωairpωi≤∑e∈Erxe∀ω∈Ω,r∈Rω,(8)∑i∈Mωpωi=1∀ω∈Ω,(9)pωi∈0,1∀ω∈Ω,i∈Mω.The objective is now to select the interdiction plans maximising the system’s cost. Constraints (7) enforce that a plan i where the path r is disrupted can be selected only if one of the edges along the path r is removed. Constraints (8) state that only one plan must be selected for any pair. Selection variables pωi are defined as binary by (9).


## 3.3. Hedging against Disruptions

In this section, we model a second benevolent actor, referred to as a network planner. Its aim is to select the right strategies to mitigate worst-case disruptions. Mathematically, this results in a bilevel model where the two levels represent the actions of the network planner and the interdictor. These actors have opposite aims, leading to a bilevel min-max model. The network planner can combine two mitigation strategies: fortifying edges and adding new edges.The assumptions of the model are given as follows:(i)
A fortified edge is immune to disruptions(ii)
The network planner can fortify a limited predefined number of edges(iii)
The network planner can add to the network a limited predefined number of new edgesThe model uses the following additional notation:Parameters.(i)
s is the maximum number of edges that can be fortified(ii)
t is the maximum number of new edges that can be built(iii)
E′ is the set of candidate new edges, indexed by e(iv)
Rω′ is the set of additional paths resulting by adding every edge in E′, connecting pair ω, indexed by r(v)
Eω′ is the set containing only the new edges being part of the path r(vi)
Mω′=0,1Rω∪Rω′Decision Variables.(i)
ze is equal to 1 if the edge e is fortified; 0 otherwise(ii)
we is equal to 1 if the edge e is added; 0 otherwiseThe bilevel model to hedge against interdictions (HRICS) is described as follows:(10)HRICSminz,wfx,y,z,w,(11)s.t.∑e∈E∪E′ze≤s,(12)∑e∈E′we≤t,(13)ze∈0,1∀e∈E∪E′,(14)we∈0,1∀e∈E′,(15)fx,y,z,w=maxx,y∑l∈L∑ω∈Ω∑r∈Rω∪Rω′dωlcrPrly+ϕldωlP¯ωly,(16)s.t.∑e∈E∪E′xe≤u,(17)yr≤∑e∈Erxe∀ω∈Ω,r∈Rω∪Rω′,(18)yr≥1−∑e∈Er′weEr′∀ω∈Ω,r∈Rω′,(19)xe≤1−ze∀e∈E∪E′,(20)xe∈0,1∀e∈E∪E′,(21)yr∈0,1∀ω∈Ω,r∈Rω∪Rω′.The objective (10) is to minimize the total expected system’s cost. Up to s edges can be fortified, and t new edges can be added, as respectively stated by constraints (11) and (12). Variables ze and we must be binary (13) and (14). The inner level is the interdiction problem as modelled by RICS, with an additional set of constraints to model the impact of the network planner actions on the interdictor. Constraints (16) and (17) share the same interpretation with (2) and (3), respectively, while accounting for new sets of edges and paths resulting from the network planner’s decisions. Constraints (18) state that a new path is unavailable unless every edge composing it is added. Constraints (19) enforce that fortified edges cannot be disrupted. Finally, as in RICS formulation, xe and yr are binary variables (20) and (21). A linear version of the model can be obtained by adopting the same strategy devised for RICS-L. We refer to this model as HRICS-L, while omitting the formulation for brevity.


## 4. Solution Approaches

### 4.1. Benders’ Decomposition for HRICS-L

In this section, we propose an optimal solution approach based on Benders decomposition. While the Benders partitioning technique [73] was initially developed to tackle complex mixed-integer linear problems, it has been extended and generalized to work with nonlinear [74] and bilevel problems [75]. Focusing on bilevel optimization, Benders decomposition has been successfully applied to a wide range of problems. For example, Fontaine and Minner [76] focused on traffic network design and Kazempour and Conejo [77] on strategic generation investment. This decomposition approach is also very frequently used to solve interdiction [40] and fortification [45] problems and to provide a benchmark for heuristic techniques [38, 78, 79].The solution framework relies on first decomposing the two levels into two separate single-level optimization problems and then iteratively solving these problems until a stopping criterion is met. The inner model is built by fixing the fortification and expansion decisions. Solving this problem provides an upper bound to the system’s cost. Each time the inner model is solved, and the resulting interdiction plan is incorporated into a master problem via a set of constraints. This master problem is built to identify the fortification and expansion strategies to thwart all the interdiction plans discovered in previous iterations. At any iteration, the master problem computes the system’s cost accounting for a subset of the entire set of interdiction plans. Therefore, its solution is a lower bound on the cost. A pseudocode description of the algorithm is provided in Appendix B. The formulations for the inner and master problems are discussed as follows.Additional notation is necessary to formulate the two problems. The superscriptk is used to identify values obtained after iteration k, and the parameter it tracks the number of iterations performed by the algorithm. Another set of selection variables qik, similarly to pi variables, is used by the network planner to identify least-cost interdiction plans. For a fixed network expansion plan identified by w^, a new set of paths and edges needs to be considered. This is reflected by the additional notation provided as follows:(i)
Rωw^=Rω∪r∈Rω′:∀e∈Er′,w^e=1(ii)
Ew^=E∪e∈E′:w^e=1(iii)
Mωw^=0,1Rωw^The inner model after iterationk is defined as follows:(22)Innerz^k,w^kmaxx,p∑l∈L∑ω∈Ω∑i∈Mωw^Pωi∑r∈Rωw^dωlcωrPωril+dωlϕlP¯ωil,(23)s.t.∑e∈Ew^xe≤u,(24)∑i∈Mωw^airpωi≤∑e∈Erw^xe∀ω∈Ω,r∈Rωw^,(25)∑i∈Mωw^pωi=1∀ω∈Ω,(26)xe≤1−z^e∀e∈Ew^,(27)xe∈0,1∀e∈Ew^,(28)pωi∈0,1∀ω∈Ω,i∈Mωw^.Innerz^k,w^k is obtained by fixing the network planner decisions in HRICS-L. Solving this model provides an upper bound to HRICS-L, and the solution can be used to generate cutting planes to add to the following Master model. Let us refer to x^k as the optimal solution returned by Innerz^k,w^k and use the binary parameters ρωrk set to 1 if the path r connecting ω is disrupted at iteration k. The master model is aimed at finding the network planner decisions to react to the cutting planes found by the inner subproblem. The formulation is given as follows:(29)Masteritminz,w,qδ,s.t.11−14,(30)δ≥∑l∈L∑ω∈Ω∑i∈Mω′qωik∑r∈RωdωlcωrPril+dωlϕlP¯ωil∀k=1,…,it,(31)∑i∈Mωqωik=1∀k=1,…,it,ω∈Ω,(32)qωik≤∑e∈Erx^ekze∑e∈Erx^ek∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω∪Rω′:air=0∧ρωrk=1,(33)qωik≤∑e∈Er′weEr′∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω′:air=0,(34)qωik≤|Er′|−∑e∈Er′we+∑e∈Erx^ek−∑e∈Erx^ekze∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω∪Rω′:air=1,(35)qωik∈0,1∀k=1,…,it,ω∈Ω,i∈Mω.The objective of the model is to find the combination of protection and network design strategies so as to react to the disruptions selected by the inner subproblem. These disruption strategies are referred to as Benders’ cuts (30). They enforce lower bounds on the objective function δ. For any iteration k and pair ω, only one plan must be selected through variables qωik (31). Constraints (32) enforce that if path r is available under plan i and is disrupted at the kth iteration, then it can be selected only if every edge disrupted (as indicated by x^ek) is fortified. Moreover, if a new path r is available under plan i, then every additional edge must be added (33). Constraints (34) disable any plans where a path r is unavailable as soon as such path is entirely built, fully fortified against x^k or both. Finally, variables qωik are forced to take binary values by (35).


### 4.2. Nested GRASP

While the Benders’ decomposition procedure can find the optimal solution, it also comes with two drawbacks: first, it relies on the linearisation approach discussed in the previous section, and therefore, it requires to precompute a number of scenarios which grows exponentially with the number of paths. Furthermore, as it will be demonstrated in the following section, this procedure poorly scales with the problem’s size. To overcome these issues, a heuristic solution approach is introduced in this section. The choice of adopting a greedy randomized adaptive search procedure (GRASP) by Feo and Resende [80] is motivated by empirical analyses of the optimal solutions which frequently are close to the greedy ones. A nested GRASP (NGRASP) is developed to tackle both network expansion and fortification decisions.A crucial step to enable NGRASP to be efficient is to quickly solve the interdiction level RICS with fixed expansion and fortification strategies. An empirical analysis of the solutions highlights the nonsurprising tendency of the interdictor to fully disrupt the connectivity of high-demand pairs so as to gain on the penalty costs. This motivates the development of a simplified interdiction problem RICS-s to be used as a proxy to RICS. Letσω be a variable taking 1 if the pair ω is disconnected, and 0 otherwise. The formulation for RICS-s is given as follows:(36)RICS−sz^,w^maxx,σ∑ω∈Ωσω∑l∈Lϕldωl,s.t.23,26,27,(37)σω≤∑e∈Ew^xe∀ω∈Ω,r∈Rωw^,(38)σω∈0,1∀ω∈Ω.The model replaces the probabilistic objective of RICS with a simpler function (36) where penalties are paid on the entire demand for every segment once a pair is disconnected. Constraints (37) enforce that a pair is disconnected only if every path serving it (including the newly built) is disrupted. Finally, constraints (38) state that σω are continuous variables taking values between 0 and 1. The rationale is that, although the optimal objectives of RICS and RICS-s are generally quite different, their interdiction strategies are often the same or similar. Solving RICS-s with a MIP solver can then be enhanced by a simple single-edge local search to identify an interdiction strategy which can then be plugged into RICS to compute the actual cost.Letfort, exp, and disr be the sets holding fortified, added, and disrupted edges, respectively, as the algorithm progresses. We assume the notation RICS−sz^,w^ to be equivalent to RICS−sexp,fort, enforcing fort=e∈M:z^e=1 and exp=e∈M′:w^e=1. The pseudocode of NGRASP is described in Algorithm 1.Algorithm 1: NGRASP procedure.
NGRASPu,s,t:exp∗=∅,fort∗=∅,cost∗=∞while Stopping condition not met doexp=∅,fort=∅,cost=∞,improved=truewhile|exp|<tdofor alle∈M′\expdocall GRASPFexp∪e,u,s to obtain coste,forteend forBuild a restricted listRCL adding the bottom μ edges according to costeLete′ be randomly selected from RCLUpdateexp=exp∪e′,fort=forte′,cost=coste′end whilewhileimproveddofor alle∈exp and ∀e′∈M′\extdoSwap inext,e with e′. Call GRASPF exp,u,s to obtain costee′,fortee′Swap back inext,e′ with eend forLete,e′ be the swap yielding the minimum costee′ifcostee′<costthen:improved=true,exp=exp\e∪e′,fort=fortee′,cost=costee′end ifend whileifcost<cost∗then:exp∗=exp,fort∗=fort,cost∗=costend ifend whileAt each iteration, the algorithm performs two main steps. First, it builds a randomised greedy solution and then aims to improve this solution via a local search. The first step sequentially builds the set of edgesexp by randomly choosing from a list of promising candidates. Such a list is obtained by solving an inner GRASP procedure (GRASPF) in which for a given set of expansions, the best found set of fortifications fort and its corresponding cost is returned. The local search described in the code performs 1-edge swaps to further improve the solution. For every odd iteration of the algorithm, this search is replaced by a 2-edges local search. Once the stopping condition is met, best expansion exp∗, fortification sets fort∗, and system’s cost cost∗ are returned. The pseudocode of GRASPF is shown in Algorithm 2.Algorithm 2: GRASPF procedure.
GRASPFext,u,s:fort∗=∅,cost∗=∞while Stopping condition not met dofort=∅,cost=∞,improved=trueSolveRICS−sexp,fort to obtain disrwhile|fort|<sdofor all∀e∈disrdoSolveRICS−sexp,fort, perform 1-edge local search to obtain disre,costeend forBuild a restricted listRCL adding the bottom ν edges according to costeLete′ be randomly selected from RCLUpdatefort=fort∪e′,cost=coste′end whilewhileimproveddofor alle∈fortdoSwap oute and solve RICS−sexp,fort to obtain disrfor alle′∈disrdoSolveRICS−sexp∪e′,fort, perform 1-edge local search to obtain disree′,costee′Swap back inext,e′ with eend forend forLete,e′ be the swap yielding the minimum costee′ifcostee′<costthen:improved=true,fort=fort\e∪e′,fort=fortee′,cost=costee′end ifend whileifcost<cost∗then:fort∗=fort,cost∗=costend ifend whilereturnfort∗,cost∗The algorithm has the same structure as NGRASP with the most noticeable difference in that it makes use of the set of disruptionsdisr to reduce the number of edges to be assessed. This is done by exploiting the obvious fact that in order to thwart a disruption, a fortification strategy should include at least one of the disrupted edges. The list of promising candidates is built by computing the cost by first solving RICS−sexp,fort to obtain disr; then, having all the decision variables fixed, the system’s cost can be computed and used to rank the given edge. Again, a 1-edge swap local search is performed.


#### 4.2.1. Improving NGRASP Efficiency

The overall NGRASP procedure described above faces computational challenges given its nested structure. Each time a candidate expansion link is evaluated multiple GRASPF calls are required. For each of these calls, several mixed-integer problems are solved. A few improving steps can be implemented to reduce the computational burden.NGRASP through its iterations is bound to evaluate combinations of variables that have been already assessed. Storing information in hash tables so that it can be efficiently retrieved leads to significant improvements. The algorithm uses two hash tables. The first takesexp as an input and returns corresponding cost and fort if that set of expansions has been already assessed by previous iterations. Similarly, a second hash table is used with exp and fort as input, returning the corresponding cost and disruption set disr. In other words, the algorithm will call the GRASPF routine only if exp is not in the first hash table. After that, the results will be used to populate the table and avoid recomputation. Similarly, RICS−sexp,fort and the following local search will be performed only if exp,fort is not already in the second hash table.Another simple strategy to enhance the computational efficiency of NGRASP is to limit the neighborhoods that are being investigated by the local searches. In order to do that, simple metrics are computed to rank all the edges, and only a subset of them are swapped in the tentative solutions. For a given (exp, fort, disr), NGRASP computes a score to rank every e \in M'. The score for an edge e is computed summing the demand (i.e., \sum_l d ^ l_\omega) of all the pairs \omega where $e$ is part of an entirely built and nondisrupted path (as defined by exp, fort and disr). As for the local search to improve the solution obtained by solvingRICS−s, the neighborhood is again restricted by computing a simple score for each e∈M∪ext. The score of e is obtained by summing the entire demand of pairs when e is part of any path connecting them and therefore has the potential to disrupt.


## 4.1. Benders’ Decomposition for HRICS-L

In this section, we propose an optimal solution approach based on Benders decomposition. While the Benders partitioning technique [73] was initially developed to tackle complex mixed-integer linear problems, it has been extended and generalized to work with nonlinear [74] and bilevel problems [75]. Focusing on bilevel optimization, Benders decomposition has been successfully applied to a wide range of problems. For example, Fontaine and Minner [76] focused on traffic network design and Kazempour and Conejo [77] on strategic generation investment. This decomposition approach is also very frequently used to solve interdiction [40] and fortification [45] problems and to provide a benchmark for heuristic techniques [38, 78, 79].The solution framework relies on first decomposing the two levels into two separate single-level optimization problems and then iteratively solving these problems until a stopping criterion is met. The inner model is built by fixing the fortification and expansion decisions. Solving this problem provides an upper bound to the system’s cost. Each time the inner model is solved, and the resulting interdiction plan is incorporated into a master problem via a set of constraints. This master problem is built to identify the fortification and expansion strategies to thwart all the interdiction plans discovered in previous iterations. At any iteration, the master problem computes the system’s cost accounting for a subset of the entire set of interdiction plans. Therefore, its solution is a lower bound on the cost. A pseudocode description of the algorithm is provided in Appendix B. The formulations for the inner and master problems are discussed as follows.Additional notation is necessary to formulate the two problems. The superscriptk is used to identify values obtained after iteration k, and the parameter it tracks the number of iterations performed by the algorithm. Another set of selection variables qik, similarly to pi variables, is used by the network planner to identify least-cost interdiction plans. For a fixed network expansion plan identified by w^, a new set of paths and edges needs to be considered. This is reflected by the additional notation provided as follows:(i)
Rωw^=Rω∪r∈Rω′:∀e∈Er′,w^e=1(ii)
Ew^=E∪e∈E′:w^e=1(iii)
Mωw^=0,1Rωw^The inner model after iterationk is defined as follows:(22)Innerz^k,w^kmaxx,p∑l∈L∑ω∈Ω∑i∈Mωw^Pωi∑r∈Rωw^dωlcωrPωril+dωlϕlP¯ωil,(23)s.t.∑e∈Ew^xe≤u,(24)∑i∈Mωw^airpωi≤∑e∈Erw^xe∀ω∈Ω,r∈Rωw^,(25)∑i∈Mωw^pωi=1∀ω∈Ω,(26)xe≤1−z^e∀e∈Ew^,(27)xe∈0,1∀e∈Ew^,(28)pωi∈0,1∀ω∈Ω,i∈Mωw^.Innerz^k,w^k is obtained by fixing the network planner decisions in HRICS-L. Solving this model provides an upper bound to HRICS-L, and the solution can be used to generate cutting planes to add to the following Master model. Let us refer to x^k as the optimal solution returned by Innerz^k,w^k and use the binary parameters ρωrk set to 1 if the path r connecting ω is disrupted at iteration k. The master model is aimed at finding the network planner decisions to react to the cutting planes found by the inner subproblem. The formulation is given as follows:(29)Masteritminz,w,qδ,s.t.11−14,(30)δ≥∑l∈L∑ω∈Ω∑i∈Mω′qωik∑r∈RωdωlcωrPril+dωlϕlP¯ωil∀k=1,…,it,(31)∑i∈Mωqωik=1∀k=1,…,it,ω∈Ω,(32)qωik≤∑e∈Erx^ekze∑e∈Erx^ek∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω∪Rω′:air=0∧ρωrk=1,(33)qωik≤∑e∈Er′weEr′∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω′:air=0,(34)qωik≤|Er′|−∑e∈Er′we+∑e∈Erx^ek−∑e∈Erx^ekze∀k=1,…,it,ω∈Ω,i∈Mω,r∈Rω∪Rω′:air=1,(35)qωik∈0,1∀k=1,…,it,ω∈Ω,i∈Mω.The objective of the model is to find the combination of protection and network design strategies so as to react to the disruptions selected by the inner subproblem. These disruption strategies are referred to as Benders’ cuts (30). They enforce lower bounds on the objective function δ. For any iteration k and pair ω, only one plan must be selected through variables qωik (31). Constraints (32) enforce that if path r is available under plan i and is disrupted at the kth iteration, then it can be selected only if every edge disrupted (as indicated by x^ek) is fortified. Moreover, if a new path r is available under plan i, then every additional edge must be added (33). Constraints (34) disable any plans where a path r is unavailable as soon as such path is entirely built, fully fortified against x^k or both. Finally, variables qωik are forced to take binary values by (35).


## 4.2. Nested GRASP

While the Benders’ decomposition procedure can find the optimal solution, it also comes with two drawbacks: first, it relies on the linearisation approach discussed in the previous section, and therefore, it requires to precompute a number of scenarios which grows exponentially with the number of paths. Furthermore, as it will be demonstrated in the following section, this procedure poorly scales with the problem’s size. To overcome these issues, a heuristic solution approach is introduced in this section. The choice of adopting a greedy randomized adaptive search procedure (GRASP) by Feo and Resende [80] is motivated by empirical analyses of the optimal solutions which frequently are close to the greedy ones. A nested GRASP (NGRASP) is developed to tackle both network expansion and fortification decisions.A crucial step to enable NGRASP to be efficient is to quickly solve the interdiction level RICS with fixed expansion and fortification strategies. An empirical analysis of the solutions highlights the nonsurprising tendency of the interdictor to fully disrupt the connectivity of high-demand pairs so as to gain on the penalty costs. This motivates the development of a simplified interdiction problem RICS-s to be used as a proxy to RICS. Letσω be a variable taking 1 if the pair ω is disconnected, and 0 otherwise. The formulation for RICS-s is given as follows:(36)RICS−sz^,w^maxx,σ∑ω∈Ωσω∑l∈Lϕldωl,s.t.23,26,27,(37)σω≤∑e∈Ew^xe∀ω∈Ω,r∈Rωw^,(38)σω∈0,1∀ω∈Ω.The model replaces the probabilistic objective of RICS with a simpler function (36) where penalties are paid on the entire demand for every segment once a pair is disconnected. Constraints (37) enforce that a pair is disconnected only if every path serving it (including the newly built) is disrupted. Finally, constraints (38) state that σω are continuous variables taking values between 0 and 1. The rationale is that, although the optimal objectives of RICS and RICS-s are generally quite different, their interdiction strategies are often the same or similar. Solving RICS-s with a MIP solver can then be enhanced by a simple single-edge local search to identify an interdiction strategy which can then be plugged into RICS to compute the actual cost.Letfort, exp, and disr be the sets holding fortified, added, and disrupted edges, respectively, as the algorithm progresses. We assume the notation RICS−sz^,w^ to be equivalent to RICS−sexp,fort, enforcing fort=e∈M:z^e=1 and exp=e∈M′:w^e=1. The pseudocode of NGRASP is described in Algorithm 1.Algorithm 1: NGRASP procedure.
NGRASPu,s,t:exp∗=∅,fort∗=∅,cost∗=∞while Stopping condition not met doexp=∅,fort=∅,cost=∞,improved=truewhile|exp|<tdofor alle∈M′\expdocall GRASPFexp∪e,u,s to obtain coste,forteend forBuild a restricted listRCL adding the bottom μ edges according to costeLete′ be randomly selected from RCLUpdateexp=exp∪e′,fort=forte′,cost=coste′end whilewhileimproveddofor alle∈exp and ∀e′∈M′\extdoSwap inext,e with e′. Call GRASPF exp,u,s to obtain costee′,fortee′Swap back inext,e′ with eend forLete,e′ be the swap yielding the minimum costee′ifcostee′<costthen:improved=true,exp=exp\e∪e′,fort=fortee′,cost=costee′end ifend whileifcost<cost∗then:exp∗=exp,fort∗=fort,cost∗=costend ifend whileAt each iteration, the algorithm performs two main steps. First, it builds a randomised greedy solution and then aims to improve this solution via a local search. The first step sequentially builds the set of edgesexp by randomly choosing from a list of promising candidates. Such a list is obtained by solving an inner GRASP procedure (GRASPF) in which for a given set of expansions, the best found set of fortifications fort and its corresponding cost is returned. The local search described in the code performs 1-edge swaps to further improve the solution. For every odd iteration of the algorithm, this search is replaced by a 2-edges local search. Once the stopping condition is met, best expansion exp∗, fortification sets fort∗, and system’s cost cost∗ are returned. The pseudocode of GRASPF is shown in Algorithm 2.Algorithm 2: GRASPF procedure.
GRASPFext,u,s:fort∗=∅,cost∗=∞while Stopping condition not met dofort=∅,cost=∞,improved=trueSolveRICS−sexp,fort to obtain disrwhile|fort|<sdofor all∀e∈disrdoSolveRICS−sexp,fort, perform 1-edge local search to obtain disre,costeend forBuild a restricted listRCL adding the bottom ν edges according to costeLete′ be randomly selected from RCLUpdatefort=fort∪e′,cost=coste′end whilewhileimproveddofor alle∈fortdoSwap oute and solve RICS−sexp,fort to obtain disrfor alle′∈disrdoSolveRICS−sexp∪e′,fort, perform 1-edge local search to obtain disree′,costee′Swap back inext,e′ with eend forend forLete,e′ be the swap yielding the minimum costee′ifcostee′<costthen:improved=true,fort=fort\e∪e′,fort=fortee′,cost=costee′end ifend whileifcost<cost∗then:fort∗=fort,cost∗=costend ifend whilereturnfort∗,cost∗The algorithm has the same structure as NGRASP with the most noticeable difference in that it makes use of the set of disruptionsdisr to reduce the number of edges to be assessed. This is done by exploiting the obvious fact that in order to thwart a disruption, a fortification strategy should include at least one of the disrupted edges. The list of promising candidates is built by computing the cost by first solving RICS−sexp,fort to obtain disr; then, having all the decision variables fixed, the system’s cost can be computed and used to rank the given edge. Again, a 1-edge swap local search is performed.


### 4.2.1. Improving NGRASP Efficiency

The overall NGRASP procedure described above faces computational challenges given its nested structure. Each time a candidate expansion link is evaluated multiple GRASPF calls are required. For each of these calls, several mixed-integer problems are solved. A few improving steps can be implemented to reduce the computational burden.NGRASP through its iterations is bound to evaluate combinations of variables that have been already assessed. Storing information in hash tables so that it can be efficiently retrieved leads to significant improvements. The algorithm uses two hash tables. The first takesexp as an input and returns corresponding cost and fort if that set of expansions has been already assessed by previous iterations. Similarly, a second hash table is used with exp and fort as input, returning the corresponding cost and disruption set disr. In other words, the algorithm will call the GRASPF routine only if exp is not in the first hash table. After that, the results will be used to populate the table and avoid recomputation. Similarly, RICS−sexp,fort and the following local search will be performed only if exp,fort is not already in the second hash table.Another simple strategy to enhance the computational efficiency of NGRASP is to limit the neighborhoods that are being investigated by the local searches. In order to do that, simple metrics are computed to rank all the edges, and only a subset of them are swapped in the tentative solutions. For a given (exp, fort, disr), NGRASP computes a score to rank every e \in M'. The score for an edge e is computed summing the demand (i.e., \sum_l d ^ l_\omega) of all the pairs \omega where $e$ is part of an entirely built and nondisrupted path (as defined by exp, fort and disr). As for the local search to improve the solution obtained by solvingRICS−s, the neighborhood is again restricted by computing a simple score for each e∈M∪ext. The score of e is obtained by summing the entire demand of pairs when e is part of any path connecting them and therefore has the potential to disrupt.


## 4.2.1. Improving NGRASP Efficiency

The overall NGRASP procedure described above faces computational challenges given its nested structure. Each time a candidate expansion link is evaluated multiple GRASPF calls are required. For each of these calls, several mixed-integer problems are solved. A few improving steps can be implemented to reduce the computational burden.NGRASP through its iterations is bound to evaluate combinations of variables that have been already assessed. Storing information in hash tables so that it can be efficiently retrieved leads to significant improvements. The algorithm uses two hash tables. The first takesexp as an input and returns corresponding cost and fort if that set of expansions has been already assessed by previous iterations. Similarly, a second hash table is used with exp and fort as input, returning the corresponding cost and disruption set disr. In other words, the algorithm will call the GRASPF routine only if exp is not in the first hash table. After that, the results will be used to populate the table and avoid recomputation. Similarly, RICS−sexp,fort and the following local search will be performed only if exp,fort is not already in the second hash table.Another simple strategy to enhance the computational efficiency of NGRASP is to limit the neighborhoods that are being investigated by the local searches. In order to do that, simple metrics are computed to rank all the edges, and only a subset of them are swapped in the tentative solutions. For a given (exp, fort, disr), NGRASP computes a score to rank every e \in M'. The score for an edge e is computed summing the demand (i.e., \sum_l d ^ l_\omega) of all the pairs \omega where $e$ is part of an entirely built and nondisrupted path (as defined by exp, fort and disr). As for the local search to improve the solution obtained by solvingRICS−s, the neighborhood is again restricted by computing a simple score for each e∈M∪ext. The score of e is obtained by summing the entire demand of pairs when e is part of any path connecting them and therefore has the potential to disrupt.


## 5. Computational Experience

### 5.1. Networks

To provide a thorough analysis, we build a dataset of square grid networks of different sizes. Figure1 shows an example of a 25-node grid. The existing edges (solid lines) are selected in a way that would replicate urban mass rail transportation systems. Specifically, inner nodes (representing central urban areas) are more interconnected than outer nodes (representing urban outskirts). Candidate edges (dotted lines) are added to allow trading off between connecting peripheral locations and increasing the connectivity of central areas.Figure 1
25-node grid.Traveling times and passengers for edge are uniformly drawn according to Table2.Table 2
Traveling times.
Edges connecting two inner nodesAny other edgeTraveling time (min)U∼2,5U∼5,10Number of passengersU∼500,10000U∼100,1000These values are set in line with the numbers of passenger using underground systems of the large metropolis. For example, the London Underground shows daily passengers on its links ranging from a few hundreds to several thousands. Four network sizes are considered as follows: 16, 25, 36, and 49 nodes. For each size, two instances are built.


### 5.2. Customer Segments and Choice Model

We choose the popular multinomial logit model (MNL) to build travel choice models. LetVrl define the utility that a customer of the segment l associates to the path r, then the probability that such path is chosen is given as follows:(39)Pωrly=1−yreVrl∑r∈Rω∪Rω′1−yreVrl+eV¯l.Therefore, the probability that a customer from the segmentl abandons the trip is given by(40)Pω0ly=1−∑∈Rω∪Rω′=eV¯l∑r∈Rω∪Rω′yreVrl+eV¯l.The utility value of a given pathr is the sum of a base utility β0 and the relative traveling time increment γr between r and the shortest path connecting the same pair. Formally,(41)Vrl=β0+βlγr,where βl is a weight quantifying segment l traveling time sensitivity. In the computational analysis, we consider two customer segments, referred to as commuters and leisure travelers. The first segment represents users that are strongly committed in making the trip and that by abandoning it would incur a significant social-economic damage (e.g., people commuting to work and to educational buildings). The second segment represents users traveling motivated by nonurgent reasons (e.g., shopping); hence, if they were to abandon the trip, the damage would be small. Commuters are assumed to be well informed about their options; hence, they mostly convey on the least-cost available paths. Conversely, we assume higher shares of casual travelers to choose more expensive journeys as a consequence of their reduced knowledge. Detailed choice model parameters are given in Table 3.Table 3
Choice model parameters.
Segmentβ0βlV¯lcωϕlCommuters5−215100Leisure4−41510To clarify the parameters’ choice, Figure2 shows the plot of the probability that a customer would abandon the trip as a function of the γr of two paths. Results are shown for both commuters (2a) and leisure segments (2b). The plots show that commuters are willing to accept longer delays in order to make the trip, while the probability of abandoning the trip for leisure travelers increases more rapidly.Figure 2
Probability of abandoning the trip when pathsr1 and r2 are available.
(a)(b)


### 5.3. Computational Analyses

In this section, the proposed solution approaches are evaluated and compared. Experiments are implemented using Visual Studio and Cplex 12.10 callable libraries. The PC uses an AMD Ryzen 7–3700X processor and 16 GB of RAM. Each combination of budget parametersu,s, and t is tested, letting u to range between 1 and 3, s and t between 1 and 4. For each of the following analyses, a 1-hour limit is enforced on the computing time of both solution algorithms. Moreover, a maximum number of 10 iterations is set for NGRASP. A summary of the results of the analyses on 16- and 25-node grids is shown in Figure 3, More detailed results are reported in Appendix C. while more detailed tables are reported in Appendix C. Each symbol in the picture represents the percentage gap between the objectives returned by Benders and NGRASP, for a specific combination of u,s,t values.Figure 3
Objective gap between NGRASP and Benders. 16- and 25-node networks.Smaller dots denote gaps for 16-node grids, while larger ones show the gaps for 25-node grids. For simplicity, the same symbols are used for the same sized grids so that it is easier to assess the impact of the size of the grids on the algorithms’ gaps. Results are arranged in a matrix of charts, where the 3 rows correspond to the disruption sizes considered, and each column is a combination ofs and t. The top left quadrant, for example, reports on the gaps when u=1,s=1, and t=1,…,4. If less than 4 dots appear for a u,s,t combination, it indicates that instances are reporting the same gaps and their dots overlap. For example, when u=1,t=1, and s=2, all four dots overlap because all gaps are 0%.The figure shows that most of the time, NGRASP retrieves the same solution as Benders (0% gap) which is the optimal one, unless Benders did not converge in time. The dots above0% line indicate that, at times, NGRASP performs worse than Benders. However, the highest gap is 4.9%, and in the majority of cases, the gap is less than 2%. On the other hand, dots below the 0% line show that when Benders fails to converge within the time limits, NGRASP retrieves better solutions. As expected, this happens more frequently once the problem’s complexity increases.With larger grids, the advantage of NGRASP over Benders becomes even more evident. Figure4 reports the gaps on 36 and 49 nodes instances. Similarly to the previous analysis, detailed results are found in Appendix C. The figure shows that when NGRASP underperforms, the gaps are typically small. The increased number of data points below the 0% line suggests that as s and t increase, Benders fails more frequently to converge to the optimal solution, and NGRASP performs significantly better.Figure 4
Objective gap between NGRASP and Benders. 36- and 49-node networks.To compare the computational efficiency of both algorithms, Figure5 and 6 show the average computing times on 16 and 49 nodes instances, respectively. On small networks, the performances are very much alike. Both algorithms become slower when s and t increase. On the other hand, a stark difference arises in the analysis of large networks. In the majority of cases, Benders reaches the 1-hour limit before convergence. Conversely, NGRASP executes all the planned iterations within the time limit in all but 4 cases.Figure 5
Average computational time on 16-node networks.Figure 6
Average computational time on 49-node networks.Overall, the picture drawn by these analyses suggests that, with an equal computational budget, Benders’ algorithm is a suitable tool only for small problems. As the size and complexity of the problem grow, NGRASP is more efficient in retrieving good quality solutions.


## 5.1. Networks

To provide a thorough analysis, we build a dataset of square grid networks of different sizes. Figure1 shows an example of a 25-node grid. The existing edges (solid lines) are selected in a way that would replicate urban mass rail transportation systems. Specifically, inner nodes (representing central urban areas) are more interconnected than outer nodes (representing urban outskirts). Candidate edges (dotted lines) are added to allow trading off between connecting peripheral locations and increasing the connectivity of central areas.Figure 1
25-node grid.Traveling times and passengers for edge are uniformly drawn according to Table2.Table 2
Traveling times.
Edges connecting two inner nodesAny other edgeTraveling time (min)U∼2,5U∼5,10Number of passengersU∼500,10000U∼100,1000These values are set in line with the numbers of passenger using underground systems of the large metropolis. For example, the London Underground shows daily passengers on its links ranging from a few hundreds to several thousands. Four network sizes are considered as follows: 16, 25, 36, and 49 nodes. For each size, two instances are built.


## 5.2. Customer Segments and Choice Model

We choose the popular multinomial logit model (MNL) to build travel choice models. LetVrl define the utility that a customer of the segment l associates to the path r, then the probability that such path is chosen is given as follows:(39)Pωrly=1−yreVrl∑r∈Rω∪Rω′1−yreVrl+eV¯l.Therefore, the probability that a customer from the segmentl abandons the trip is given by(40)Pω0ly=1−∑∈Rω∪Rω′=eV¯l∑r∈Rω∪Rω′yreVrl+eV¯l.The utility value of a given pathr is the sum of a base utility β0 and the relative traveling time increment γr between r and the shortest path connecting the same pair. Formally,(41)Vrl=β0+βlγr,where βl is a weight quantifying segment l traveling time sensitivity. In the computational analysis, we consider two customer segments, referred to as commuters and leisure travelers. The first segment represents users that are strongly committed in making the trip and that by abandoning it would incur a significant social-economic damage (e.g., people commuting to work and to educational buildings). The second segment represents users traveling motivated by nonurgent reasons (e.g., shopping); hence, if they were to abandon the trip, the damage would be small. Commuters are assumed to be well informed about their options; hence, they mostly convey on the least-cost available paths. Conversely, we assume higher shares of casual travelers to choose more expensive journeys as a consequence of their reduced knowledge. Detailed choice model parameters are given in Table 3.Table 3
Choice model parameters.
Segmentβ0βlV¯lcωϕlCommuters5−215100Leisure4−41510To clarify the parameters’ choice, Figure2 shows the plot of the probability that a customer would abandon the trip as a function of the γr of two paths. Results are shown for both commuters (2a) and leisure segments (2b). The plots show that commuters are willing to accept longer delays in order to make the trip, while the probability of abandoning the trip for leisure travelers increases more rapidly.Figure 2
Probability of abandoning the trip when pathsr1 and r2 are available.
(a)(b)


## 5.3. Computational Analyses

In this section, the proposed solution approaches are evaluated and compared. Experiments are implemented using Visual Studio and Cplex 12.10 callable libraries. The PC uses an AMD Ryzen 7–3700X processor and 16 GB of RAM. Each combination of budget parametersu,s, and t is tested, letting u to range between 1 and 3, s and t between 1 and 4. For each of the following analyses, a 1-hour limit is enforced on the computing time of both solution algorithms. Moreover, a maximum number of 10 iterations is set for NGRASP. A summary of the results of the analyses on 16- and 25-node grids is shown in Figure 3, More detailed results are reported in Appendix C. while more detailed tables are reported in Appendix C. Each symbol in the picture represents the percentage gap between the objectives returned by Benders and NGRASP, for a specific combination of u,s,t values.Figure 3
Objective gap between NGRASP and Benders. 16- and 25-node networks.Smaller dots denote gaps for 16-node grids, while larger ones show the gaps for 25-node grids. For simplicity, the same symbols are used for the same sized grids so that it is easier to assess the impact of the size of the grids on the algorithms’ gaps. Results are arranged in a matrix of charts, where the 3 rows correspond to the disruption sizes considered, and each column is a combination ofs and t. The top left quadrant, for example, reports on the gaps when u=1,s=1, and t=1,…,4. If less than 4 dots appear for a u,s,t combination, it indicates that instances are reporting the same gaps and their dots overlap. For example, when u=1,t=1, and s=2, all four dots overlap because all gaps are 0%.The figure shows that most of the time, NGRASP retrieves the same solution as Benders (0% gap) which is the optimal one, unless Benders did not converge in time. The dots above0% line indicate that, at times, NGRASP performs worse than Benders. However, the highest gap is 4.9%, and in the majority of cases, the gap is less than 2%. On the other hand, dots below the 0% line show that when Benders fails to converge within the time limits, NGRASP retrieves better solutions. As expected, this happens more frequently once the problem’s complexity increases.With larger grids, the advantage of NGRASP over Benders becomes even more evident. Figure4 reports the gaps on 36 and 49 nodes instances. Similarly to the previous analysis, detailed results are found in Appendix C. The figure shows that when NGRASP underperforms, the gaps are typically small. The increased number of data points below the 0% line suggests that as s and t increase, Benders fails more frequently to converge to the optimal solution, and NGRASP performs significantly better.Figure 4
Objective gap between NGRASP and Benders. 36- and 49-node networks.To compare the computational efficiency of both algorithms, Figure5 and 6 show the average computing times on 16 and 49 nodes instances, respectively. On small networks, the performances are very much alike. Both algorithms become slower when s and t increase. On the other hand, a stark difference arises in the analysis of large networks. In the majority of cases, Benders reaches the 1-hour limit before convergence. Conversely, NGRASP executes all the planned iterations within the time limit in all but 4 cases.Figure 5
Average computational time on 16-node networks.Figure 6
Average computational time on 49-node networks.Overall, the picture drawn by these analyses suggests that, with an equal computational budget, Benders’ algorithm is a suitable tool only for small problems. As the size and complexity of the problem grow, NGRASP is more efficient in retrieving good quality solutions.


## 6. Case Study

In this section, HRICS is applied to a case study based on the central London Underground network. TFL APIs and open data are used to estimate traveling times, the origin-destination demand matrix, and the share of commuters and casual travelers. The same choice model introduced in Section5 is used. Given the lack of information regarding new candidates’ links, they have been selected with the purpose of filling the geographical gaps in the network. Their traveling times are set using the average traveling time per unit of distance computed using the existing links. The network has 51 stations and 82 links plus 17 additional links for expansion.The case study network is drawn in Figure7. Thicker (gray) lines are used to draw existing connections, while thinner (orange) lines identify the candidate links for expansion. Details of the different tube lines are given in Figure 8.Figure 7
Central London Tube network.Figure 8
Central London Tube network (Lines).The first analysis is aimed at assessing the impact of increasing the budget for protection and to evaluate the benefit of combining fortifications with network expansions. To this end, the graph in Figure9 shows the relative reduction in the objective function as resources b increase from 0 to 10. The concept of resource is maintained generic as it changes across the 3 cases reported in the analysis. Specifically, with the fortification-only case, b represents the number of fortified links (i.e., b=s). Similarly, with the expansion-only case, b is the number of links to add (i.e., b=t). In the combined case, b limits the sum of fortifications and expansions (i.e., b=s+t). Given that for a given b, there are several combinations of s and t such that b=s+t, the graph only shows the objective of the best combination.Figure 9
Relative cost reduction with an increasing budget. Comparison between fortifications only, expansions only, and both.Results show that when budget is limited, adding new links has a better impact than fortifying existing ones. For example, whenb=2, in all 3 disruptive scenarios, adding two links results in cost reductions 5%, 3%, and 4% greater than those obtained by fortifying two links. As the budget increases, the benefit of adding new links becomes marginal and eventually is outpaced by the benefit of fortifications. It is also evident the advantage of incorporating both protection strategies. This is particularly true for small disruptions (i.e., u=1), when, on multiple budget points, cost reductions about 10% greater are achieved. Smaller yet significant improvements can also be observed when u=2,3.The trade-off between fortifications and expansions is further evident in Table4, where the best distribution between the two protection strategies of the budget resources is reported. For small values of b, the majority of the resources are devoted to adding links to the network. This trend changes for larger values of b, where most of the resources are directed to fortifications. The disruption intensity also appears to play a role. When u=1, the budget distribution is more balanced, with an average gap of 0.7 (3.1 fortifications and 2.4 expansions). This gap increases to 1.6 and 3, when u=2 and u=3, respectively.Table 4
Best budget distribution between fortifications and expansions.
ub12345678910Average1s00123335773.1t12222343232.42s00122447793.6t12223231211.93s01233567784.2t11112111211.2Both analyses discussed suggest that the intensity of disruptions leads to different protection strategies. To further analyse this issue, Figure10 highlights how frequent links are chosen for fortification and expansion. The heatmap figures are obtained by solving HRICS with s,t=0,…,4, and darker and thicker lines (green for fortification and orange for expansions) are used to draw the links that are more frequently part of the protection plans. For small disruptions (i.e., u=1), Figure 10(a) shows that fortifications are focused on the Victoria line which is the longitudinal line connecting stations King’s Cross to Victoria. Moreover, a few more links approaching Bank/Monument station along the circle and district line are targeted for fortification. Not surprisingly, links more frequently added are near fortification hotspots, with the sole exception of the link connecting Marble Arch to Edgware Road. When u=2 and u=3, fortification sets become more diverse. A wider focus is given to the central point of the network (i.e., Oxford Circus) being a crossroads of multiple lines. Interestingly, in both cases, when focusing on expansions, there is clear a pattern of providing alternative north-to-south connections, by adding new paths from Angel to Temple (u=2, Figure 10(b)) and from Edgware Road to Victoria (u=3, Figure 10(c)). While differences are evident, there are also links that are consistently selected across the 3 scenarios, suggesting that some portions of the networks are vulnerability hotspots independently on how intense a disruption might be. One clear example is the new link connecting St. Pauls to Mansion House which is almost always selected. It is also clear that the demand is the strongest factor, with popular stations like Bank/Monument, King’s Cross and Victoria, and their neighbours acting like gravity points for protective measures.Figure 10
Links most frequently selected in expansion and fortification plans. (a)u = 1. (b) u = 2. (c) u = 3.
(a)(b)(c)The following analysis is aimed at showing the benefit of incorporating more complex route choice models and multiple travel segments. To this aim, we compare the HRICS model with commuters and leisure travelers (referred to for clarity asHRICScl) against HRICS with only commuters HRICSc and only leisure HRICSl. Moreover, a fourth comparative model HRICSd is used assuming a single travel behaviour where passengers abandon the trip only when a pair is fully disconnected. A cross-comparison of the protection and expansion strategies is performed by first identifying the solutions of each model and then testing these solutions on other models to measure the cost impact.Results Table5 show that HRICSc, HRICSl, and HRICSd all identify solutions that are suboptimal when both customer segments are considered. Although results are not too bad on average, occasionally, these models lead to significant cost increments.Table 5
Relative cost increment when solutions fromHRICSc, HRICSl, and HRICSd are applied to HRICScl.
Protection strategy obtained fromApplied touHRICScHRICSlHRICSdAvg obj incrMax obj incrAvg obj incrMax obj incrAvg obj incrMax obj incrHRICScl10.12%1.10%0.07%0.53%0.11%0.72%20.51%3.69%0.68%4.53%0.51%3.69%30.29%1.98%0.33%1.98%0.26%1.98%Average0.31%2.25%0.36%2.35%0.29%2.13%Finally, we narrow the focus on a single pair (from OxfordCircus to Cannon Street), chosen as it exemplifies the advantages of joint fortification and expansion strategies.Figure11(a) shows 2 paths available under normal conditions. The shortest one starts with the Bakerloo line and switches to the District line at Embankment. As shown in Figure 11(b), both paths are disrupted under the worst-case scenario when u=1, leading to the loss of the entire demand. Fortifying only the disrupted link is not enough as the same disruption outcome can be achieved by disrupting any of the 4 links shared by the 2 paths. If we assume that there are resources to fortify 2 links and add 2 new links (i.e., s=t=2), the best solution will connect Hyde Park Corner with Victoria and St. Pauls with Mansion House and fortify a small path linking Mansion House to Bank via Cannon Street (Figure 11(c)). The expansion strategy generates a new path sharing the last fortified link with the other paths. Combining these protection strategies creates a more robust connection between Oxford Circus and Cannon Street, forcing the interdictor to focus its resources elsewhere.Figure 11
(a) Bond Street to Mansion House,u = 0; s = 0, t = 0. (b) Bond street to mansion house, u = 1; s =0, t=0. (c) Forti_cation and protections, u = 1; s =2; t = 2. (d) Bond street to mansion house, u = 1; s = 0, t = 0.
(a)(b)(c)(d)


## 7. Conclusions and Future Work

This paper introduces a bilevel optimization model to find optimal fortification and design strategies to mitigate the impact of worst-case disruptions on a transportation network. The formulation accounts for different customer segments and travel behaviours. Two solutions approaches are discussed, one based on Benders decomposition, and the other is a nested GRASP heuristic. A computational analysis is run on a set of artificial grids. The decomposition algorithm identifies optimal solutions but poorly scales with the size of the networks. The heuristic is capable of finding good solutions in larger instances. The model is also applied to a case study based on the Central London underground, using real data, where possible.Overall, the analyses point out a few interesting strategic managerial insights: the amount of budget available to hedge against disruptions plays a role in what mitigation strategies should be pursued. When the budget is small and only a few connections can be fortified, then a better strategy is to devote the resources to provide additional paths instead. As the budget increases, the benefit of combining fortifications and expansions becomes clear. These results highlight that mitigation strategies and network planning decisions, where possible, should not be considered as separate problems. Moreover, travel segments and their choice behaviours have a direct impact on where and how the protective resources are focused. Thus, suggesting that developing accurate demand models is an instrumental step for any disruption management problem. Failing to do so is bound to result in suboptimal defensive plans.Several lines of research can be pursued to further extend this work. Accounting for multimodal transportation can allow to capture more realistically the travel. Applying the model to realistic case studies with more detailed travel segments can provide valuable insights. This is particularly true in a postcovid society where the traveling patterns and behaviour are likely to change due to the wide adoption of work-from-home policies. Expanding protection and design decisions to nodes as well is an interesting direction. For example, the current model can be combined with bike-sharing station decisions.Finally, further efforts in developing efficient solution algorithms are needed. This is particularly important if the focus is moved to larger and more complex networks, where the number of paths increases drastically and multiple travel segments are involved.



---

*Source: 9158784-2022-06-06.xml*