Robust post-disaster route restoration

Route restoration is considered to be a task of foremost priority in disaster relief. In this paper, we propose a robust optimization approach for post-disaster route restoration under uncertain restoration times. We present a novel decision rule based on restoration time ordering that yields optimal restoration sequencing and propose conditions for complexity reduction in the model and prove probability bounds on the satisfaction of these conditions. We implement our models in a realistic study of the 2015 Gorkha earthquake in Nepal.


Introduction
In recent years, the increasing intensity and duration of natural disasters have made it primordial for decision-makers, stakeholders and communities to enhance the safety and security of critical infrastructure systems during and after disruptions. Typically disaster-prone critical infrastructure systems, road networks spatially connect one location to another and are enablers of disaster relief by providing accessibility to shelters, hospitals, emergency management centres and so on. These networks are complex and it is a challenge of utmost importance in disaster response for decision-makers to rebuild their functionalities to acceptable limits.

3
Previous events have shown the consequences of natural disasters on road networks. In the aftermath of the 1994 Northridge earthquake four freeway routes were closed due to failures. It is estimated that one fifth of the $6.5 billion total loss in regional economic activity caused by the earthquake could be attributed to transportation disruptions alone (Gordon et al. 1998). In the 1995 Hyogoken-Nanbu earthquake, significant damages to the Hanshin Expressway reduced traffic volumes to 30-55% of pre-earthquake levels, severely impacting disaster relief activities (Chang and Nojima 2001). Following the Haiti earthquake in 2010, relief distribution proved almost impossible, despite the abundance of supplies, due to road infrastructure damages (Çelik 2016). More recently, the 2015 earthquake in Gorkha, Nepal, where the earthquake and a major aftershock triggered landslides around the steep and mountainous regions of Nepal, blocked important road segments, preventing relief supplies from reaching communities (Collins and Jibson 2015).
Due to these potential far-reaching consequences, route restoration is considered to be a task of foremost priority in disaster relief. The route restoration problem (RRP) is the problem of choosing which roads to repair/restore/clear and make traversable following a disaster so that responders can gain access to demand points. In its Public Assistance Debris Management Guide, FEMA pens down debris clearance for emergency route restoration as the first step in disaster response and recovery (Federal Emergency Management Agency 2007). Indeed, adequately planned route restoration can be the precursor of a successful disaster response effort. In the aftermath of the 2011 Tohoku earthquake in Japan, the Tohoku Regional Development Bureau proceeded with urgency to restore roads with a strategy referred to as the "teeth of a comb strategy." The strategy was to first open an inland route along the backbone of the Tohoku Region and then to restore routes to the disaster-stricken areas along the coast (Kazama and Noda 2012). The Tohoku Regional Development Bureau reports the success of prioritizing route restoration in their Earthquake Memorial Museum (see Tohoku Regional Bureau 2014). With the "teeth of a comb strategy", by the day following the earthquake, 11 routes were cleared, allowing the flow of emergency traffic, including ambulances, the police, and medical teams.
The optimization of post-disaster route restoration suffers from three main challenges. The first challenge is the modelling of uncertainties. It is well-established that fitting underlying probabilities to post-disaster damage scenarios is one of the most challenging tasks in disaster management. Disasters interact and are often the result of dynamic historical processes involving multiple stresstriggering events. For example, a main seismic shock can cause aftershocks and induce stresses on nearby fault lines, which leads to increasing earthquake probabilities in nearby region, a phenomenon called earthquake clustering. Successive fissures along the North Anatolian Fault have resulted in a cluster of large earthquakes in the 1900s. The dynamic interactions that lead to complications in estimating earthquake probabilities is well-known (Parsons et al. 2000) and such complications extend beyond seismic events, such as in the Tohoku disaster, where a tsunami caused a nuclear meltdown. The second challenge in the optimization of post-disaster route restoration is model tractability. Optimization models for post-disaster route restorations often involve huge networks of damaged road segments and these can lead to intractable models. In the aftermath of the great Gorkha earthquake in Nepal, planning route restoration for a small rural district north of the Araniko highway required constructing a network of 457 nodes and 555 edges, 66 of which were damaged (Aydin et al. 2018). The third challenge in post-disaster route restoration is in the interpretability of solutions. Disasters are chaotic events and no amount of modelling precision will ever capture the true extent of the damages and perfectly reflect on-site realities. Building overly complex models can lead to optimal decision generation from black boxes, which limits the ability of disaster responders to understand the decision process and react accordingly when unexpected events arise.
In this paper, we tackle these three challenges by proposing a robust optimization approach to optimize the route restoration strategy under uncertain restoration times. Restoration time uncertainty originates from uncertainty in damage or blockage estimation. Studies which consider road damage or blockage uncertainty in disasters abound outside the optimal route restoration context (see Caunhye and Nie 2018;Aydin et al. 2018, for example). Robust optimization is an approach that works by immunizing decision-making against all uncertain data realizations within a deterministic uncertainty set. It is attractive in that it only needs moderate information (specifically, the range of restoration times) about the underlying uncertainty and it is nonparameteric, meaning that it does not require any probability specification. Recent years have seen a tremendous growth in the application of robust optimization for decision-making under uncertainty (Zhen et al. 2018) mainly due to its computational viability and its distribution-free characteristic. In robust optimization, probability distributions are replaced with uncertainty sets (typically conic representable), omitting any information about the distribution except for its support. The method is especially enticing to model extreme events and situations for which parametric probability estimates are hard to obtain. Several seminal papers are highly cited in the robust optimization literature, see (El Ghaoui and Lebret 1997;El Ghaoui et al. 1998;Ben-Tal and Nemirovski 1998;Bertsimas and Sim 2004). Robust optimization models are generally semi-infinite problems that are solved with approximations known as decision rules. A multitude of decision rules, such as the linear (Goh and Sim 2010;Chen and Zhang 2009;Caunhye and Cardin 2018), binary Georghiou 2015, 2018) and finite adaptability (Hanasusanto et al. 2015;Caunhye and Cardin 2017;Bertsimas and Caramanis 2010) decision rules, have been proposed to approximate these models with solvable forms. Even though decision rules are known to be optimal for a few problems, such as some variants of vehicle routing problems (see Gounaris et al. 2013), they potentially sacrifice a significant amount of optimality (Bertsimas and Goyal 2012). In addition, they often require large numbers of variable and constraint additions that can lead to intractability in big post-disaster route restoration models. In this paper, we propose a novel decision rule, based on restoration time ordering, that is optimal for the RRP. We also show that this decision rule can be implemented with little sacrifice to tractability in several variants of the RRP. Furthermore, owing to its simplicity, this decision rule is easily interpretable and offers a practical rule of thumb for route restoration.

