Travel time resilience of roadway networks under disaster

Transportation Research Part B 70 (2014) 47–64

Contents lists available at ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

Travel time resilience of roadway networks under disaster Reza Faturechi, Elise Miller-Hooks ⇑ Department of Civil and Environmental Engineering, Univ. of Maryland, 1173, Glenn Martin Hall, College Park, MD 20742, USA

a r t i c l e

i n f o

Article history: Received 6 December 2013 Received in revised form 26 June 2014 Accepted 16 August 2014

Keywords: Disaster management lifecycle Mitigation Multi-stage stochastic programming Stochastic Mathematical Programs with Equilibrium Constraints (SMPEC) Partial User Equilibrium Piecewise linearization

a b s t r a c t A bi-level, three-stage Stochastic Mathematical Program with Equilibrium Constraints (SMPEC) is proposed for quantifying and optimizing travel time resilience in roadway networks under non-recurring natural or human-induced disaster events. At the upper-level, a sequence of optimal preparedness and response decisions is taken over pre-event mitigation and preparedness and post-event response stages of the disaster management life cycle. Assuming semi-adaptive user behavior exists shortly after the disaster and after the implementation of immediate response actions, the lower-level problem is formulated as a Partial User Equilibrium, where only affected users are likely to rethink their routing decisions. An exact Progressive Hedging Algorithm is presented for solution of a singlelevel equivalent, linear approximation of the SMPEC. A recently proposed technique from the literature that uses Schur’s decomposition with SOS1 variables in creating a linear equivalent to complementarity constraints is employed. Similarly, recent advances in piecewise linearization are exploited in addressing nonseparable link travel time functions. The formulation and solution methodology are demonstrated on an illustrative example.  2014 Elsevier Ltd. All rights reserved.

1. Introduction More than 4 million miles of U.S. public roads serve approximately 90% of passenger transport in the country (BTS, 2013). Natural and human-caused hazards threaten this roadway network and create the possibility for significant economic loss due to damage. Damage caused by Hurricane Irene to the Vermont transportation network amounted to $65 million (Lunderville, 2012). The collapse of the I-35 W Bridge over the Mississippi River interrupted more than 140,000 daily vehicular trips causing more than $0.4 million increase in daily passenger trip costs due to traffic rerouting (Zhu et al., 2010). The repair and reconstruction costs of transportation infrastructure systems after Hurricane Katrina were estimated to have exceeded $32 billion (Sundeen and Reed, 2006). Transportation infrastructure systems are also attractive targets for malicious acts. Recent examples include bombings of passenger rail systems in London (2005), Madrid (2004), and Mumbai (2006). Since the 1990s, more than 25% of terrorist attacks have either targeted surface transportation systems or used them to provide access to other targets (MurrayTuite, 2008). In addition to resulting physical damage, these events have long-term socio-economic and psychological impacts. Furthermore, they affect traveler decisions. Gordon et al. (2007) identified a 6% reduction in passenger trips and a sizable shift from public transit services to private automobiles during a two-year period following the 9/11 attack. To prevent significant loss from disaster events, whether caused by a malicious act, an accident or technology failure, or nature, the transportation system must be resilient, and thus able to cope with disaster impact. This notion of resilience was initially conceptualized in relation to stability of ecological systems under a disruptive event (Holling, 1973). Holling ⇑ Corresponding author. Tel.: +1 301 405 2046. E-mail address: [email protected] (E. Miller-Hooks). http://dx.doi.org/10.1016/j.trb.2014.08.007 0191-2615/ 2014 Elsevier Ltd. All rights reserved.

48

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

included in the paper’s resilience definition the system’s ability to bounce back to a state of equilibrium post-event. This definition of resilience as including the system’s ability both to resist and adapt to disruptions has been applied in many settings since then. These settings include transportation and civil infrastructure systems. Quantitative approaches that adopt this definition of resilience in the context of transportation applications include: (Nair et al., 2010; Cox et al., 2011; Chen and Miller-Hooks, 2012; Faturechi et al., 2014). Faturechi and Miller-Hooks (2014) provide a synthesis of approximately 200 works in the literature on resilience and other related measures, including coping capacity, robustness, flexibility and recovery, in the context of transportation. In addition to works that focus on resilience estimation, there are works that determine optimal pre-event mitigation or preparedness strategies (Losada et al., 2012), post-event response actions (Chen and Miller-Hooks, 2012; Vugrin et al., 2010) or both (Miller-Hooks et al., 2012) with the goal of maximizing resilience. A single paradigm for understanding and optimizing resilience and related measures that builds on the existence or nonexistence of possible actions that can be taken pre- or post-disaster is provided in (Faturechi and Miller-Hooks, 2014). All prior works related to the maximization of resilience consider only applications in which resilience enhancing actions are chosen with the aim of achieving a system optimal solution. Such solutions inherently assume that the users of the system will follow these system optimal directives. For example, traffic might be centrally directed to use predetermined routes seeking a system optimum implementation. This is appropriate in many applications, such as in freight networks where the goods to be moved are not cognizant. Several relevant works involving network design under supply or demand uncertainty explicitly recognize the ability of people to make their own decisions regarding their path choice, often with the goal of selfishly maximizing their own utility functions. This is discussed in detail in (Nagurney and Qiang, 2012). These works generally involve a bilevel program structure, where design decisions, such as capacity expansion of a network link, are taken at the upper level, while the response of travelers to the supply offerings is assessed at the lower level. Supply uncertainty typically arises from day-to-day incidents, like traffic accidents, that may cause degradation in network performance. The impact of demand uncertainty is typically measured through variations in travel speeds and, thus, travel times. Chen et al. (2011) provide an extensive review of this literature. A few works in the literature employ a similar bilevel structure in addressing network design or enhancement problems in the context of disaster mitigation. Specifically, these works consider retrofit (Fan and Liu, 2010) and expansion (Lo and Tung, 2003; Dimitriou and Stathopoulous, 2008) actions with the aim of reducing the impact of potential disaster events on network performance. Link capacities are only known with certainty post-disaster. These works build in the capacity uncertainty within the lower-level problem, where the system users take decisions only after the disaster scenario is realized. They propose inexact solution techniques in which the complicating complementarity constraints are relaxed or other heuristics. Through such mitigative actions, system reliability (the probability of continued functionality post-event) or robustness (a post-event functionality measure) are improved. Both are measures of inherent coping capacity and omit from the evaluation the system’s ability to adapt. Faturechi and Miller-Hooks (2014) provide a framework for understanding these and other measures in the context of resilience. In the earlier works where uncertainty in supply (e.g. link capacities) was considered in the lower level, a User Equilibrium (UE) is determined for each potential event scenario given upper-level decisions, and upper-level decisions are taken deterministically. Achieving a UE assumes fully-adaptive behavior by system users. The users are presumed to have perfect information about the state of the roadway network in choosing their paths from their historical travel experiences. Despite a rich literature on travel behavior, modeling such behavior under disruption has received little attention. This lack of attention to post-disaster travel behavior is mentioned in (Zhu et al., 2010). In this context, this assumption of perfect information by all users might be valid only long after the event’s initial occurrence at which time system users have enough information to adapt their travel behavior to the new situation, and a new UE is established. However, shortly after the event occurrence, such an assumption is likely erroneous. The subject of this lower-level problem is the period arising shortly after the occurrence of a disaster event in which short-term, contingency plans can be implemented. According to a user behavior survey of De Palma and Rochat (1999), users have high flexibility in their route choice shortly after the occurrence of an event. That is, user behavior is characterized as being semi-adaptive given limited information on network conditions, including information on damage and completion of repair (Iida et al., 2000). Thus, the lower-level problem is formulated as a Partial UE (PUE) traffic assignment problem. This concept of a PUE was introduced in Sumalee and Watling (2008) for modeling user behavior in post-disaster circumstances. Recently, He and Liu (2012) proposed a prediction–correction traffic assignment model to capture day-to-day travel behavior post-disaster that may also have some relevance. Their model incorporates gradually improving information about network conditions over time. This paper incorporates user behavior in the measurement and maximization of travel time resilience for roadway networks given under a set of possible disaster scenarios. The problem of quantifying and optimizing travel time resilience (i.e. the Travel Time Resilience Problem (TTRP)) is formulated as a bilevel, three-stage stochastic program with lower-level equilibrium (PUE) constraints. Both upper- and lower-level problems involve capacity uncertainty. The upper-level includes a three-stage decision making process in which both pre- and post-event resilience enhancing actions may be taken. The decision process is informed as information is revealed at each stage, and is compatible with the Disaster Management Life-Cycle (DMLC) (Waugh, 2000): (1) pre-event expansion and retrofit as mitigation options to enhance the coping capacity of the road network, (2) pre-event preparedness where resources are acquired and prepositioned shortly in advance of a predicted event occurrence to facilitate response actions, and (3) post-event short-term response actions taken post-disaster to

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

49