3 2 Literature review
In this section, we provide a broad review of the literature on route restoration, starting with general network science methodologies to assess restoration strategies and ending with a specific review of optimization methodologies to plan restoration strategies. Network science methods to tackle the assessment of route restoration strategies typically focus on identifying critical road segments with connectivity measures. The work in Aydin et al. (2019) develops an origin-destination betweenness metric to evaluate post-disaster road performance to assist in decision-making during the recovery process. A different evaluation strategy is proposed in Zhou et al. (2019) with a percolation-based connectivity model used to assess post-earthquake recovery of road networks. In addition to the giant connected component to represent global connectivity, the authors also introduce a metric called local connectivity which is quantified using the number of neighbouring nodes. GIS-based accessibility modelling is also used to assess post-disaster route restoration planning, where the idea is to view route restoration as a means to restore access to services. In Kim et al. (2018a), a scenario-based system dynamics approach is proposed to evaluate the performances of post-disaster debris clearance strategies. The work of Ertugay et al. (2016) and Toma-Danila (2018) also assess the performances of route restoration strategies using accessibility modelling, with the novelty being in the consideration of road closure probabilities.
In optimization modelling for post-disaster route restoration, the focus is on choosing best restoration strategies. Numerous deterministic optimization models have been proposed to such effect. In Moreno et al. (2019), a mixed integer programming problem is proposed to integrate crew scheduling and routing decisions for route restoration. Crew scheduling is also the underlying problem tackled in Kim et al. (2018b), where dynamic damages are addressed and an ant colony heuristic is proposed to solve the resulting mixed integer programming model. In Shin et al. (2019), a similar modelling and solution paradigm is followed, but with the addition of goods delivery to the crew scheduling for road repairs problem. An adaptation of the ant colony heuristic is proposed in Yan and Shih (2012) for the time-space network formulation of the crew scheduling and route restoration problem. The work of Aksu and Ozdamar (2014) approaches route restoration from the perspective of a path-based network optimization problem with the aim to restore access in the network. A decomposition approach is used to improve the tractability of the model. The network-based approach is also used in Akbari and Salman (2017b) and Akbari and Salman (2017a) where mixed integer programming models are used to restore subsets of road segments and reconnect the network. In Akbari and Salman (2017b), the focus is on (1) restoring network connectivity in the shortest possible time and (2) tractability enhancement, via the proposition of a metaheuristic based on relaxation and local search. The focus of Akbari and Salman (2017a) is on an effective decomposition-based solution method for the RRP. In Perrier et al. (2008), the problem is approached as a vehicle routing problem for snow plowing equipment in urban areas to clear roads. A general routing problem is considered for debris removal in Sahin et al. (2016) and a constructive heuristic based on Dijkstra's shortestpath algorithm is proposed to close optimality gaps. The work in Yan and Shih (2009) takes a multi-objective approach to optimize route restoration. Two objectives are considered, namely, minimizing the time spent on restoration and on relief distribution and a heuristic is also proposed as a solution method. In Ajam et al. (2019), latency, which is defined as the time elapsed until a node is visited, is minimized for a debris clearing problem. A meta-heuristic that includes a greedy randomized adaptive search procedure and variable neighbourhood search is proposed to solve the problem. Novel objectives are also proposed in Kasaei and Salman (2016), where two problems are used to plan debris cleaning operations. The first one minimizes the shutdown time of the road system, while the second one focuses on increasing the overall benefits of the reconnecting the network components in a timely manner.
Given the widely recognized intractability of models involving route restoration, as evidenced by the number of reviewed work that propose heuristic solution methods, uncertainty has only rarely been factored into consideration. The work proposed in Çelik et al. (2015) adopts a stochastic approach via a partially observable Markov decision process model to sequence road clearance. The authors recognize the difficulty of incorporating uncertainty considerations and devise a heuristic solution approach as well.

Methodology
This section develops a comprehensive framework, comprising of 6 mathematical programming models, to optimize post-disaster route restoration. The first model is a deterministic model which we use as a baseline for comparisons, especially regarding the impact of uncertainty on route restoration decisions. In the second model, we add restoration time uncertainties to the route restoration problem, via a two-stage robust optimization framework with polyhedral uncertainty set. The first stage, which makes decisions prior to uncertainty realizations, selects the routes to be restored and the second stage sequences route restoration after uncertainty realizations. The two-stage framework generates route selection decisions that are robust to uncertainty and restoration sequencing decisions that are dependent on, and therefore flexible to, uncertainty realizations. The polyhedral uncertainty set has the advantage of (1) being distribution-free and (2) permitting a level of control on the conservatism of our solutions, which is appropriate for disaster settings where probability distributions are hard to estimate and over-and under-conservatism can yield sub-optimal planning. Our third model is a single-stage solvable equivalent of the robust optimization model, which is produced using a novel decision rule based on restoration time ordering.
The first three models implicitly assume that routes can be restored concurrently. Concurrent route restorations require the availability of multiple restoration crews and equipment and in situations where resources are severely limited, this may be impractical. Our fourth model relaxes the concurrent restoration assumption by 1 3 considering sequential job starts. Whilst this would be traditionally achieved using constraint addition, doing so in our case would introduce further complexity in solving the robust counterpart. To make sure that our decision rule remains optimal for the model with sequential job starts and hence, that the tractability and solvability are not impeded, we propose a novel method of relaxing the concurrent restoration assumption using time period enlargement/contraction. Another assumption which may hamper the practical validity of our models is a fundamental one in two-stage robust optimization, which is that restoration time realizations are available prior to route sequencing. In practice, having such restoration time information involves precise post-disaster damage assessment, which is resource-intensive and may be unrealistic in some cases. Our fifth model relaxes this fundamental assumption by optimizing sequencing decisions without a priori restoration time knowledge. Our sixth model is a solvable form of the fifth model under our restoration time ordering decision rule. Altogether, our 6 models follow a methodological underpinning, which is the restoration time ordering decision rule, to propose tractable formulations of the robust route restoration problem for a variety of disaster situations.

Route restoration problem
In the RRP, a disaster has affected a network, resulting in blocked edges and a pressing need for supplies to be delivered to specific locations. Let G = (N, E) be an undirected graph that models the affected network. The set N of nodes consists of a supply node, denoted as 0, a set N d of demand nodes, and a set N ⧵ {{0} ∪ N d } of transshipment nodes. The set E of edges contains a set B of blocked edges, which are edges not traversable post-disaster, and a set U of edges that are still traversable after the disaster.
In our RRP, repair teams are dispatched to restore access from the supply node to the demand nodes. The RRP is the problem of sequencing route restorations over a finite time horizon of T time periods so as to minimize the makespan of access restoration from the supply node to all demand nodes. In this paper, the restoration of access from the supply node to all demand nodes is referred to as "network restoration." Note that the single supply node is not restrictive in where supplies can originate. In the case where there are multiple supply points, they can be represented as nodes connected via dummy edges to a source node 0.
Let [T] be the set of running indices from 1 to T, |S| be the cardinality of set S , and A i be the set of nodes connected to node i. Denoting the number of time periods (integer) needed to restore a blocked edge (i, j) as r ij , we can define our RRP as follow:

Decision variables
x t ij : Is 1 if restoration for edge (i, j) is started at the beginning of period t, 0 otherwise f t ij : Dummy integer flow variable from node i to j in period t z t : Is 1 if graph is fully restored at time t.
where and A i is the set of adjacent nodes to node i.
The objective of Model (P1) is to minimize the makespan of network restoration, reducing degeneracy with the small value . Constraint (1) allows only one restoration to be started in a time period, for simplification. Constraint (2) ensures that a blocked edge may only be unblocked once. Note that we use the terms unblock and restore interchangeably here. Constraint (3) allows flow along an edge if and only if the edge has been unblocked. In this model, flow variables are unrestricted in sign so as to reflect the undirectedness of the graph in a more efficient way. Edges are only defined once, with the sign of the flow variables indicating the direction of the flow. For example, edge (1,2) with a flow of +1 indicates a flow from node 1 to node 2. The same edge with a flow of -1 indicates a flow from node 2 to node 1. If non-negativity restrictions were imposed on flow variables, we would need to define both edge (1,2) and edge (2,1) in the model. Constraint (4) indicates that there is no flow on an edge before its restoration is complete. Constraint (5) is necessary for flow conservation. It also specifies the conditions for the completion of network restoration, which are (1) all demand points have unit net inflows, (2) there is a net (2) outflow of |N d | units from the supply node, and 3) all transshipment nodes have zero net inflows and outflows, meaning that there is a path from the supply node to every demand node. Constraint (6) stipulates that the network must be restored by the end of the planning horizon. Note on Constraint (1) While the constraint allows for a single restoration job to start in any time period, it does not restrict the number of concurrent restoration jobs that can be carried out.
Proposition 1 Given a feasible set of restored edges, B * , sequencing of edge restoration in descending order of the restoration time is optimal.
be a sequence of edge restorations given set B * , where (i n , j n ) is the edge which is chosen for restoration at time n, implying that x n i n j n = 1 . It is clear that the makespan of full restoration, ∑ t∈ [T] tz t , for this sequence can be calculated as {n + r i n j n } . If S is the set of all possible sequences given set B * , the optimal makespan is min It is clear that any other sequence where edge (i n * , j n * ) is chosen for restoration at time n * will have a makespan greater than or equal to n * + r i n * j n * . If edge (i n * , j n * ) is not chosen for restoration at time n * , but is rather interchanged, in the sequence, with edge (i m , j m ) : 1. for m > n * , the makespan increases since m + r i n * j n * > n * + r i n * j n * . 2. for m < n * , the makespan either remains the same or increases since n * + r i m j m ≥ n * + r i n * j n * from the knowledge that r i m j m ≥ r i n * j n * . This shows that for any possible other sequence, the makespan will either remain the same or increase, meaning that a sequence in descending order of restoration time gives the minimum makespan. ◻ An interesting note about Proposition 1 is that it remains valid even if vector r contains real-valued, instead of integer, restoration times. Proposition 1 mainly indicates that there is a simple rule governing the optimal sequencing of edge restorations. It also tells us that if we know the optimal selection of edges to restore, we can sequence edge restorations optimally without running a model.