restore network capacity, minimize the extent of damage, and/or protect the remaining facilities. A multi-hazard perspective is taken, whereby actions that may be effective in one scenario may be ineffective in another. An exact Progressive Hedging Algorithm (adapted from Rockafellar and Wets (1991)) is presented for solution of a single-level equivalent to this bilevel, three-stage stochastic program. Whether addressing day-to-day incident-induced traffic congestion or disaster events, including pre-event or both preand post-event actions for enhancing system performance, or employing a UE or PUE, these problems involving uncertainty in available system capacity in which user response to network supply decisions is captured can be mathematically modeled as Stochastic Mathematical Programs with Equilibrium Constraints (SMPECs). Thus, they are a type of Stochastic Network Design Problem (SNDP). A general discussion on the properties of SMPECs and potential solution techniques can be found in (Patriksson, 2008). In addition to its contributions to resilience measurement, this paper extends the study of SNDPs (SMPECs) generally. Key to its contributions are its consideration of supply uncertainty in both upper- and lower-level problems, incorporation of a three-stage stochastic program in the upper-level to capture key relevant DMLC stages, use of a PUE rather than a UE in traffic assignment as is appropriate for the disaster-context, and application of cutting-edge linearization methodologies for dealing with complementarity constraints and nonlinear, nonseparable travel time functions, as well as nonlinear design decision terms. The contributions of this work are derived from: (a) development of a SMPEC-based framework for conceptualizing a measure of resilience compatible with the stages of the DMLC, (b) using a Partial UE to model post-disaster behavior of roadway network users, (c) transformation of the proposed SMPEC into a Stochastic Mixed Integer Program (SMIP) through the application of cutting-edge piecewise linearization techniques, and (d) adaptation of an exact solution technique based on the Progressive Hedging Algorithm to solve the transformed problem. It must be noted that a number of assumptions were required in developing the proposed model and solution schemes. First, only supply-oriented uncertainties were taken into account. That is, travel demand is assumed to be inelastic to the disaster impact. Second, a PUE adequately captures driver route-choice behavior post-disaster. The next section introduces the problem formulation. This is followed by description of the solution method in Section 3 and application on an illustrative example in Section 4. 2. Problem formulation At the upper level of the proposed TTRP, mitigation, preparedness and response actions are chosen with information from the lower level about the resulting total travel time for all O–D pairs that can be expected given upper level choices. The upper level acts as the leader, determining the optimal supply decisions. The lower level acts as the follower, wherein affected system users selfishly determine their paths with knowledge of the upper-level decisions. The optimal solution to the bilevel problem results at a Stackelberg equilibrium (Gibbons, 1992). The upper-level problem is a three-stage stochastic program that accounts for the occurrence of one of a set of potential disaster events, as well as the information that is revealed about these events and their consequences over time. In the first stage, possible disaster events and their consequences are known probabilistically. At the end of this stage, after some passage of time, certain attributes of the disaster event may be revealed, informing the second stage or equivalently creating a second information state. Here, this information includes the disaster event type and its temporal and spatial properties/distributions. The third information stage arises, revealing a final information state, once the event has occurred and quick assessment of the disaster region has been completed. This information process can be captured through the concept of a scenario tree (Table 1) within which each node represents an information state and each link carries with it the probability of transitioning from one information state to the next over time (or stages). See (Heitsch and Römisch, 2009) for additional background on scenario tree development. Decisions are taken at each stage and are, thus, captured at the nodes as depicted in Table 1. The tree captures all possible outcomes, and each path from the root of the tree to a leaf (i.e. from the first stage to the last) gives a possible scenario. Where a finite set of possible disaster scenarios is considered, the scenarios can be completely enumerated and are known a priori. Notation employed in describing the travel time resilience problem is given next, followed by the problem formulation. A network representation of the roadway system is exploited. The network’s topology is given by G = (N, A), where N is the set of nodes and A is the set of links. Associated with each link is its travel time and capacity limitation, both of which are random variables. Network performance is measured under a set of possible disaster scenarios, each of which is defined by a disaster event type, affected links and its impact on the travel time and capacity of these links. A network representation of the roadway system is exploited. The network’s topology is given by G = (N, A), where N is the set of nodes and A is the set of links. Associated with each link is its travel time and capacity limitation, both of which are random variables. Network performance is measured under a set of possible disaster scenarios, each of which is defined by a disaster event type, affected links and its impact on the travel time and capacity of these links. 2.1. A measure of travel time resilience in roadway networks Total travel time, tt, to serve a given O–D demand is chosen as the system-level measure of performance. The disruption profile given in Fig. 1, graphically captures the variation of a roadway network’s tt over the DMLC, from pre-event conditions

50

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

Table 1 DMLC stages and evolution in uncertainty.

⁄

Pr(npjnp1) = probability of an event occurring at stage p given information revealed in stage p  1.

Fig. 1. Travel time disruption profile for passenger traffic in roadway network, tto, ttd, and ttr are total travel times at end of pre-event, confusion, and response stages.

in which a UE is reached until recovery is complete and a new UE state is established given new network conditions. Changes in the network conditions may result from long-term activities, such as reconstruction. Users adapt their travel behaviors to this new situation, creating a new UE. Immediately after the occurrence of a disaster event, i.e. post-event (confusion), capacity is degraded and users may not be able to ascertain the disaster’s impact on potential routes. Thus, they may be confused or indecisive, creating inefficiencies in the use of a network that is simultaneously performing below its norm. Post-response, travel times improve, reaching a PUE, as the network capacity is partially restored to a satisfactory level through the implementation of short-term repairs, and users have received limited information on network conditions and improvements. Note in the figure travel demand is presumed to be inelastic. The resilience of passenger traffic in roadway networks is defined as the network’s ability to resist and adapt to disruption (adapted from Faturechi and Miller-Hooks, 2014). Total travel time is employed in assessing resilience, RT,B, under a given budget, B, for taking mitigation, preparedness and response actions and given time allotted for response action implementation, T. Given the negative relationship between the network’s system-level performance and its total travel time, as in Eq. (1), the reciprocal of total travel time achieved in reaching a PUE at the end of the response stage divided by the reciprocal of the total travel time achieved in a UE pre-event and pre-action is taken to quantify resilience.

51

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

RT;B ¼

ttr

1

o1

tt

¼

hxo ; to i hxr ; tr i

ð1Þ

2.2. Notation The notation employed within the mathematical program is as follows: Sets A N W Kw S P

set of links, a, in the roadway network set of network nodes, n, representing roadway intersections and points of demand set of origin–destination (OD) pairs w set of paths k between OD pair w 2 W set of disaster types s 2 S set of stages, p = 1, 2, 3, over which information is gained: p = 1 refers to a pre-event stage where the current state is known deterministically, but possible future disaster events in terms of type, location and consequences are known only probabilistically; p = 2 refers to a later point in time when the event type and location are known deterministically (either pre- or post-event), but the impact on the system is unknown, and the probability distribution of the event’s potential or perceived consequences is updated; and p = 3 refers to the point in time after the event has occurred and all event characteristics, as well as the system state, are known deterministically Modeling parameters np information process capturing the state of knowledge, n, about the system’s current and future states at stage p 2 P: n1 captures current conditions deterministically. Future conditions are known only probabilistically in stage 1. n2 specifies event type and location. In stage 2 probability distributions of future conditions are updated conditioning on the event type and location. n3 specifies conditions of the network once the event is fully realized D vector of OD travel demand, D = [. . ., Dw, . . .]T "w 2 W h iT o o fo(n1) known vector of pre-event path flows, f ðn1 Þ ¼ . . . ; f k;w ðn1 Þ; . . . 8k 2 K w ; w 2 W h iT xo(n1), known vectors of pre-event link flows and capacities, xo ðn1 Þ ¼ . . . ; xoa ðn1 Þ; . . . and co ðn1 Þ ¼ o 1 c (n ) h iT . . . ; coa ðn1 Þ; . . . 8a 2 A, respectively h iT vector of pre-event link travel times, to ðn1 Þ ¼ . . . ; t oa ðn1 Þ; . . . "a 2 A to(n1) h iT cd(n3) vector of post-event link capacities under information state n3 ; cd ðn3 Þ ¼ . . . ; cda ðn3 Þ; . . . 8a 2 A

D

link-path incidence matrix, D = [Da,k,w] "a 2 A, k 2 Kw, w 2 W (Da,k,w = 1 if path k 2 Kw uses link a, and =0 otherwise)

K

OD pair-path incidence matrix, K = [Kk,w] for "k 2 Kw, w 2 W(Kk,w = 1 if path k connects OD pair w and Kk,w = 0 otherwise) post-event damage state matrix of paths under information state n3, P(n3) = [pk,w(n3)] "k 2 Kw, w 2 W(pk,w(n3) = 1 if the path k 2 Kw is affected under n3, and pk,w(n3) = 0 otherwise) disaster type matrix H(n3) = [. . .,Hs(n3), . . .]T, where Hs(n3) = 1 if when reaching information state n3 the disaster event that has occurred is of type s, Hs(n3) = 0 otherwise