Illustrative example: baseline case
To showcase the application of Model (P1), an illustrative network is used. Suppose a planner is given a planning horizon of 50 h to restore the network of 13 nodes and 20 edges, of which 14 are blocked, pictured in Fig. 1. The demand nodes are shaded and the number on a traversable edge represents the edge serial number. On blocked edges, the notation a : bh expresses a as the serial number and b as the number of hours required to restore the edge. The value of is chosen to be 1 × 10 −5 . Model (P1) is solved for the illustrative example using ILOG CPLEX. The solution time is 0.09 seconds and the optimal results are shown in Table 1.

Robust counterpart
In practice, restoration times are uncertain. We model the uncertainty in edge restoration times r ∈ ℝ |B| + using the polyhedral uncertainty set where , which varies in the range, [ , 1] , is a parameter, called the "budget of uncertainty", used to control the size of the uncertainty set. When = ∈B r ij and since r ∈ [r,r] , R becomes a singleton wherein R = {r} . As increases, the size of the uncertainty set enlarges. As such, larger total variations in restoration times are allowed, resulting in higher degrees of uncertainty. Notice that we have relaxed the integrality requirement on r to allow more flexibility in defining uncertainty ranges. Our uncertainty set is general in that it allows the planner to define r ij = 0 when doubts exist over whether edge (i, j) has been affected at all.
The robust counterpart of Model (P1) under uncertain restoration times is formulated using a two-stage robust optimization approach as follows:

Decision variables
a ij : Is 1 if edge (i, j) is selected for restoration, 0 otherwise g ij : Dummy integer variable to model flow from node i to j G: Completion time of network restoration x t ij : Is 1 if restoration for (i, j) is started at the beginning of period t, 0 otherwise.
where In the first stage, proactive edge selection is made. A fixed set of edges is selected for restoration after the disaster has struck. This selection is made before edge restoration time realizations. It allows repair teams to start preparations such as gathering necessary equipment, manpower, and other resources. In the second stage, Sequencing of edge restorations is done reactively, that is, subject to restoration time realizations. This means that the planner waits for information on restoration time realizations to decide the sequence in which edge restorations are to be conducted. The robust counterpart selects edges for restoration in such a way that they yield the best worst-case sequence. The term ∑ (i,j)∈B a ij is employed for degeneracy reduction. Constraint (7) ensures that flow is possible on an edge only if it has been chosen for restoration. Constraint (8) ensures if all the edges in the first-stage selection are restored, the network is restored, meaning that the edge selection is a feasible one. Constraint (9) is necessary to calculate the makespan of an edge restoration sequence. Constraint (10) makes it possible to sequence an edge if and only if it has been selected for restoration in the first stage. Constraint (11) allows only one restoration to be started in a time period, consistent with the deterministic model. Notice that the model has complete recourse, meaning that its second stage is always feasible.
In the robust optimization literature, two-stage models are typically solved with decision rules, which are constructs that approximate the solution space for solvability. When second-stage decisions are real-valued, decision rules in affine or polynomial forms are generally employed. When these decisions are integer-valued, finite adaptability or binary decision rules are usually used. Whilst decision rules enable solvability, they are, nonetheless, approximations that generate sub-optimal solutions in most cases. Interestingly, for our case, the two-stage model can be solved exactly, as shown in the following proposition. (10) Otherwise.