P(n3) H(n3)

available budget B T response time 1st stage variables T c1(n1) 1 vector of first-stage action variables, c1(n1) = [c1,e(n1), c1,h(n1)]T where c1;e ðn1 Þ ¼ ½. . . ; c1;e a ðn Þ; . . . is the T

1 vector of link capacity expansion levels "a 2 A, and c1;h ðn1 Þ ¼ ½. . . ; c1;h a;s ðn Þ; . . . is the vector of link retrofit levels "a 2 A, s 2 S. Since n1 is revealed from the start of the decision horizon, c1(n1) is written as c1 for simplicity. 2nd stage variables c2(n2) vector of disaster-specific link preparedness (resource availability – second stage) action levels given h iT information state n2, c2 ðn2 Þ ¼ . . . ; c2a;s ðn2 Þ; . . . 8a 2 A; s 2 S, where c2a;s ðn2 Þ ¼ 1 if appropriate preparedness

action for disaster type s is taken in link a, and c2a;s ðn2 Þ ¼ 0 otherwise (continued on next page)

52

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

3rd stage variables c3(n3) vector of disaster-specific link response (third-stage) levels under information state n3, h iT c3 ðn3 Þ ¼ . . . ; c3a;s ðn3 Þ; . . . 8a 2 A; s 2 S, where c3a;s ðn3 Þ ¼ 1 if appropriate response action for disaster type s is taken in link a, and c3a;s ðn3 Þ ¼ 0 otherwise

xr(n3), cr(n3) tr(n3)

h iT vectors of post-response link flows and capacities under information state n3, xr ðn3 Þ ¼ . . . ; xra ðn3 Þ; . . . and h iT cr ðn3 Þ ¼ . . . ; cra ðn3 Þ; . . . 8a 2 A, respectively h iT vector of post-response link travel times as a function of link flows and capacities, tr ðn3 Þ ¼ . . . ; tra ðn3 Þ; . . .

8a 2 A

fr(n3)

sr(n3)

h iT r r vector of post-response path flows under information state n3 ; f ðn3 Þ ¼ . . . ; f k;w ðn3 Þ; . . . 8k 2 K w ; w 2 W h iT vector of post-response path travel times under information state n3 ; sr ðn3 Þ ¼ . . . ; srk;w ðn3 Þ; . . .

8k 2 K w ; w 2 W ur(n3)

h iT vector of post-response shortest travel times under information state n3 ; ur ðn3 Þ ¼ . . . ; urw ðn3 Þ; . . .

8w 2 W

2.3. TTRP Next, the TTRP is formulated as a bi-level, three-stage, stochastic, nonlinear program for maximizing travel time resilience of roadway networks. The program involves stochasticity in both upper- and lower-levels. Upper-level problem

( max

c1 2C1

En2 jn1

" max En3 jn2 ;n1

c2 2C2 ðc1 ;n2 Þ

#) max RT;B ðn3 Þ

ð2Þ

c3 2C3 ðc2 ;n3 Þ

( ) i X   d 3  o h 1;h 3 3 d 3 3 ; 8a 2 A s:t: cra ðn3 Þ ¼ 1 þ c1;e ðn Þ þ H ðn Þ c  c ðn Þ c þ c ðn Þ c s a a a a a;s a;s

ð3Þ

s

q3a;s ðn3 Þ 6 T; 8a 2 A; s 2 S ( ) h i X 1;e X 1;h 2 3 ba þ Hs ðn3 Þ ba;s þ ba;s ðn2 Þ þ ba;s ðn3 Þ 6 B

ð4Þ ð5Þ

s2S

a2A

n   1;e 1;e In the upper-level formulation, the first-stage feasibility set for mitigation actions is given by C1 ¼ c1 c1;e a 6 ca 6 ca ; 1;h 1;e 1;e a for link a as lower- and upper-bounds on capacity expansion level, respectively. 0 6 ca;s 6 1; 8a 2 A; s 2 S:g with ca and c The retrofit variable c1;h a;s sets the desired fortification level for link a for each disaster type s. More than one retrofit action can 1;h be taken on the same link. c1;h a;s ranges between 0 and 1, where ca;s ¼ 1 refers to the highest fortification level obtainable for link a such that no damage will be incurred as a consequence of the occurrence of a disaster of type s. The range on second  stage preparedness levels is defined by the feasibility set C2 ðc1 ; n2 Þ ¼ fc2 ðn2 Þ0 6 c2a;s ðn2 Þ 6 1; 8a 2 A; s 2 Sg, where c2a;s ðn2 Þ ¼ 1

means all resources required to repair damage following a disaster of type s are provided in advance upon realizing infor n o  mation state n2. C3 ðc2 ; n3 Þ ¼ c3 ðn3 Þ0 6 c3a;s ðn3 Þ 6 1; 8a 2 A; s 2 S is the third-stage response level feasibility set for information state n3, where c3a;s ðn3 Þ ¼ 1 infers that capacity along link a is restored to its pre-event capacity. The objective function (2) seeks to maximize the expectation of network resilience over all possible scenarios given by each possible information state n3. This objective function form takes advantage of the fact that the numerator in Eq. (1) of the resilience measure RT,B(n3), xo, to, is constant. Therefore, objective function (2) can be replaced by Eq. (6) where a minimization is employed in place of maximization:

( min En2 jn1

c1 2C1

" min

c2 2C2 ðc1 ;n2 Þ

En3 jn2 ;n1

#) min

c3 2C3 ðc2 ;n3 Þ

r

3

r

3

hx ðn Þ; t ðn Þi

ð6Þ

Throughout, the realization of n2 is conditioned on n1 and similarly n3 is conditioned on and n2 and n1. To simplify the notation, with the exception of the expectation operator, these conditional statements are assumed to be implicit. Postresponse capacity of link a, cra ðn3 Þ, is defined in Eq. (3) as a function of the link’s pre-action, pre-event capacity, coa , and its post-event capacity under information state n3 ; cda ðn3 Þ, accounting for link expansion, retrofit and response decisions that

53

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

are taken in the first and second stages. The effects of decisions are presumed to be linear to the original link capacities. Inclusion of parameter Hs (n3) ensures that the effects of specialized link retrofit and response actions have effects that are consistent with disaster type s and associated information state n3. Second-stage link preparedness actions do not directly affect link capacity, and are not included in the equation. q3a;s ðn3 Þ in constraints (4) is defined as the implementation time of response action of type s along link a under information state n3. These constraints guarantee that, for each information state n3, all response actions that are to be taken are completed before the end of the response period, i.e. by time T. The budget limitation is assured through constraint (5). Link expansion, retrofit, preparedness, and response costs used in constraint (5), as well as response implementation times of constraints (4) are assumed to be functions of action level as described through constraints (7)–(11). 1;e 1;e  c1;e ; 8a 2 A ba ¼ b a a 1;h 1;h  c1;h ; 8a 2 A; s 2 S b ¼ ð1 þ c1;e Þ  b

ð7Þ

2 ba;s ðn2 Þ

ð9Þ

a;s

a

a;s

ð8Þ

a;s

2 ðn2 Þ  c2 ðn2 Þ; 8a 2 A; s 2 S ¼b a;s a;s " # n h i o co  cd ðn3 Þ 3 3 3 3 3 2 1;e 3 3 3 2 a a   ba;s ðn Þ ¼ ð1 þ ca Þ ba;s ðn Þ  ba;s ðn Þ  ba;s ðn Þ ca;s ðn Þ c3a;s ðn3 Þ; 8a 2 A; s 2 S coa " # n h i o co  cd ðn3 Þ 3 3 3 3 2 3 1;e 3 3 3 2 a a a;s ðn Þ  q a;s ðn Þ  qa;s ðn Þ ca;s ðn Þ qa;s ðn Þ ¼ ð1 þ ca Þ q c3a;s ðn3 Þ; 8a 2 A; s 2 S coa

ð10Þ ð11Þ

1;e and b 1;h are first-stage unit costs of expanding link a or, for a given disaster event of type s, retrofitting link a, where b a a;s respectively. 2 ðn2 Þ The implications for first-stage retrofit costs of link expansion are captured in constraints (8). In constraints (9), b a;s denotes second-stage unit costs of link preparedness actions for a given disaster event of type s and information state n2. Third-stage unit costs and times required for implementing response actions are defined in constraints (10) and (11), respectively. They also account for link expansion when undertaken in the first stage. Both are functions of response and preparedness levels, wherein the effects of preparedness in advance of an event affect the efficiency of post-event response actions. 3 ðn3 Þ½q 3a;s ðn3 Þ are minimum (maximum) post-disaster costs [times] of complete reconstruction of link a b3a;s ðn3 Þ½q3a;s ðn3 Þ and b a;s having realized information state n3. In the former, it is presumed that no preparedness actions were taken ðc2a;s ðn2 Þ ¼ 0Þ, while in the latter all appropriate preparedness actions were taken ðc2a;s ðn2 Þ ¼ 1Þ. Taking a preparedness action in advance 3 ðn3 Þ and of an event can enable response actions at lower cost and implementation time. Thus, b3 ðn3 Þ 6 b a;s

a;s

3a;s ðn3 Þ. Incurred costs or implementation times are also a function of damage-level. Thus, an additional term, q3a;s ðn3 Þ 6 q coa cda ðn3 Þ , coa

is included in constraints (10) and (11) to capture the increasing expense and effort required to address situations

with higher damage levels. Implementation time limitations may be required in the second stage, where it may be the case that preparedness actions must be implemented within only a couple of days or less. The amount of time depends on the impending hazard event classification and forecasting capabilities. In the formulation as stated, it is assumed that there are sufficient resources for the estimated costs to complete chosen preparedness actions within this short time frame. If desired, one can add additional constraints to the second stage to restrict the preparedness action implementation times as is done in the third stage with respect to response actions. This inclusion of preparedness implementation time restrictions will not impact the general form of the formulation and will not affect the applicability of proposed solution methods. The link flows needed to compute the objective function of the upper-level problem are determined through solution of the lower-level problem (12) for each information state n3 arising post-event. The lower-level problem seeks a flow pattern that achieves a PUE given actions taken in solution of the upper-level problem. In a PUE, user behavior is characterized as semi-adaptive and assumes that only those who are affected are likely to reconsider their original route decisions. Lower-level problem:

min

xr 2Xrx ðn3 Þ

XZ a

0

xra ðn3 Þ

t ra ðv ; cra ðn3 ÞÞdv ; r

ð12Þ o

r

o

where Xrx ðn3 Þ ¼ fxr ðn3 Þjxr ðn3 Þ ¼ Df ðn3 Þ þ D½I  diagPðn3 Þf ; f ðn3 Þ ¼ diagPðn3 Þf , Kfo = D, fr(n3) P 0} is a feasibility vector 2 3 .. 0 0 7 6 . 3 7 set of post-response link flows in which the diagonal matrix diagPðnÞ ¼ 6 4 0 pk;w ðn Þ 0 5, "k 2 Kw, w 2 W, and .. . 0 0 2 3 .. 6 . 0 0 7 7 I¼6 4 0 1 0 5 is the identity matrix of the same size. The formulation is path-based and is adopted from Sumalee and .. . 0 0

54

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

Watling (2008). xr(n3) = Dfr(n3) + D[I  diagP(n3)]fo defines the link flow as the summation of the post-response flows along affected paths and pre-event flows along unaffected paths using that link. Only affected users will reroute post-response, i.e. fr(n3) = diagP(n3)fo. The affected links are determined in the scenario generation process. Affected users are those whose original paths use affected links. The affected paths are input to the model through the post-event damage state matrix of paths under information state n3, P(n3). Total post-response flows on affected paths between an OD pair w equals the total pre-event flow between that OD pair, preserving the original demand. Inelasticity in the demand assumption is also captured through flow conservation constraints, Kfo = D, wherein pre-event flows are set in accordance with a fixed demand vector, =[. . ., Dw, . . .]T "w 2 W. This assumption is consistent with expectations described in (Iida et al., 2000) in which traffic demand patterns are said to adjust to pre-event patterns regardless of repair state in the settlement stage (immediately after the disaster event). 3. Solving the TTRP The bilevel TTRP (2)–(12) is solved by first reducing it to a single-level problem as is often done in the literature. To accomplish this, the lower-level problem is eliminated and corresponding Karush–Kuhn–Tucker (KKT) conditions are embedded within the upper-level problem. Larsson and Patriksson (1995) showed, in a similar context involving a bilevel program with a related UE formulation in the lower level, that use of the KKT conditions in place of the lower-level problem is both necessary and sufficient for optimality. Their proof can be directly extended to this application, and applies equivalently for the PUE principle employed in the formulation herein. This is because the PUE formulation in (12) differs from the comparable UE formulation only in the specification of the feasibility set and the feasibility set in (12) remains convex. The resulting single-level program is a three-stage stochastic program with nonlinear, but convex objective and nonconvex constraints. Obtaining a globally optimal solution to such a program is formidable. Thus, linear approximations are employed. With these approximations, a single-level, three-stage Stochastic Mixed Integer Linear Problem (SMILP) results. Steps taken to obtain the SMILP are presented in Sections 3.1, 3.2, and 3.3. An exact solution technique based on concepts of the Progressive Hedging Algorithm initially introduced in (Rockafellar and Wets, 1991) is presented in Section 3.4. 3.1. Single-level TTRP The single-level TTRP can be reformulated to include the upper-level problem taken over the lower-level feasibility set

Xrx ðn3 Þ, along with additional UE constraints (14) and (15): ( " min En2 jn1

min

c1 2C1

c2 2C2 ðc1 ;n2 Þ

En3 jn2 ;n1

min

c3 2C3 ðc2 ;n3 Þ;xr 2Xrx ðn3 Þ

r

3

#) r

3

hx ðn Þ; t ðn Þi

ð13Þ

s.t. (3)–(5)

 T r f ðn3 Þ sr ðn3 Þ  ur ðn3 Þ ¼ 0;

ð14Þ

sr ðn3 Þ  ur ðn3 Þ P 0;

ð15Þ T

T

where sr ðn3 Þ ¼ D t ðn3 Þ ¼ ½. . . ; srk;w ðn3 Þ; . . . is the vector of post-response path travel time, and ur ðn3 Þ ¼ ½. . . ; urw ðn3 Þ; . . . is the vector of post-response shortest paths, "k 2 Kw, w 2 W, given information state n3. 1 r

3.2. Linear approximations 3.2.1. Complementarity constraints UE constraints (14), which model the restriction that a path take flow only if it is the shortest path, are transformed into mixed integer constraints through Schur’s decomposition (Horn and Johnson, 1985) using Special Ordered Sets of Type 1 (SOS1) variables (Siddiqui and Gabriel, 2013). Alternative methods often use a disjunctive constraint approach (e.g. Wang and Lo (2010)) employing binary variables and a large constant. However, this approach requires extensive computational resources and its results are highly sensitive to the constant’s value. Specific to the TTRP, applying Schur’s decomposition on constraints (14) results in Eqs. (16)–(18). r

f ðn3 Þ þ ½sr ðn3 Þ  ur ðn3 Þ ; 2 r 3 f ðn Þ  ½sr ðn3 Þ  ur ðn3 Þ Hr ðn3 Þ ¼ ; 2

Gr ðn3 Þ ¼

3 T

ð16Þ ð17Þ

3 T

Gr ðn3 ÞGr ðn Þ  Hr ðn3 ÞHr ðn Þ ¼ 0; r

3

r

3

ð18Þ r

3

r

3

r

3

r

3

where G (n ) and H (n ) are Schur’s decomposition vector functions. Since f (n ), s (n )  u (n ) P 0, G (n ) P 0. Thus, only the positive square root of Gr(n3)  Gr(n3)T is feasible and bilinear constraints (18) can be reformulated as in (19).

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

Gr ðn3 Þ  jHr ðn3 Þj ¼ 0:

55

ð19Þ r

To eliminate the absolute value function, jH (n)j, constraints (19) are transformed through the introduction of additional þ  SOS1 variables Hr ðnÞ and Hr ðnÞ, at most one of them can be non-zero. þ



Gr ðnÞ  ½Hr ðnÞ þ Hr ðnÞ ¼ 0:

ð20Þ

3.2.2. Objective function The objective (13) of the TTRP seeks to minimize the expectation of total travel time incurred along the shortest time paths over all O-D pairs. The objective requires the product of flow and travel time variables, and thus, is nonlinear. A technique for linearizing the objective function introduced by Wang and Lo (2010) is employed. This technique exploits common travel time properties of active paths for each O–D pair under UE conditions. Given fixed demand vector, D = [. . ., Dw, . . .]T, and information state n, objective (13) can be replaced by (21).

( min En2 jn1

c1 2C1

" min

c2 2C2 ðc1 ;n2 Þ

En3 jn2 ;n1

#) min

c3 2C3 ðc2 ;n3 Þ;xr 2Xrx ðn3 Þ

hD; ur ðn3 Þi

:

ð21Þ

3.2.3. Link travel time function The link travel time is estimated using the Bureau of Public Roads (BPR) function (Eq. (22)). Given post-response link flow and capacity random variables for a given link a 2 A, this function is a two-dimensional, nonseparable function.

tra ðn3 Þ

¼

t 0a ðn3 Þ

þ ma

" #na xra ðn3 Þ cra ðn3 Þ

; a 2 A;

ð22Þ

where t 0a ðn3 Þ is the link free flow travel time, and ma and na are BPR function parameters (herein, na = 4 for all a 2 A). In this paper, a novel piecewise linearization technique introduced by Vielma and Nemhauser (2011) for general multidimensional functions is applied in linearizing this link travel time function. It has been shown to outperform other existing piecewise linearization techniques (Vielma et al., 2010). The following describes the application of this technique for the TTRP. The procedure begins by partitioning link flow and capacity variable domains. The travel time domain is thus defined by a two-dimensional flow-capacity domain. In general, any point in an n-dimensional domain can be uniquely represented by a convex combination of n + 1 points (Carathéodory, 1911). For two-dimensions, three points, thus, are required, and therefore, the link travel time domain can be partitioned into non-overlapping triangles. Any flow-capacity pair falls within a single triangle and is given by a convex combination of the associated triangle’s corner-point coordinates. The link travel time values associated with the corner-points are directly calculated using Eq. (29). Vielma and Nemhauser’s (2011) method identifies the active triangle containing the flow-capacity pair under consideration through the introduction of binary variables and additional constraints, and approximates the corresponding travel time through a convex combination of the travel time values at its corner-points. € ¼ ½. . . ; € c ¼ ½. . . ; €ca;j ; . . .; 8i; Let xra ðn3 Þ and cra ðn3 Þ be represented by the vector of segments x xa;i ; . . . and € j 2 V a ¼ f0; 1; . . . ; v a g, respectively, where va is a power of two. Note that an extra constraint may be required to prevent searching partitions added to create a large enough square to cover the feasible region for the specific problem given a large enough va. The domain of the corresponding travel time function will be [0, va]2 with the segments represented within the   matrix €t ¼ . . . ; t ra ð€ xa;i ; €ca;j Þ; . . . , "a 2 A. This domain is triangulated using the J1 Union Jack triangulation approach (originally proposed by Todd (1977)). Fig. 2 graphically depicts the J1 Union Jack triangulation of the two-dimensional domain of the link travel time, where the domain [0, va]2 is covered by copies of a 2  2 size square (highlighted in Fig. 1), each encompassing eight triangles. The entire domain is partitioned into 2v 2a triangles, accordingly. As shown in Fig. 2, there are groups of white and gray triangles such that each square contains one triangle of each color. xra ðn3 Þ and cra ðn3 Þ are formulated as convex piecewise-linear functions of € xa;i and €ca;j points, respectively, in Eqs. (23)–(26).

xra ðn3 Þ ¼

X

€xa;i xa;i;j ðn3 Þ; 8a 2 A;

ð23Þ

€ca;j xa;i;j ðn3 Þ; 8a 2 A;

ð24Þ

i;j2V a

cra ðn3 Þ ¼ X

X

i;j2V a

xa;i;j ðn3 Þ 6 1; 8a 2 A;

ð25Þ

xa;i;j ðn3 Þ P 0; 8a 2 A; i; j 2 V a ;

ð26Þ

i;j2V a

where xa,i,j(n3) is a convex combination weight under information state n3 for a 2 A, i, j 2 Va. Accordingly,

tar ðn3 Þ ¼

X

i;j2V a

t ra ð€xa;j ; €ca;j Þxa;i;j ðn3 Þ; xa;i;j ðn3 Þ P 0:

ð27Þ

56

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

Fig. 2. J1 Union Jack triangulation of link travel time domain.

The active triangle is determined through a binary branching scheme of logarithmic depth in three steps. In the first step, an independent SOS1 type branching is employed to select the active column of 1  1 size squares containing the active triangle. Let V a ¼ f1; . . . ; v a g be the set of columns in the link travel time domain for a 2 A. A corresponding set is defined as La = {1, . . ., log2va} containing a logarithmic number of columns. Let Ba : V a ! f0; 1glog2 v a be a bijective function with special structure such that Ba(j) and Ba(j + 1) are allowed to be different in at most one vector element "j 2 Van{va}. Let r(Ba) be the support of vector Ba. SOS1 type branching is implemented on the logarithmic set La to find the active column.

X

X

xa;i;j ðn3 Þ 6 X a;l ðn3 Þ; 8a 2 A; l 2 La ;

j2V a i2J þ ðl;Ba ;v a Þ 2

X

X

xa;i;j ðn3 Þ 6 ½1  X a;l ðn3 Þ; 8a 2 A; l 2 La ;

ð28Þ

j2V a i2J 0 ðl;Ba ;v a Þ 2

X a;l ðn3 Þ 2 f0; 1g; 8a 2 A; l 2 La ; 0 where Vþ a ðl; Ba ; v a Þ ¼ fj 2 V a : 8i 2 V a ðjÞ; l 2 r½Ba ðiÞg and V a ðl; Ba ; v a Þ ¼ fi 2 V a : 8i 2 V a ðjÞ; l R r½Ba ðiÞg are valid sets for SOS1 constraints. See (Vielma and Nemhauser, 2011) for more background of SOS1 constraints. Next, a similar SOS1 type branching is employed to select the active row of 1  1 size squares which contains the active triangle.

X

X

xa;i;j ðn3 Þ 6 C a;l ðn3 Þ; 8a 2 A; l 2 La ;

i2V a j2J þ ðl;Ba ;v a Þ 2

X

X

xa;i;j ðn3 Þ 6 ½1  C a;l ðn3 Þ; 8a 2 A; l 2 La ;

ð29Þ

i2V a j2J þ ðl;Ba ;v a Þ 2

C a;l ðn3 Þ 2 f0; 1g; 8a 2 A; l 2 La : Given a square (i.e. equal number of columns and rows) in the link travel time domain V a ¼ f1; . . . ; v a g, a 2 A, can also be used to represent the rows. The active square is, thus, determined through the selection of active columns and rows. In the final step, the active triangle is determined. A single binary variable, Ia(n3), for a 2 A, is introduced in the following constraints to identify the color (gray or white) of the active triangle (Ia(n3) = 1 if the active triangle is white, and Ia(n3) = 0, otherwise).

XX

xa;i;j ðn3 Þ 6 Ia ðn3 Þ; 8a 2 A;

i2Ea j2Oa

XX

xa;i;j ðn3 Þ 6 ½1  Ia ðn3 Þ; 8a 2 A;

ð30Þ

i2Oa j2Ea

Ia ðn3 Þ 2 f0; 1g; 8a 2 A; where Ea = {0, 2, . . ., va}  Va and Oa = {1, 3, . . ., va  1}  Va are subsets of even and odd elements of Va, "a 2 A, respectively. A schematic of the three-step process for selecting the active triangle in a simple two-dimensional travel time domain is given in Fig. 3. Suppose that each axis is partitioned into two segments; that is, the domain contains 8(=2  22) triangles, half of which are white and the other half of which are gray as illustrated in the figure. The point P is targeted. It is located in a gray triangle as shown in Fig. 3. One binary variable (=log22) is introduced to select the active column, X1, and another one to select the active row, C1. Binary variable, I1, is added and determines the triangle’s color. First, the active column is selected by setting X1 = 0. Setting C1 = 0 in the second step, the active row is determined that when coupled with the first step column selection reveals the active square. In the third stage, I1 = 0 indicates the gray color of the active triangle and distinguishes it

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

57

for the other white triangle in the active square. Note that the black areas in this figure indicate the union of triangles forbidden to be selected in the process based on the setting of the corresponding binary variables. Upon determination of the active triangle, the targeted point P can be obtained from the unique convex combination of the triangle’s endpoints. This three-step binary branching scheme produces the set of linear constraints (23)–(30) to approximate nonlinear, nonseparable link travel time function (22). 3.2.4. Design decision terms 3 1;h 1;e 3 There are bilinear action level terms c1;e a  ca;s and ca  ca;s ðn Þ in post-response link capacity Eq. (3) and link retrofit cost 2 3 1;e 2 3 Eq. (8), as well as trilinear terms ca  ca;s ðn Þ  ca;s ðn Þ expressed in response action time and cost Eqs. (10) and (11). The bilinear terms are approximated using the LP relaxation of their convex envelopes introduced by McCormick (1976). 3 3 3 1;h 1;e 3 Let first-stage variable u1a;s ¼ c1;e a  ca;s and third-stage variable /a;s ðn Þ ¼ ca  ca;s ðn Þ 8a 2 A; s 2 S. A convex relaxation of first-stage bilinear terms is obtained through a change of variables in (3) and (8) and addition of constraints (31) that give h i 1;h 1;e  ½0; 1 with upper and lower bounds on c1;e an outer-approximation of the rectangular feasible region c1;e a ; ca a and ca;s , respectively. 1;h u1a;s P c1;e a ca;s ; 8a 2 A; s 2 S 1;h 1;e 1;e u1a;s P c1;e a ca;s þ ca  ca ; 8a 2 A; s 2 S

u1a;s 6 ca1;e c1;h a;s ; 8a 2 A; s 2 S

ð31Þ

1;e 1;e u1a;s 6 ca1;e c1;h a;s þ ca  ca ; 8a 2 A; s 2 S 3 1;e 3 1;e Similarly, c1;e a 6 ca 6 ca and 0 6 ca;s ðn Þ 6 1 produces constraints (32).

3 3 /3a;s ðn3 Þ P c1;e a ca;s ðn Þ; 8a 2 A; s 2 S 3 3 1;e 1;e 1;e /3a;s ðn3 Þ P c a ca;s ðn Þ þ ca  ca ; 8a 2 A; s 2 S 3 3 1;e /3a;s ðn3 Þ 6 c a ca;s ðn Þ; 8a 2 A; s 2 S

ð32Þ

3 3 1;e 1;e /3a;s ðn3 Þ 6 c1;e a ca;s ðn Þ þ ca  ca ; 8a 2 A; s 2 S

A generalization of McCormick’s relaxation method proposed by Misener and Floudas (1995) for trilinear terms is applied 2 3 3 3 2 3 2 3 1;e 2 3 to linearize trilinear terms c1;e a  ca;s ðn Þ  ca;s ðn Þ. Let wa;s ðn Þ ¼ ca  ca;s ðn Þ  ca;s ðn Þ, "a 2 A,s 2 S, and replace the trilinear term 2 3 1;e 2 3 1;e in constraints (10) and (11) through a change of variables. Given c1;e a 6 ca 6 ca ; 0 6 ca;s ðn Þ 6 1, and 0 6 ca;s ðn Þ 6 1, addi-

tional constraints (33) are introduced.

Fig. 3. Schematic of three-step process of active triangle selection in Vielma and Nemhauser’s (2011) piecewise linearization method.

58

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

w3a;s ðn3 Þ P 0; 8a 2 A; s 2 S 1;e w3a;s ðn3 Þ P c1;e a  ca ; 8a 2 A; s 2 S 2 3 2 1;e 3 1;e w3a;s ðn3 Þ P c1;e a ca;s ðn Þ þ ca ca;s ðn Þ  ca ; 8a 2 A; s 2 S

ð33Þ

2 3 2 1;e 3 1;e 1;e w3a;s ðn3 Þ P c a ca;s ðn Þ þ ca ca;s ðn Þ  ca ; 8a 2 A; s 2 S 2 2 1;e 3 3 1;e w3a;s ðn3 Þ P c1;e a ca;s ðn Þ þ ca ca;s ðn Þ  ca ; 8a 2 A; s 2 S

1;e 1;e 3 3 1;e w3a;s ðn3 Þ P c1;e a þ ca cas þ ca ca;s ðn Þ  2ca ; 8a 2 A; s 2 S 3.3. Resulting SMILP The original TTRP (2)–(12), a 3-stage SMPEC, is transformed into an equivalent three-stage SMILP:

( maxEn2 jn1

c1 2C1

" max

En3 jn2 ;n1

c2 2C2 ðc1 ;n2 Þ

#) max

hD; ur ðn3 Þi

c3 2C3 ðc2 ;n3 Þ;xr 2Xrx ðn3 Þ

s:t: ðaÞ Three  stage action decision constraints ð3Þ—ð5Þ:

ð34Þ

ðbÞ UE constraints ð20Þ: ðcÞ Link travel time function linear constraints ð23Þ—ð30Þ: ðdÞ LP relaxation of bilinear and trilinear action level terms ð31Þ—ð33Þ: 3.4. Progressive Hedging Algorithm Problem (34) contains pure continuous first- and second-stage variables and mixed integer third-stage variables. Benders-based decomposition methods, e.g. Disjunctive Decomposition-based Branch-And-Cut (D2-BAC) approach by Sen and Sherali (2006), are computationally intensive, and the Lagrangian-based decomposition method by Caroe and Schultz (1999), which might ordinarily be applicable, will not guarantee a globally optimal solution for this problem class. Thus, an exact solution method that is based on concepts of progressive hedging (Rockafellar and Wets, 1991) is presented. This method decomposes the problem by scenario using Lagrangian decomposition. It is particularly attractive here, because it guarantees global optimality for the transformed SMILP (34) with pure continuous first- and second-stage variables within a limited number of iterations. While convexity is not required to achieve a global optimum, global optimality of each subproblem solution is necessary. In this approach, first- and second-stage variables, c1 and c2(n2), are converted into third-stage scenario-dependent variables, c1(n3) and c2(n3), respectively. This allows decomposition of the problem by third-stage information states (i.e. ~1 over scenarios). The following so-called non-anticipativity constraints are added to force c1(n3) to take a single value c 2 3 2 3 2 ~ ðn Þ over those third-stage information all third-stage information states, n , and to force c (n ) to take identical values c states arising from the same n2, i.e. n3jn2.

c1 ðn3 Þ  c~1 ¼ 0; c2 ðn3 Þ  c~2 ðn2 Þ ¼ 0;

ð35Þ ð36Þ

~1 and c ~2 ðn2 Þ are vectors of unrestricted variables. where c The Progressive Hedging Algorithm solves each scenario-specific problem (37) separately wherein non-anticipativity constraints are relaxed.

min

hD; ur ðn3 Þi

c1 2C1 ðn3 Þ;c2 2C2 ðn3 Þ;c3 2C3 ðn3 Þ;xr 2Xrx ðn3 Þ

s:t: ðaÞ three  stage action decision constraints ð3Þ—ð5Þ:

ð37Þ

ðbÞ UE constraints ð20Þ: ðcÞ Link travel time function linearization constraints ð23Þ—ð30Þ: ðdÞ LP relaxation of bilinear and trilinear actionleveltermsð31Þ—ð33Þ: If non-anticipativity constraints (35) and (36) are met, identical solutions for all first- and second-stage variables regardless of the information state n3 will be guaranteed and the problem is solved. However, this is rarely the case. If all first-stage variables are equal, then they are also equal to their expected value. Similarly for second-stage variables. Given the solutions of (37) for all scenarios n3, the expected values of first- and second-stage variables are computed: c^1 ¼ En3 ½c1 ðn3 Þ and c^2 ðn2 Þ ¼ En3 jn2 ½c2 ðn3 Þ, respectively. The deviation in their values from the expected value is measured. ^1 k; kc2 ðn3 Þ  c ^2 ðn2 Þk 6 e, where The optimal solution is obtained when the values converge to the expected value: kc1 ðn3 Þ  c

59

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

e is referred to as the convergence factor. In future iterations, non-anticipativity constraints that do not meet convergence requirements are penalized in the objective function through Lagrangian relaxation as in (38). min

c1 2C1 ðn3 Þ;c2 2C2 ðn3 Þ;c3 2C3 ðn3 Þ;xr 2Xrx ðn3 Þ

þ

q 2

^1 i þ h2 ðn3 Þ; c2 ðn3 Þ  c ^2 ðn2 Þ þ hD; ur ðn3 Þi þ hh1 ðn3 Þ; c1 ðn3 Þ  c

q 2

^ 1 k2 kc1 ðn3 Þ  c

^2 ðn2 Þk2 ; kc2 ðn3 Þ  c

ð38Þ

where h1(n3) and h2(n3) are the vectors of dual variables corresponding to non-anticipativity constraints (35) and (36), and q P 0 is a penalty parameter. The quadratic terms kc1 ðn3 Þ  c^1 k2 and kc2 ðn3 Þ  c^2 ðn2 Þk2 are nonseparable and complicate the     ^1  and c2 ðn3 Þ  c ^2 ðn2 Þ. These absolute solution process. Thus, these terms are replaced by related absolute terms c1 ðn3 Þ  c terms are piecewise linear and can be expressed through the introduction of pairs of continuous nonnegative variable vectors: c1þ ðn3 Þ and c1 ðn3 Þ, and c2þ ðn3 Þ and c2 ðn3 Þ, respectively. Consequently, objective function (38) is replaced by (39), a linear function.

min

c1 2C1 ðn3 Þ;c2 2C2 ðn3 Þ;c3 2C3 ðn3 Þ;xr 2Xrx ðn3 Þ

þ

q 2

^1 i þ hh2 ðn3 Þ; c2 ðn3 Þ  c ^2 ðn2 Þi þ hD; ur ðn3 Þi þ hh1 ðn3 Þ; c1 ðn3 Þ  c

q 2

c1þ ðn3 Þ þ c1 ðn3 Þ



c2þ ðn3 Þ þ c2 ðn3 Þ :

 ð39Þ

Thus, problem (37) is given with its new objective function (39) and the following additional constraints.

c1 ðn3 Þ  c^1 ¼ c1þ ðn3 Þ  c1 ðn3 Þ: c2 ðn3 Þ  c^2 ðn2 Þ ¼ c2þ ðn3 Þ  c2 ðn3 Þ:

ð40Þ ð41Þ

At each future iteration l, the expected values of first- and second-stage variables are updated given the new solutions, and penalties h1(n3) and h2(n3) are revised as in (42) and (43) accordingly.

  ^1 þ h1l ðn3 Þ: h1lþ1 ðn3 Þ ¼ q c1 ðn3 Þ  c   ^2 ðn2 Þ þ h2l ðn3 Þ: h2lþ1 ðn3 Þ ¼ q c2 ðn3 Þ  c

ð42Þ ð43Þ

An overview of the Progressive Hedging Algorithm is depicted in the flowchart of Fig. 4. Global convergence of the proposed Progressive Hedging Algorithm in finite time is assured. Proof given in (Rockafellar and Wets, 1991) for similar two-stage problems with pure continuous first-stage variables can be directly extended to this problem with three stages. 4. Illustrative example The model and solution method are illustrated on a test network with 6 nodes and 16 links representing a highway system as illustrated in Fig. 5. The network was first introduced in (Suwansirikul et al., 1987) and has been used for similar purposes in many works (e.g. Li et al., 2012). Links 2, 5, 6 and 9 represent bridges. Four OD pairs are considered, and the

Fig. 4. Flowchart of the Progressive Hedging Algorithm.

60

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

5 1 7 2

11

9

3

15

4

10

13

6

16 14

12 8

Fig. 5. Test network (Suwansirikul et al., 1987).

Table 2 Origin–destination (OD) pair details. OD ID

Origin

Destination

Travel demand

1 2 3 4

1 6 4 6

6 1 1 2

10 20 5 10

Table 3 Values of link travel time function parameters. Link ID

Link type

t0a

na

ma

coa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Roadway Bridge Roadway Roadway Bridge Bridge Roadway Roadway Bridge Roadway Roadway Roadway Roadway Roadway Roadway Roadway

1 2 3 4 5 2 1 1 2 3 9 4 4 2 5 6

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

10 5 3 20 50 20 10 1 8 3 2 10 25 33 5 1

3 10 9 4 3 2 1 10 45 3 2 6 44 20 1 4.5

corresponding travel demand is reported in Table 2. The network data, including the values of link travel time function parameters from Eq. (22), are given in Table 3, and are identical to values suggested in (Suwansirikul et al., 1987). Three disaster categories (earthquake (s = 1), flood (s = 2) and malicious acts (s = 3)) are considered. A specialized version of Monte Carlo simulation by Chang et al. (1994) is employed to generate disaster scenarios while addressing spatial dependencies that relate to each hazard (see Chen and Miller-Hooks (2012) for more details on the scenario generation process). 30 random scenarios were generated corresponding to each disaster type to capture possible consequences in terms of the level of damage to network links, i.e. 90 scenarios in all were created. All network links except for bridges are candidates for capacity expansion with lower and upper bounds of 0 and 10 units, respectively. The bridge links (2, 5, 6 and 9) can be retrofitted for protection against earthquakes or specific malicious acts. Links 10, 11 and 13–16 at the eastern end of the network may be subject to flooding, and are candidates for related mitigation actions. Second-stage preparedness decisions are considered when flooding is predicted. When the event relates to an earthquake or malicious act, no preparedness actions will be available in this stage. In this case study, a possible secondstage decision may be to pre-position equipment based on information, such as predicted location and event type, received in this stage. For a flooding event, such equipment might include pumps. While these preparedness decisions do not directly impact system performance, they can make the post-event response more effective and efficient. Finally, in stage 3 response actions are considered for restoring capacity of network links that may have been affected under any one of the disaster types. The unit action costs and unit implementation times of response actions are given in Table 4.

61

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64 Table 4 Unit cost of actions. Link ID

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ⁄

Actions Expansion

Retrofit

1;e b a

1;h b a;1

1;h b a;2

1;h b a;3

2 b a;1

2 b a;2

2 b a;3

Response⁄   3 q 3a;1 b a;1

  3 q 3a;2 b a;2

  3 q 3a;3 b a;3

2 – 5 4 – – 4 3 – 5 6 8 5 3 6 1

– 6 – – 8 6 – – 8 – – – – – – –

– – – – – – – – – 3 3.5 – 4 2 4 2

– 2 – – 2 2 – – 2 – – – – – – –

– – – – – – – – – – – – – – – –

– – – – – – – – – 0.5 0.5 – 0.5 0.5 0.5 0.5

– – – – – – – – – – – – – – – –

3.5 (7) 5.5 (11) 8 (16) 7 (14) 5 (10) 5 (10) 4.5 (9) 4 (8) 6 (12) 6.5 (13) 10 (20) 12 (24) 6 (12) 5.5 (11) 7.5 (15) 2.5 (5)

– – – – – – – – – 4.5 (9) 7.5 (15) – 5 (10) 3.5 (7) 5 (10) 1.5 (3)

– 4 (8) – – 3 (6) 3 (6) – – 4 (8) – – – – – – –

Preparedness

3 ðn3 Þ and Note that the perfect preparedness in advance is presumed to reduce response cost (implementation time) by 30%, i.e. b3a;s ðn3 Þ ¼ 0:7b a;s

Resilience

3a;s ðn3 Þ; 8a 2 A; s 2 S q3a;s ðn3 Þ ¼ 0:7q

Fig. 6. Resilience indifference curves for the numerical experiment.

Solution of the SMILP (34) was obtained using the Progressive Hedging Algorithm proposed herein. The algorithm was implemented in GAMS 23.3 and solved using the CPLEX solver. Resilience indifference curves resulting from solution of this illustrative example for different combinations of limited budget B and response time T are provided in Fig. 6. In all problem instances, optimal solutions were obtained within six to ten iterations. Computation times were on the order of three to four hours using a personal computer with Intel Core 2 Duo 2.50 GHz CPU and 3.00 GB RAM running Windows XP Professional Edition. As depicted in this figure, resilience is generally more sensitive to budget than to response time. However, when response times are short, resilience is almost unaffected by budget level. Likewise, when the budget is small, resilience is almost unaffected by a response time increase. In the former case, this is because time available to take action is too restrictive regardless of budget level. In the latter case, funds are unavailable to take additional actions. Detailed results are given through plotting the cumulative distribution of network resilience in Fig. 7 for three sample strategies: (B, T) = (0, 0), (B, T) = (3, 0), and (B, T) = (3, 3). Each point in Fig. 7 represents network resilience under a particular scenario called point resilience. This concept with respect to resilience was introduced in Nair et al. (2010). This figure illustrates that the range and variance of the distribution decreases with larger budget and response time. Moreover, resilience under the worst-case scenario, which occurs at the lowest probability level for each data set, improves with increasing budget and response time. For example, given (B, T) = (0, 0), the range of resilience values is between 0.12 and 1, while with subsequent increases in budget and response time in strategy (B, T) = (3, 3) the range reduces to between 0.56 and 1. Additional runs were conducted assuming a UE could be achieved post-event after the response time has elapsed. Results of these runs are compared in Fig. 8 with those assuming that only a PUE is obtained. From this comparison of results, it is seen that the expected total travel time (the numerator of the resilience measure) is slightly larger under a UE than

62

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

1 B=0, T=0 B=3, T=0 B=3, T=3

Probability

0.8 0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

Resilience Fig. 7. Cumulative distribution function of point resilience.

E[tt]

700

PUE UE

650 600 550

(0,0)

(B,T)

(3,0)

Fig. 8. The network’s total travel time under a UE and PUE.

under a PUE. Through further investigation, it was found that for results from individual sampled scenarios, response actions were focused on different links under the UE and PUE assumptions. In particular, under a PUE the target of the response was on only impacted arcs, while under the UE assumption, unaffected arcs were also improved. 5. Conclusions This paper proposes a novel stochastic network design formulation for maximizing travel time resilience for roadway networks. In particular, it targets freeway networks. The problem explicitly addresses the first three stages of the decision processes of the disaster management life cycle, specifically pre-event mitigation and preparedness, and post-event response. Decisions are taken at each stage based upon the evolution of uncertainty over the stages. The problem is formulated as a bilevel, stochastic mathematical program with partial user equilibrium constraints. The three-stage decision process is embedded within the upper-level problem and user response to the upper-level decisions is modeled in the lowerlevel problem. This study differs from previous studies on stochastic transportation network design in which supply uncertainty is explicitly modeled in the context of long-term mitigation planning applications. In these applications, a UE is an appropriate travel behavior model for estimating network travel times. On the contrary, the PUE approach employed herein accounts for route choice decisions taken by system users assuming that only affected users will have information on the disaster event’s impacts and even these users will have limited information on network damage and completed repairs. A multi-hazard approach is employed and decisions are disaster event-dependent. Thus, mitigation actions may target different hazard scenarios even before the hazard event type is known. In fact, the model accounts for the varying benefits of any such action under different hazard classes. Preparedness decisions are taken only once the hazard class is known, but the specific event realization is uncertain. Response actions are designed for specific disaster event scenarios and are determined once the disaster scenario is realized. This mathematical framework suggests the presence of a single decision-maker with a single budget for investments with both long- and short-term consequences. In reality, multiple agencies could be involved in different aspects of such investments. This tool provides an optimal allocation assuming collaboration among decision-makers. Estimates provided from the proposed techniques, thus, may provide only an upper bound on resilience, as lack of coordination could result in sub-optimal decision-making from the system’s perspective. The developed mathematical model and proposed solution approach can be used to provide benchmark solutions against which solutions from alternative heuristics that can be applied on large-scale problem instances can be evaluated. A number of assumptions are made in conceptualizing this problem. First, problem uncertainties are presumed to be exogenous, where disaster event scenarios are generated independently from the decision process. In this study, the

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

63

mitigation actions, particularly network link retrofit actions, are assumed to only impact disaster consequences in terms of the level of damage to the link and not the disaster probability. The authors are working on extensions in which scenario probabilities are updated based on decisions made at each information state. Secondly, only network supply/capacity uncertainties are taken into account. That is, the travel demand pattern is assumed to be fixed and identical to the original pre-event pattern. One might extend this work to account for demand uncertainty and elasticity. Thirdly, only travel time is considered in modeling user behavior. In reality, other factors, such as safety, might also play a role in user route choice, particularly under disaster events. It is assumed that all affected users are homogenous with respect to the evaluation function (a utility function) used in route selection in obtaining a PUE. Moreover, it is assumed that the users have perfect information on the damaged and repaired network links. Alternative models may be of interest to capture uncertainty in user perception, heterogeneity in user route choice and other factors affecting user decisions. Finally, linear approximations were taken in dealing with nonlinear terms to achieve problem tractability. With improved computer technologies in future years, it may become possible to solve for a global optimum solution directly. Acknowledgments This work was funded by the National Science Foundation. This support is gratefully acknowledged, but implies no endorsement of the findings. References Bureau of Transportation Statistics (BTS), 2013. Pocket Guide to Transportation 2013. Research and Innovative Technology Administration, Bureau of Transportation Statistics, Washington, DC. Carathéodory, C., 1911. Über den Variabilitätsbereich der Fourierschen Konstanten von positiven harmonischen Funktionen. Rendiconti del Circolo Matematico di Palermo 32, 193–217. Caroe, C., Schultz, R., 1999. Dual decomposition in stochastic integer programming. Operations Research Letters 24, 37–45. Chang, C.H., Tung, Y.K., Yang, J.C., 1994. Monte Carlo simulation for correlated variables with marginal distributions. Journal of Hydraulic Engineering 120 (3), 313–331. Chen, L., Miller-Hooks, E., 2012. Resilience: an indicator of recovery capability in intermodal freight transport. Transportation Science 46 (1), 109–123. Chen, A., Zhou, Z., Chootinan, P., Ryu, S., Yang, C., Wong, S.C., 2011. Transport network design problem under uncertainty: a review and new developments. Transport Reviews 31, 743–768. Cox, A., Prager, F., Rose, A., 2011. Transportation security and the role of resilience: a foundation for operational metrics. Transport Policy 18 (2), 307–317. De Palma, A., Rochat, D., 1999. Understanding individual travel decisions: results from a commuters survey in Geneva. Transportation 26, 263–281. Dimitriou, L., Stathopoulous, A., 2008. Reliable stochastic design of road network systems. International Journal Industrial and Systems Engineering 3 (5), 549–574. Fan, Y., Liu, C., 2010. Solving stochastic transportation network protection problems using the progressive hedging-based method. Networks and Spatial Economics 10, 193–208. Faturechi, R., Miller-Hooks, E., 2014. A mathematical framework for quantifying and optimizing protective actions for civil infrastructure systems. Computer-Aided Civil and Infrastructure Engineering Systems 29 (8), 572–589. Faturechi, R., Miller-Hooks, E., 2014. Measuring performance of transportation infrastructure systems in disaster: a comprehensive review. ASCE Journal of Infrastructure Systems, in press. Faturechi, R., Levenberg, E., Miller-Hooks, E., 2014. Evaluating and optimizing resilience of airport pavement networks. Computers and Operations Research 43, 335–348. Gibbons, R., 1992. Game Theory for Applied Economists. Princeton University Press, Princeton, New Jersey. Gordon, P., Kim, S., Moore, J.E., Park, J., Richardson, H.W., 2007. The economic impacts of a terrorist attack on the U.S. commercial aviation system. Commercial Aviation System. Risk Analysis 27, 505–512. He, X., Liu, H.X., 2012. Modeling the day-to-day traffic evolution process after an unexpected network disruption. Transportation Research B 46 (1), 50–71. Heitsch, H., Römisch, W., 2009. Scenario tree modeling for multistage stochastic programs. Mathematical Programming 118 (2), 371–406. Holling, C.S., 1973. Resilience and stability of ecological systems. Annual Review of Ecology and Systematics 4, 1–23. Horn, R.A., Johnson, C.R., 1985. Matrix Analysis. Cambridge University Press, Cambridge. Iida, Y., Kurauchi, F., Shimada, H., 2000. Traffic management system against major earthquakes. IATSS Research 24 (2), 6–17. Larsson, T., Patriksson, M., 1995. An augmented Lagrangean scheme for capacitated traffic assignment problems. Transportation Research B 29 (6), 433–455. Li, C., Yang, H., Zhu, D., Meng, Q., 2012. A global optimization method for continuous network design problems. Transportation Research Part B 46 (9), 1144– 1158. Lo, H., Tung, Y., 2003. Network with degradable links: capacity analysis and design. Transportation Research Part B 37 (1), 345–363. Losada, C., Scaparra, M.P., O’Hanley, J.R., 2012. Optimizing system resilience: a facility protection model with recovery time. European Journal of Operational Research 217, 519–530. Lunderville, N., 2012. Irene Recovery Report, A Stronger Future. Agency of Administration, Office of Secretary, Vermont. McCormick, G.P., 1976. Computability of global solutions to factorable nonconvex programs: Part I – Convex underestimating problems. Mathematical Programming 10, 147–175. Miller-Hooks, E., Zhang, X., Faturechi, R., 2012. Measuring and maximizing resilience of freight transportation networks. Computers and Operations Research 39 (7), 1633–1643. Murray-Tuite, P.M., 2008. Transportation network risk profile for an origin–destination pair: security measures, terrorism, and target and attack method substitution. Transportation Research Record 2041, 19–28. Nagurney, A., Qiang, Q., 2012. Fragile networks: identifying vulnerability and synergies in an uncertain age. International Transactions in Operations Research 19 (1), 123–160. Nair, R., Avetisyan, H., Miller-Hooks, E., 2010. Resilience of ports, terminals and other intermodal components. Transportation Research Record 2166, 54–65. Patriksson, M., 2008. On the applicability and solution of bilevel optimization models in transportation science: a study on the existence, stability and computation of optimal solutions to stochastic mathematical programs with equilibrium constraints. Transportation Research Part B 42, 843–860. Rockafellar, R., Wets, R.B., 1991. Scenarios and policy aggregation in optimization under uncertainty. Mathematics of Operations Research 16, 119–147. Sen, S., Sherali, H.D., 2006. Decomposition with branch-and-cut approaches for two-stage stochastic mixed-integer programming. Mathematical Programming 106 (2), 203–223. Siddiqui, S., Gabriel, S.A., 2013. An SOS1-based approach for solving MPECs with a natural gas market application. Networks and Spatial Economics 13, 205–227.

64

R. Faturechi, E. Miller-Hooks / Transportation Research Part B 70 (2014) 47–64

Sumalee, A., Watling, D., 2008. Partition-based algorithm for estimating transportation network reliability with dependent link failures. Journal of Advanced Transportation 42 (3), 213–238. Sundeen, M., Reed, J.B., 2006. Surface Transportation Funding. Options for States. At the direction of the NCSL Transportation Funding Partnership Committee and NCSL Transportation Committee, National Conference of State Legislatures. Suwansirikul, C., Friesz, T.L., Tobin, R.L., 1987. Equilibrium decomposed optimization: a heuristic for the continuous equilibrium network design problems. Transportation Science 21 (4), 254–263. Todd, M.J., 1977. Union jack triangulations. In: Karamardian, S. (Ed.), Fixed Points: Algorithms and Applications. Academic Press, New York, pp. 315–336. Vielma, J.P., Nemhauser, G.L., 2011. Modeling disjunctive constraints with a logarithmic number of binary variables and constraints. Mathematical Programming, Series B 128, 49–72. Vielma, J.P., Ahmed, S., Nemhauser, G.L., 2010. A note on a superior representation method for piecewise linear functions. INFORMS Journal of Computing 22, 493–497. Vugrin, E., Turnquist, M., Brown, N., 2010. Optimal Recovery Sequencing for Critical Infrastructure Resilience Assessment. Sandia National Laboratories, Report No. SAND2010-6237, Sandia National Laboratories, Albuquerque, MN. Wang, D.Z.W., Lo, H., 2010. Global optimum of the linearized network design problem with equilibrium flows. Transportation Research B 44 (4), 482–492. Waugh, W.L., 2000. Living with Hazards, Dealing with Disasters: An Introduction to Emergency Management. M.E. Sharpe, Armonk, NY. Zhu, S., Levinson, D., Liu, H.X., Harder, K., 2010. The traffic and behavioral effects of the I-35W Mississippi River bridge collapse. Transportation Research Part A 44 (10), 771–784.