Proposition 2 Model (P2) is equivalent to mixed integer programming model
Proof From Proposition 1, we know that the minimum makespan, given a restoration time realization and a feasible set of restored arcs, results when the restoration sequence is in decreasing order of restoration times. By adding constraints to satisfy that decreasing-order condition, we can therefore convert the optimization model Q(a, r) into a feasibility problem. The constraints that model the decreasing-order condition would be Therefore max r∈R Q(a, r) becomes the model where d t ∈ {0, 1} is a binary variable indicating whether network restoration is complete at the beginning of time t or not. The binary variable helps establish the definition of G, the completion time of network restoration. Indeed, . It is clear that the above model can also be expressed as The classical way in robust optimization to proceed from this point is by dualizing the inner maximization problem and inferring the final mixed integer model from strong duality. However, the dualization yields a non-linear model that requires linearization with additional variables. In our model, there is a way to dualize the model that requires fewer additional variables by noting that the decreasing-order condition is expressed as constraints on the objective function of the inner maximization problem. Thus, the model can be re-written as

3
Dualizing the inner maximization problem, and because of strong duality, the above model becomes Combining with the first-stage model, we obtain the model in the proposition. ◻

Illustrative example with uncertain restoration times
Suppose that for the network shown in Fig. 1, the restoration times are uncertain, with the uncertainty expressed as ranges in hours. The ranges are chosen randomly in such a way that they contain the nominal values used for the deterministic network. A value of = 0.6 is used.
The worst-case makespan obtained from the robust optimization model is 14.1 − 1 = 13.1 h, whereas that obtained from the deterministic model is 21.2 − 1 = 20.2 h. Even though the deterministic model yields a smaller number of routes to restore, its worst-case network completion time from optimal sequencing is min G around 7 h longer than that of the robust optimization model. One can easily check that both edge selections make all demand points fully accessible.

Discussions on concurrent restoration jobs
In the models described so far, restoration jobs can be conducted in parallel, meaning that the restoration of the next edge can be started before the completion of the restoration of the previous edge. Table 3 shows the number of concurrent jobs in every time period in our worst-case makespan solutions. As discussed before, the term ∑ (i,j)∈B a ij and in Model (P3) ensures that, among multiple optimal solutions, the one with the smallest number of routes restored is chosen. Starting a restoration job involves equipment and crew setups and these can be avoided by choosing fewer routes. However, as Table 3 shows, fewer routes may also mean more concurrent jobs and the number of jobs is limited by the number of crews and equipment available. By varying , the planner can vary the number of routes restored and as a result, the number of concurrent jobs implied by the solution. For instance, a small negative would maximize the number of routes restored without sacrificing the makespan, if multiple optimal solutions exist. This method, however, is arbitrary in terms of how many concurrent jobs will be carried out. With this method, the planner cannot specify the number of concurrent jobs that his/her resources permit. Such a specification would require the addition of constraints, which will alter the optimality of the restoration time ordering decision rule and thus impact the tractability and solvability of the robust counterpart. In addition, a pre-defined number of concurrent jobs is hard to calculate, given that it depends on the number of equipment available, the manpower at the planner's disposal, the distance of crews and equipment from the disaster area, the contribution of resources from agencies and third parties, and so on.

Sequential restoration jobs
In numerous cases though, simultaneous restorations are impossible to conduct, such as if one road can only be accessed if another has been restored, or if there is a single restoration team at work, or if only one set of equipment is available. Interestingly, in that particular case, it is possible to reformulate the model with the same level of tractability using time period enlargement/contraction. The underlying principle is to redefine time periods such that every job that is started in a time period is completed by the end of it.
Mathematically, this means that r ≤ 1 . If this condition is not met in the original problem description, such as in the network used for our illustrative example, it can be enforced by altering the actual time window that a time period represents. Any RRP can be converted into one with sequential restoration jobs by defining time period t to be a window wherein any job that starts can be completed. In this section, it is helpful to understand r ij as the number of time periods necessary to complete the restoration of edge (i, j). A conservative way to alter the time period definition would be to define it (taken to be in hours for clarity here) as a max (i,j)∈B {r ij }-hour time window. The number of time periods necessary to complete the restoration of edge (i, j) then becomes r ij max (i,j)∈B {r ij } , which is always less than or equal to 1, ∀(i, j) ∈ B.
From a model formulation viewpoint, the robust optimization model (P3) with sequential restoration jobs has an alternative formulation with fewer decision variables and constraints. This simplification is especially useful when networks with large numbers of nodes and edges (especially of blocked edges) are considered. The alternative formulation is derived in the following proposition.
{t + r i t j t } = |B * | + r i |B * | j |B * | . As such, Model (P2) can be expressed as Dualizing the maximization problem max r∈R Q(a, r) , we obtain

and add the constraints
The third constraint, together with the non-negativity of E ij , ensure that when a ij = 1 , E ij = 0 and the first two constraints guarantee that when a ij = 0 , E ij = D ij . With the linearized dual, Model (P2) becomes a single-stage mixed integer programming model. ◻ To apply Model (P4) to the network in Fig. 2, we define our time period as max (i,j)∈B {r ij } = 24.7 hours. The optimal edge selection generated is (9,12), (0,5), (3,6), (1,11), (0,2). The worst-case optimal sequencing is min ∑ (3, 6) → (0, 5) → (0, 2) → (1, 11) → (9, 12) , where → represents "followed by". The worst-case makespan is (5.132 − 1) time periods, which is 4.132 × 24.7 = 102.1 h. This is significantly larger than worst-case makespan obtained when simultaneous job starts are allowed. While it is not a surprise that sequential restoration jobs lead to a larger makespan than concurrent jobs, the size of the difference is noticeably large. The main reason for this big jump is the amount of idle time in between jobs, a direct result of increasing the hours contained per time period from 1 to 24.7. This increase is conservative in that it forces r ij to be less than or equal to 1, for all The conservatism can be reduced by decreasing the time window that a time period represents. Let be the time window that a time period represents. In the above example, = 24.7 h and the optimal objective value (representing the worstcase completion time of network restoration) is , where it is easy to see The worst-case completion time for the same edge selection is therefore We can show that a * is also optimal under the new time window definition (max {(i,j)∈B∶a * ij =1} {r ij } × 24.7) ≤ < 24.7 . We first note that the number of edges selected for restoration under the 24.7-hour time window is the smallest that can be chosen. If a smaller number of selected edges existed, it would have been optimal. Since the worst-case completion time with a * is less than or equal to ∑ (i,j)∈B a * ij + 1 , a larger number of selected edges is also not possible. It is also clear that there is no possible replacement of selected edges that will improve min {(i,j)∈B∶a * ij =1} {r � * ij } and thus, the optimal solution remains the same. For our illustrative example, max {(i,j)∈B∶a * ij =1} {r ij } × 24.7 = 20.2 h is the smallest time window that gives the same optimal edge selection solution. The new optimal objective value is 5 + 0.132 × 24.7 20.2 = 5.161 time periods, which gives a makespan of (5.161 − 1) time periods, which is 4.161 × 20.2 = 84.1 h, a significant improvement from the 102.1 h obtained with the 24.7-h time window. The time window definition is impactful in practice. The re-definition to 20.2 h means that the restoration team moves to the next job after 20.2 h instead of 24.7 h. It also means smaller idle time in between jobs.

Sequencing without a priori restoration time knowledge
In the previous sections, sequencing is done in response to restoration time realizations, with the objective of optimizing the worst-case, over all possible realizations, of sequencing strategy. This a priori knowledge of uncertainty realizations may not necessarily be available for some disaster situations. For instance, resources may not be available to conduct precise damage assessments and produce exact restoration time realizations, or the disaster may be unprecedented. In the chaotic postdisaster environment, unexpected events are bound to happen, such as traffic slowing access to routes, aftershocks from an earthquake derailing response strategies and complicating damage assessments, and interactions of disasters worsening damages, such as the nuclear meltdown that followed the tsunami during the Tohoku disaster. While the two-stage framework is very helpful in guiding decision-making, it may face some practical challenges in more complex situations because of its reliance on restoration time realizations. From a modelling perspective, the lack of a priori restoration time knowledge means that sequencing edge restorations in decreasing order of restoration times is not possible and therefore, that our decision rule is not implementable. The decision-maker has to envisage the possibility of the chosen restoration sequence being such that the edge with the highest restoration time realization is started last in the sequence, delaying the completion time as much as possible. The decision-making problem thus becomes Model (P5) finds the best worst-case sequence, with the understanding that the worst makespan possible by any sequence happens when the edge with the highest restoration time is started last. In traditional robust optimization methods, the inner maximization problem would have been dualized and the resulting single-stage model linearized to obtain the final mixed integer programming problem. However, this augments the complexity of the model with dual variables and additional linearization variables, thus decreasing its scalability. In the following proposition, we prove that there is a more scalable solvable formulation of Model (P5).

Proposition 4 Model (P5) is equivalent to
Proof Let F be a family of feasible sets of restored edges and S B * be the set of all possible edge restoration sequences for B * ∈ F . Define ((i t , j t )) t∈[|B * |] ∈ S B * as an edge restoration sequence, where (i t , j t ) is the edge whose restoration starts in the t th period. Model (P5) can be rewritten in a concise manner as where equality (i) follows from the definition of our uncertainty set and for simplicity, we denote as {(i t , j t )} t>|B * | , as the set of edges not in B * . It is clear that r i t j t = r i t j t , ∀t ∈ [|B|] ⧵ {|B * |} . Therefore, Model (P5) becomes min B * ∈F {�B * � + min (i,j)∈B * {r ij , ∑ (k,l)∈B (r kl − r kl ) + r ij }} (the proof the same logic as that of Proposition 4). This implies that |B w | + W B w ≥ |B a | + W B a and thus, W B w ≥ W B a . From condition 3., W B w = max (i,j)∈B {min{r ij , ∑ (k,l)∈B (r kl − r kl ) + r ij }} , implying that W B w ≥W B a . Since it was shown that W B w ≤W B a , we can conclude that Therefore, B w is also an optimal edge selection set for Model (P2). ◻ Conditions 2., 3., and 4. in Proposition 5 are automatically satisfied in the case where the restoration of all network edges is required. While hitherto our models have focused on disaster response, where the priority is to restore access to demand points, the restoration of all edges is a case of longer-term disaster recovery, where the network is restored back to the state it was in before the occurrence of the disaster. Condition 1. is always satisfied in the case of sequential restoration jobs ( r ≤ 1 ). Therefore, optimizing the restoration of all edges of a network in sequential restoration jobs is similar to optimizing the restoration of all edges of the network without a priori restoration time knowledge.
The purpose of Proposition 5 is to show the conditions under which the planner can solve the reduced Model (P6), which is the solvable form of (P5), instead of the less scalable problem (P3), which is the solvable form of (P2). Without concurrent jobs, Model (P3) reduces to a simpler model anyway, as shown in proposition 3. However, Condition 1. is met in other circumstances as well. It is useful then, to define a probability bound on the makespan position being at |B|, without the restriction r ≤ 1.
. For independent graphs with |B|(≥ 2) blocked edges and independent restoration times where the theoretical range of restoration time upper bounds is given by W , the probability of the optimal worst-case makespan occurring at position |B| has a lower bound of Proof Defining ̄ t be the t th smallest r ij over all (i, j) ∈ B . For independent graphs with |B|(≥ 2) blocked edges, ̄ t is a random variable. Since we know that in Model (P3), the optimal edge restoration sequence is in decreasing order of restoration times and the objective is to optimize the worst-case makespan over all uncertainty realizations, |B| is the makespan position (restoration of all |B| edges required) when |B| +̄ |B| ≥ t +̄ t , ∀t ∈ [|B| − 1] . This can be clearly seen because when meaning that r ij is achievable at the optimal makespan position. We therefore want to find ℙ(|B| +̄ |B| ≥ t +̄ t , ∀t ∈ [|B| − 1]) . Due to graph independence, and assuming the restoration times are independent, this is equivalent to ∏ �B�−1 t=1 ℙ(�B� +̄ �B� ≥ t +̄ t ) . We therefore have where (i) is obtained from Hoeffding's inequality. ◻ The variations of the lower bound on the probability of the optimal worst-case makespan occurring at position |B| with different |B| and W are illustrated in Fig. 3.
We first observe that there is an asymptotic behaviour when |B| increases. The curves become closer to each other because the extra term in ∏ �B� a=2 � 1 − e −(a∕ W) 2 � tends to 1. Figure 3 provides an interesting guide to choosing time windows. When a time window is chosen such that r ≤ 1 , we know that Model (P3) reduces to the simpler Model (P4). Figure 3 shows that for the case where network recovery is required, even when a larger time window is chosen, there is still a minimum probability-which may be high-of a simpler model (notably Model (P6)) providing the 1 − e −(a∕1.5) 2 � = 0.815 probability that its solution will be optimal for Model (P3). This is irrespective of network topology. Such probability bounds are helpful in guiding decision-makers towards the tractable solution of large models. When the restoration of large networks is required, such as following disasters over large areas, decision-makers are often faced with the prospect of solving intractable combinatorial problems to obtain optimal edge sequencing for robustness. By combining the use of both Proposition 5 and Proposition 6, the tractability of the robust optimization model can be significantly improved. The decision-maker first verifies if Condition 2., 3., and 4. in Proposition 5 are satisfied for the network. Thereafter, (s)he can use restoration time bounds to obtain a lower bound on the probability of satisfaction of Condition 1 from Proposition 6. Notice that by varying the time period size, (s)he can guarantee an adequately high lower bound on that probability. The bound represents a confidence level with which the decision-maker can assume that the model without a priori restoration time knowledge has the same solution as the robust optimization model with descending-order decision rule. Since the model without a priori restoration time knowledge has better tractability and scalability than the robust optimization model with decision rule, it means that if the lower bound on the probability of satisfaction of Condition 1. is high, the decision-maker can be confident that the optimal solution of the two models will be equivalent.

Computational experiment: impact of robustness on an Erdős-Rényi random graph
We first investigate the impact robustifying a random network using our robust optimization model with decision rule (P3) and our robust optimization model without a priori restoration time knowledge (P6). We generate the random network using a G(N, p) Erdős-Rényi graph model, where N = 10 is the number of nodes and p = 0.5 is the probability that an edge is included in the graph. The network in Fig. 4 is our generated graph, with the pale-shaded nodes representing demand nodes. All edges are considered to be blocked and their randomly-generated nominal and ranges of restoration times are shown in Table 4. The random generation is carried out as follows: a random number in the range [0, 50] is used as the nominal restoration time. Then random numbers are generated from 0 to the nominal value and from the nominal value to 75, to be used as lower and upper bounds, respectively. We take a time period of 1 hour, a planning horizon of 100 h, = 0.8 , and = 1 × 10 −4 . Node 0 is the supply node. The optimal edge selections from the deterministic model (P1), the robust optimization model (P3), and the robust optimization model without a priori restoration time knowledge (P6) are shown in Fig. 5.
Under nominal restoration times, the edge selections in (P1), (P3), and (P6) can be sequenced for restoration with minimum makespans of 25 h, 41 h, and 32 h, respectively. This shows that our robust optimization models sacrifice nominal performances to hedge against worst-case uncertainty realizations. To compare worstcase performances, we will compare two situations, which are: (1) when there is a priori knowledge of restoration time realizations, meaning that edges can be sequenced according to the descending-order decision rule, and (2) when there is no a priori knowledge of restoration time realizations and therefore, the decision rule is not implementable. For the case where the decision-maker has a priori knowledge of restoration time realizations, the worst-case makespans obtained from the edge selections in (P1), (P3), and (P6) are 62 h, 46 h, and 49 h, respectively. When the decision-maker does not have a priori knowledge of restoration time realizations, the worst-case makespans obtained from the edge selections in (P1), (P3), and (P6) are 70 h, 53 h, and 53 h, respectively This shows the significant improvement in worstcase performances brought about by robust optimization. It also shows that in this particular graph the optimal edge selection for Model (P3) is also optimal for Model (P6) (the converse is not true).
To test the performances of our models outside the nominal and worst-case realizations, we randomly generate 1000 scenarios of restoration time realizations and find the optimal edge sequencing makespans given the edge selections in Fig. 5. The descriptive statistics of the makespans are shown in Table 5. The results show significant improvements in worst-case performances for our robust optimization models, achieved by sacrificing best-case performances. Moreover, our robust optimization models produce lower average makespans and lower makespan variances compared to the deterministic model, showing that by robustifying the network, we have achieved better average performances together with better performance stability, as compared to nominal planning. level of landslides. These surveys, together with landslide magnitude and intensity maps were used to assign minimum and maximum restoration times to closed road segments and the full results are reported in Aydin et al. (2018).
This study focuses on the rural district of Sindhupalchok located on the north side of Kathmandu City. It comprises 457 nodes and 555 edges, 66 of which were rendered impassable due to landslide debris. Access needs to be restored to 81 demand points, which represent settlements in Sindhupalchok. All resource/recovery teams to restore road segments must come from the largest highway in the district, called Araniko Highway. The affected region is show in Fig. 6.
In Aydin et al. (2018), restoration time ranges are taken to be ±2 and this yields non-overlapping ranges for five different classes of landslides in terms of the size of debris, namely, very small, small, medium, high, very high. In our study, we  relax this assumption by taking restoration time ranges to be ±3 , thus allowing greater uncertainty, as well as restoration time overlaps for the different classes of landslides. We also relax the integrality restrictions imposed on the restoration time ranges. The restoration time ranges are shown in Table 6. Taking a time period as 1 hour, a planning horizon of 100 h, = 0.8 , and = 1 × 10 −4 , the robust counterpart (P3) restores access to all demand points by restoring 18 out of the 66 blocked roads. The optimal worst-case restoration completion makespan is 31.2 h. Note that the solution times for all models in this section are less than 10 seconds. Figure 7 maps the worst-case restoration sequence of the 18 restorations.
With sequential restoration jobs with a time period representing 76.8 h, Model (P4) also prescribes the restoration of 18 out of the 66 blocked roads, but the optimal  set of restored edges is different from that of Model (P3). The optimal worst-case makespan is 1386 h and a significant worst-case theoretical idle time of 1292 h. Note: Throughout this work, equal time windows are used, meaning that every time period represents the same time interval. This is to conform to the conventions of optimization modelling. Having variable lengths of time negatively impact optimization models. Because of the equal time windows, the makespan for sequential jobs is much higher than expected. In practice, the restoration crew can move to the next job immediately after finishing the first one and therefore, avoid idle time altogether. In this situation, the makespan becomes 69.3 h (the sum of restoration times of restored edges).

Performance on random scenarios
In this sub-section, we investigate the performances of the deterministic model (P1), the robust optimization (P3) and the model without a priori information (P6) under randomly generated restoration time scenarios. We generate 1000 scenarios of restoration time realizations that are within our uncertainty set and test each model by computing the optimal sequencing and makespan (using the descending order decision rule) for every scenario, given the optimal edge selection of the model as input. The deterministic model is run using the mean restoration times to obtain its optimal edge selection. The results of the experiment are shown in Table 7.
The robust optimization models (P3) and (P6) sacrifice average performances for greater performance stability. We can see that the standard deviations are smaller for these models as compared to the deterministic model. This is because the underlying Fig. 7 Worst-case restoration sequence with our robust optimization approach principle of robust optimization is the search for the best worst-case performance and indirectly, for a reduction in the variation of performances. The robust counterpart (P3) plans for the best worst-case performance assuming that restoration time realizations are available. The model without a priori information (P6) does not use this assumption and therefore suffers from further drops in average performances in the tests. However, because restoration time realizations are not available when model (P6) plans for the best worst-case performance, it generates more conservative plans and therefore offers greater performance stability than all the other models.

Conclusions
In this paper, we propose a robust optimization approach to optimize post-disaster route restoration under uncertain restoration times. We present a novel decision rule based on restoration time ordering that yields optimal restoration sequencing. Under this decision rule, the two-stage robust optimization model becomes a single-stage mixed integer programming problem. We analyze the robust counterpart under two main conditions: 1) restoration can only be performed sequentially and 2) restoration must be performed without a priori restoration time knowledge. We show that under some conditions, the robust counterpart can be reduced to the more tractable and scalable model without a priori restoration time knowledge. These conditions are interestingly, conditions under which 1) restoration of all network edges is required (a network recovery case) and 2) network recovery is complete when the restoration of the edge which is sequenced last is complete. We also prove probability bounds on the satisfaction of the second condition. We implement our models in a realistic study of the 2015 Gorkha earthquake in Nepal and show that with less than a third of blocked roads restored, full access to demand points can be achieved.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.