Data-driven estimation of interdependencies and restoration of infrastructure systems

Reliability Engineering and System Safety 181 (2019) 167–180

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Data-driven estimation of interdependencies and restoration of infrastructure systems

T

Mauricio Monsalve ,a, Juan Carlos de la Lleraa,b ⁎

a b

Research Center for Integrated Disaster Risk Management (CIGIDEN), CONICYT/FONDAP/15110017, Chile Department of Structural Engineering, Pontifical Catholic University of Chile, Chile

ARTICLE INFO

ABSTRACT

Keywords: Infrastructure systems Interdependent systems Service restoration Resilience

Modern urban systems contain intricate interconnected networks whose components depend on each other to operate, provide value, and sustain a functional society. However, this interconnectedness increases the fragility of these systems by allowing the propagation of disruptions through their interdependencies, which may result in large cascades of failures that can cause severe loss of functionality and recovery capability. Furthermore, the resilience of these systems does not only depend on the individual components, but on their combined ability to recover promptly. With the aim of quantifying the interdependence between these systems, this work introduces a new statistical model for evaluating and simulating the restoration of complex interdependent systems, while modeling their restoration as interdependent processes. The statistical model is introduced along with a custom calibration algorithm that fits the model to observed time series data of infrastructure restoration of functionality. Data from six iconic earthquakes are used to fit and test the model against a suite of service restoration curves associated with different infrastructures. It is concluded that the model may be used to simultaneously estimate the restoration and resilience of an infrastructure system after the disruptions caused by a mainshock. Limitations, possible extensions, and improvements of the model are discussed.

1. Introduction Urban systems and complex networks, such as utilities, transportation, and healthcare, consist of social and technical subsystems embedded in intricate relations of interdependencies [1]. In the case of a natural event, such as an earthquake, these interdependencies allow the propagation of disruptions, increasing the overall fragility of the combined systems [1–5]. Consequently, it is necessary to study such processes to understand how disruptions propagate across systems, and how to better mitigate and reduce their fragility. Ideally, these systems should be studied as large and detailed complex networks of components whose interconnections represent their interdependencies [2–8]; however, this would require considerable amount of data and knowledge about them (see the catalogs of interdependencies of [9–14]), which are still quite limited in practice. Therefore, several researchers have approached this problem previously from a simpler perspective, i.e., aggregating the interdependencies, and leading to a broad collection of approaches for studying cascading failures [15] aimed to bypass detailed modeling of these systems. Simpler methods include network flow methods [16,17], system dynamics models such as stock and flow [18] and influence diagrams [19], simplified weighted network models

⁎

[20], and matrix-based methods [21–25]. In particular, the Inoperability Input-Output model (IIO) [22–25] has become especially popular because of its interpretability and mathematical simplicity. The quantification of interdependencies among systems is a central element in these models, where it is often referred to as coupling strength [26]. However, since the measurement of coupling strength is not part of these models, several researchers have resorted to different methods for their estimation. Expert judgement has been used by several investigators. For instance, Correa et al.proposed eliciting risk maps and scoring [27], Muller also proposed providing estimates using expert judgement to be used in fuzzy logic [28], and The Lifelines Council of the city and county of San Francisco developed a report containing tables and matrices describing the dependencies between the critical infrastructures of the city [29]. Others have relied on quantitative methodologies for estimating the coupling strength between infrastructures. In alignment of the IIO model, Barker and Santos used commercial data (consumption relations) to provide coupling strength estimates [30], Fioriti et al.proposed estimating coupling strength by measuring the synchronization of nodes [31], and Chou and Tseng proposed estimating dependencies from recorded sequences of failures in infrastructure systems [10].

Corresponding author. E-mail addresses: [email protected] (M. Monsalve), [email protected] (J.C. de la Llera).

https://doi.org/10.1016/j.ress.2018.10.005 Received 6 March 2018; Received in revised form 31 July 2018; Accepted 7 October 2018 Available online 09 October 2018 0951-8320/ © 2018 Elsevier Ltd. All rights reserved.

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

Lately, a new literature that addresses the issue of interdependencies during service restoration has become available. In particular, Dueñas-Osorio and Kwasinski proposed a method for quantifying interdependencies between infrastructures based on restoration data, which they applied to study the recovery phase of the 2010 Maule, Chile earthquake [32]. This method has been used and extended by a number of works [33–35]. Independently, Sharkey et al. investigated the restoration interdependencies of the recovery process after hurricane Sandy, using a detailed and extensive approach [14], while Barker and Haimes proposed the IIO model for analyzing infrastructure restoration [24]. Studying interdependencies is especially important because they influence the restoration process and its speed, and have a direct effect on the robustness and resilience of the system. This work introduces a new statistical model for simulating the restoration of an entire infrastructure system after a severe shock (e.g., an earthquake), while at the same time quantifying the dependencies between the infrastructures observed in the restoration data. The model is meant to be calibrated with actual infrastructure restoration data, represented as time series of functionality, or serviceability, of each infrastructure. Afterwards, the model can generate times series of functionality that represent the restoration of the same infrastructure system after a different initial state of damage or disruption. The adequacy of the proposed statistical model is evaluated by its ability to reproduce the simultaneous restoration processes that occurred in several infrastructures after the following six major earthquakes: (i) the 1995 January 17, Mw 6.9, Hanshin-Awaji, Japan earthquake [36]; (ii) the 2004 October 24, Mw 6.6, Mid-Niigata, Japan earthquake [37]; (iii) the 2010 February 27, Mw 8.8, Maule, Chile earthquake [32,38]; (iv) the 2011 February 22, Mw 6.2, Christchurch, New Zealand earthquake [34]; (v) the 2011 March 11, Mw 9.0, Tohoku, Japan earthquake [33]; and (vi) the 2016 April 16, Mw 7.0, Kumamoto, Japan earthquake [39]. The simulation results show that the model is able to reproduce these curves both stochastically (full model) and deterministically (meanfield approximation). This work considerably elaborates on a previous exploratory investigation on deterministic restoration models of interdependent infrastructures [40].

Unlike detailed interdependency data, aggregate data about infrastructure systems is much easier to collect and maintain. Naturally, cascading failures have been studied using coarser data through network flow models [16,17], system dynamics [18,19], aggregate networks [20], and matrix-based methods [21–25]. Of these, the Inoperability Input-Output model (IIO) [22–25], a method adapted from Macroeconomics, has become especially popular. The IIO model is a matrix-based model that attempts to estimate the final size of cascades of failure by using only the initial disruption as input. Let d0 be the column-vector describing the initial disruption to the system, where its ith component represents the initial disruption to the ith infrastructure (component) of the system. Let J be the matrix that quantifies the dependencies of the system, such that the entry in row i and column j represents the loss of functionality that component i suffers when j is inoperable. Then, the initial disruption d0 propagates in one step causing an additional disruption (1)

d1 = J d 0.

Furthermore, consecutive propagations cause disruptions, which at the kth step become

dk = J dk

1

= J k d 0.

(2)

The total disruption c is represented by the sum of all disruptions, i.e.,

c=

J k d0,

dk = k 0

k 0

(3)

which leads to the relation

c = d 0 + J c.

(4)

Thus, as long as the sum converges, the estimated total disruption c to the system is

c = (I

J ) 1d 0.

(5)

Besides estimating total cascade size, the IIO model has also been extended to model system restoration through the inclusion of a resilience matrix K [24]. Nevertheless, the IIO model is limited to representing only linear relationships between the infrastructures. Please note that the above computations also apply to weighted networks; if J corresponds to the weighted adjacency matrix of the network and if the initial disruption is d 0 = (1, …, 1), then the total disruption vector c corresponds to the Katz centrality of the nodes (infrastructures) with discount parameter 1 [42]. The Katz centrality is a social network analysis score that intends to measure the importance of a node within a network; therefore, the IIO model is implicitly performing an analogous to social network analysis but on a system of interdependent infrastructures.

2. Models of interdependent infrastructures 2.1. Cascading failure It has been recognized that interconnected infrastructure systems are susceptible to failures, or disruptions, which propagate through the dependencies between the elements of one or more systems, potentially resulting in large cascading failures [1]. The California power outages in early 2001 are a clear example of a cascading failure: the blackouts affected oil and gas production, refinery operations, jet fuel transport, air travel, and agriculture, just to mention a few impacts [1]. Italy’s nationwide blackout in September 2003 represents another example of a cascading failure, where the mutual dependence between the power and communications infrastructures allowed the disruption to propagate through these systems, until leaving almost the entire country without power and communications [2]. Like these, there are many examples of relatively small disruptions that propagate until becoming of considerable size [1–3,5]. Unfortunately, as human-driven systems become increasingly interconnected and complex, their susceptibility to large-scale cascades of failure increases as well. Since interdependent infrastructure systems can be modeled as networks, researchers have studied the topological properties of these networks to understand cascading failures [4–6,41]. In particular, analytical and simulation studies have shown that interdependent systems exhibit increased fragility, in the sense that the joint fragility of a system is greater than the sum of the fragilities of its components [4,5]. However, prospective studies on realistic networks need access to complete, detailed, and up-to-date data to build the network models, which requires a considerable data collection effort.

2.2. Restoration processes Most research on interdependent infrastructure systems has focused on the fragility of these systems and on the phenomenon of cascading failures, employing a variety of approaches to study these topics [43]. However, despite the early works on infrastructure restoration [44,45], data-centered research is still scant and emergent. In fact, of the types of interdependencies found in the review of Ouyang [43], only one of them, budgetary, is related to restoration or maintenance. The rest of the interdependencies explain how infrastructures depend on each other to provide service or how disruptions propagate from one to another (e.g., geographical interdependency). A number of works have relied on exploratory data analysis (EDA) to measure the coupling strength among infrastructures using service restoration data [32–35]. These works employ time series of functionality as input, which describe the degree of functionality (e.g., percentage) or service that each infrastucture had over the subsequent days after a specific shock (e.g., an earthquake). To exemplify this, let xi(t) ∈ [0, 1] denote the degree of functionality of infrastructure i, t days after the 168

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

shock, with 1 representing full functionality. Then, the time series is subjected to standard EDA treatment, resulting in a treated time series yi(t),

yi (t ) = log xi (t + 2)

2 log xi (t + 1) + log x i (t ),

3.2. Statistical model The proposed model for the restoration process of each infrastructure i from day t to day t + 1 is,

(6)

sij =

| ij |* 1+

|t|*

,

(8)

x i (t + 1) = xi (t ) + ri (t ),

that is generally suitable for exploratory analyses. Then, to quantify the coupling strength between infrastructures i and j, the absolute value of the lagged cross-correlation function between yi and yj is maximized. Let |ρij|* be the maximum cross-correlation obtained, with |t|* the lag which produces it. Then, the coupling strength between infrastructures i and j is quantified as

where ri(t) is the daily restoration rate of infrastructure i on day t, and is defined as

ri (t ) =

i (1

x i (t ))

x j (t ) ij +

i

j

i i (t ),

(9)

were constants αi > 0, βi > 0, and γij ≥ 0 are model parameters; σiϵi(t) is an additive random term, where σi ≥ 0 and ϵi(t) is a standard normal and t denote the infrastructure and random variable; and i time (day) under consideration, respectively. In particular, it is argued that terms γij reflect the degree of dependence of the pair of infrastructures i, j.

(7)

or by other similar measures. Although it is not supported by a model, the above method could be used as an input to, say, the IIO model. In particular, the IIO model has been extended to model the recovery of a system of interdependent infrastructures [24]. In this extension, recovery is modeled as a recurrence relation, where the restoration rate of an infrastructure is a linear combination of the current degree of functionality of all the infrastructures of the system. Similar works include the recovery operator approach of González et al. [46], which states how to find the matrix of interdependencies, or the differential equations of Nguyen et al. [47], which use continuous time. Other works proposed infrastructure restoration estimation methods that are not data-driven. Instead, these methods first require that the analyst defines the interdependencies among the infrastructures. For instance, Lee et al. propose a network flow model for infrastructure system restoration, in which a linear programming model (the flow network) is implied by the types of interdependencies considered [48]. In other cases, the infrastructure restoration models propose restoration strategies to follow [49]. Finally, other works have estimated the recovery times and repair rates of infrastructures neglecting the possible interdependencies between them [50–53]. In these works, survival analysis models (accelerated failure times and proportional hazards) and data mining methods (regression trees) have been fitted to infrastructure restoration data, using as covariates the characteristics of the natural hazard, the damage caused, and the type of affected infrastructure. This work extends a previous exploratory investigation on infrastructure restoration models [40], which already took interdependencies into account, but neglected the stochastic nature of restoration processes and the necessity of defining an adequate calibration algorithm to fit the models.

3.3. Interpretation The model introduced above in Eqs. (8) and (9) embodies the intuition that the restoration rate of an infrastructure does not depend on the restoration rate of the other infrastructures, but on their current level of functionality or serviceability. For example, the restoration rate of the electric distribution depends on the current functionality or serviceability of the road network (travel times, reachability, safety), since it determines the ability to timely reach places where repairs might be needed. Indeed, travel times of the repair teams do not depend on the pace at which roads are being restored, but rather on the current connectivity and quality of the road network. Note that this diverts from the intuition behind the time series cross-correlation models [32–34], which quantify interdependencies from the time lag correlation between twice differenced time series. In Eq. (9), term j x j (t ) ij represents the effect of dependencies on the daily restoration rate of infrastructure i . Herein, terms γij ≥ 0, i, j , reflect the degree or intensity of the dependency of infrastructure i on infrastructure j . If ij = 0, for some pair of infrastructures i, j, then Eq. (9) assumes that infrastructure i is independent from infrastructure j. The other intuition that Eq. (9) embodies is that of diminishing marginal returns. For instance, technicians should prioritize repairs that provide the greatest benefit or impact, say number of clients serviced, at the least cost, time, or difficulty. Hence, as xi(t) becomes closer to 1, its rate of change should slow down. This intuition is encoded in the term (1 x i (t )) i , which decreases as xi approaches 1, and hence ri(t) tends to 0. Constant βi > 0 controls the importance of term (1 x i (t )) i in the rate of change of xi. If βi ≈ 0, then the term (1 x i (t )) i 1 until xi(t) ≈ 1, i.e., when the structure is nearly restored.

3. Model-based interdependence estimation

3.4. Error model

3.1. Formal definitions

The additive error term, σiϵi(t), is composed of a constant standard , and of a random term deviation σi associated with infrastructure i ϵi(t), which follows a standard normal distribution, i.e., ϵi(t) ∼ N(0, 1). It is assumed that ϵi(t) is independent of every i (t ), except for i = i and t = t. Since an additive Gaussian error term is used, the model allows values xi(t) < 0 and xi(t) > 1 in principle, which contradict intuition. Hence, two rules were added to the model to correct this situation: xi(t) ← max (xi(t), 0) and xi(t) ← min (xi(t), 1). These rules are applied in succession. The choice of the error model, which is additive, is inspired by actual restoration data. The restoration curves show that functionality, say the number of customers receiving the service, increases over time, often reproducing a sigmoidal shape. Yet, functionality does not strictly increase over time; it sometimes stalls or decreases as well. The simplest

This work considers the problem of estimating the restoration of a system of interdependent infrastructures after a shock. Let = {0, 1, …, nI 1} be the set of infrastructures under consideration, where nI is their number. Let = {0, 1, …, nT 1} be the time window of restoration in days, where nT corresponds to the last recording (note that nT ). The initial shock to the system, which triggers the restoration process, occurs at time t = 0 . Finally, let xi(t) denote the degree of functionality of infrastructure i at time t . The objective of this research is to propose a model to compute simultaneously the restoration of the system of interdependent infrastructures, i.e., the , t , after an initial state (shock) xi(0), dynamics of xi(t), i i . Moreover, it is also of interest to estimate the degree of in, and the resilience of the terdependency between infrastructures i , j infrastructure system. 169

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

choice that allows ri (t ) = x i (t + 1) x i (t ) to be negative sometimes is the use of an additive Gaussian error. On the multiplicative errors, the simplest is the lognormal, which unfortunately would never allow x i (t + 1) xi (t ) to be negative.

as it would be normally done. Furthermore, the model will be trained so that the observed daily restoration rate approximates the model’s daily restoration rate given by Eq. (9) under the least-squares criterion, i.e.,

3.5. Simulating the restoration after a shock

where ri (t ) = x i (t + 1) x i (t ) is given by Eq. r˜i (t ) = x˜i (t + 1) x˜i (t ) is the observed restoration [ri (t )|xj (t ) = x˜j (t ), j ], i.e., i (t ) =

To reproduce the infrastructure restoration process after a shock, we define the initial state of the system and then apply Eq. (8) for each infrastructure i . The initial state of the system is specified through the values of xi(0), for every infrastructure i at time t = 0 . Then, iterated application of Eq. (8) allows the estimation of xi(1), then xi(2), and so . forth, for all i Since Eq. (8) represents a stochastic process, repeated realizations are necessary for estimating the restoration dynamics of the system of interdependent infrastructures, including the mean daily restoration rate and its variability over time. However, in spite of the above, a mean-field approximation can be , and converting Eq. (8) into estimated easily by enforcing i = 0, i a deterministic model.

i (t ;

( i (t ;

After the initial shock, the restoration of each infrastructure might have different durations. For instance, it is common to see during earthquakes that power distribution is usually restored in a few days while water distribution takes longer to fully recover. In this case, the time window of the dataset should consider the longer duration, and if the model were fitted to the data, power distribution would have several days of complete functionality; thus, a model that just predicts x i (t ) = 1 for the power distribution would have perfect accuracy most of the time, without even making use of the network interdependencies. This example shows that it is necessary to discriminate which data points are useful to fit the model, and which ones are not. The solution to the issue presented above is resolved by weighting the data, as explained next. Let x˜i (t ) denote the observed functionality of infrastructure i on day t. Then, its associated weight wi(t) is determined as:

i

x˜j (t ) ij ,

i

(13)

j

(14)

r˜i (t ))2 .

i , i *)

i , i *)

1 2

=

wi (t )( i (t ;

i,

t

i , i *)

r˜i (t ))2 ,

(15)

The gradient descent algorithm is used to minimize differentiable functions when they are not subject to constraints. Let f(Θ) be the function to minimize, ε ≳ 0 be a stopping threshold, and δ ≳ 0 be a stepsize. Then, the gradient descent algorithm is defined as follows: 1. 2. 3. 4. 5.

Let = 0 be the initial solution. Compute g = f ( ). If ||g|| < ε, then return Θ and terminate; otherwise, g . Update = Go to step 2 and repeat. Thus, minimization of f only requires estimating the gradient ∇f. To fit Si, its gradient ∇Si must be computed for every infrastructure . The derivatives of Si (Eq. (15)) are as follows:

i

4.2. Least-squares criterion

Si

Since the proposed model comprises an additive random Gaussian term, it makes sense to fit the model using the least-squares criterion as suggested by the Gauss–Markov theorem [54]. Under this criterion, the model parameters are chosen so that they minimize the sum of the squared differences between the expected value of the model and the observations (the data). Unfortunately, Eq. (9) does not exhibit the required form of a deterministic term plus a normal error, since the terms xi(t) are all random. To overcome this limitation and stick to the use of least-squares, the approach taken here consists in fitting the conditional probability distribution

=

i

=

i

i (t )

wi (t )( i (t )

r˜i (t ))

wi (t )( i (t )

r˜i (t )) i (t )log(1

wi (t )( i (t )

r˜i (t )) i (t )log(x˜j (t )).

t

Si

,

(16)

i

x˜i (t )),

t

Si

=

ij

t

(17) (18)

The following constraints are also enforced during gradient descent, i

> 0,

i

> 0,

ij

0,

through the application of the following rules,

(10)

instead of fitting the full probability distribution

(xi (t + 1)),

x˜i (t ))

4.3. Gradient descent algorithm

i

)

i (1

for all i . Note that, as explained before, the weights wi(t) have been incorporated in the definition of the objective function Si by weighting the squared differences between ϕi and r˜i .

This weighting scheme only states that, the closer the infrastructure is to full recovery, the less relevant the data becomes for fitting the model.

(xi (t + 1)|xj (t ) = x˜j (t ), j

=

(9) and rate. Let

. Then, minimizing Eq. (12) is equivalent to mimizing

i,

Si ( i,

• if x˜ (t ) 0.95, then w (t ) = 1; else • if x˜ (t ) > 0.95, then w (t ) = 20(1 x˜ (t )). i

i , i *)

(12)

],

The reduction of Eqs. (12)–(14) is supported by the following simple observation. Let h(θ) be a determinist function evaluated on parameter θ, and let ε be a random variable of zero mean whose distribution does not depend on θ. Then, the observation is that the derivatives of [(h ( ) + µ )2] and (h ( ) µ )2 with respect to θ are identical. This follows because of the linearity of the expectation operator, since [(h ( ) + µ )2] [2(h ( ) + µ) leads to h ( )] = 2(h ( ) µ ) h ( ), which in turn is the derivative of (h ( ) µ ) 2 . Following Eq. (14), the least-squares objective function Si to minimize is, then,

4.1. Weighting of observations

i

i,

for all i

4. Model fitting

i

r˜i (t ))2|x j (t ) = x˜j (t ), j

[(ri (t )

(11) 170

i

<

i

= ,

(19)

i

<

i

= ,

(20)

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

ij

<

ij

= 0.

Fig. 1 summarizes the restoration curves for the 1996 Hanshin, 2004 Mid-Niigata, 2011 Christchurch and 2016 Kumamoto earthquakes. The restoration curves for the 2011 Tohoku earthquake are shown in Fig. 2 while the restoration curves for the 2010 Maule earthquake are shown in Fig. 3. The data on the latter two earthquakes are presented separately because not all of the data available in the original works were considered here. Several restoration curves were recorded in the 2011 Tohoku earthquake [33], but only the restoration data that spanned the entire area of interest were considered in this work. Thus, district-level restoration data were not considered herein. Fig. 2(a) shows the restoration curves that are considered in this paper. It is apparent that a sharp loss of functionality occurred on day 28 due to an aftershock that occurred the previous day. To avoid this singularity, the restoration data was cropped until day 26. The result, which is shown in Fig. 2(b), contains the curves that were actually used to test the proposed model. The restoration curves collected for the 2010 Maule earthquake contain information on two different regions of Chile: the Maule and Bio-Bio regions [32]. These two regions are adjacent but administratively independent so that their recovery could be considered as independent. As an example, the distribution of electricity is handled by different companies in these two regions. The collected curves shown in Fig. 3 represent the restoration of the power distribution and the availability of fixed (landlines) and mobile telephony. More restoration curves were available for the Bio-Bio region, but they only regarded the cities of Concepción and Talcahuano (independently) [32], and hence, these data were not used in this work. Some of the empirical restoration curves collected had missing entries, which the model cannot handle directly. The workaround chosen to resolve this issue was to linearly interpolate the time series. Introducing interpolation should not be an problem since only a few values are missing.

(21)

The last rule permits the dropping of potential interdependencies, allowing the model to become simpler. In the task of prediction, models that are less complex, i.e., which consider less variables, tend to outperform complex models [55]. When the algorithm converges, it will be assumed that the parameters αi, βi, γij obtained minimize Si (Eq. (15)). It is at this moment that the σi parameter can be estimated. The standard deviation of the model is simply the weighted average of the squared residuals, i

=

Si ( i, t

i , i *)

wi (t )

.

(22)

5. Experiments 5.1. Evaluation protocol To demonstrate the adequacy of the proposed model (Eq. (8)) for simulating the restoration of interdependent infrastructure systems, the model must demonstrate the ability to reproduce actual restoration data, be competitive against other models in the literature, and be easy to use in practice. Therefore, this section is organized to first present empirical infrastructure restoration curves of six major events, all of which were collected from the literature, as described in Section 5.2. Afterwards, the model was evaluated using these data. More precisely, the stochastic model, its mean-field approximation, and the dynamic IIO model were benchmarked on the collected data. The dynamic IIO model was used as a baseline representing the state of the art. The results of this benchmark are described in Section 5.3. After showing that the proposed model tends to outperform the dynamic IIO model, as thorough discussion of the performance of the stochastic simulations is done in Section 5.4. Furthermore, Section 5.5 contains an analysis of the model residuals and provides justification for the choice of an additive error. Finally, Subsection 5.6 demonstrates how the proposed model can be interpreted and used in practice.

5.3. Performance of the models The proposed model (Eq. (8)) and the IIO model were fitted to the empirical restoration curves using the algorithms introduced earlier. Then, the stochastic simulations, the mean-field approximation, and the dynamic IIO model were compared in the following estimation tasks:

5.2. Infrastructure restoration data Infrastructure restoration curves were collected from academic and technical articles regarding major earthquakes from around the world. Table 1 describes the different earthquakes considered in this work, together with the source of the infrastructure restoration data. The functionalities of the restoration curves, depicted in Figs. 1–3, represent the proportion of served households in the impacted area, except in the case of the Christchurch 2011 earthquake, where functionality is measured either in terms of connectivity or in terms of service provision in the case of hospitals. The references mentioned in Table 1 contain additional information on how these functionality curves were originally collected.

1. Whole curve: the whole restoration curves were simulated from the original (loss of) functionalities measured immediately after each event; the purpose is to evaluate how well the estimation methods reproduce these empirical curves. 2. One day: starting from any day in the empirical data, the functionalities of the next day were estimated. Note that the models were fitted by minimizing the squared errors of these predictions (this is known as training set). 3. Two days: starting from any day in the empirical data, the functionalities of two days away were estimated.

Table 1 Summary of earthquake information and infrastructure restoration data. Name†

Country

Date

Mw

Restoration data

1996 Hanshin (Kobe earthquake) 2004 Mid-Niigata (Chuetsu earthquake) 2010 Maule (Chile 2010 earthquake) 2011 Christchurch (Canterbury earthquake) 2011 Tohoku (Great East Japan earthquake) 2016 Kumamoto

Japan

1996 January 17

6.9

Gas, Power, Water [36]

Japan

2004 October 24

6.6

Gas, Power, Water [37]

Chile

2010 February 27

8.8

New Zealand

2011 February 22

6.2

Japan

2011 March 11

9.0

Fixed and Mobile phones, Power [32] Gas, Power, Water, Hospital, Telecom [34] Gas, Power, Water [33]

Japan

2016 April 16

7.0

Gas, Power, Water [39]

† Alternative names are presented in parenthesis. 171

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin earthquake

(b) 2004 Mid-Niigata earthquake

(c) 2011 Christchurch earthquake

(d) 2016 Kumamoto earthquake

Fig. 1. Infrastructure restoration data on the (a) 1996 Hanshin, (b) 2004 Mid-Niigata, (c) 2011 Christchurch, and (d) 2016 Kumamoto earthquakes.

4. Three days: starting from any day in the empirical data, the functionalities of three days away were estimated.

predicting the restoration curves of the 2010 Maule earthquake, in the Maule region; 2011 Tohoku; 2011 Christchurch; and 2016 Kumamoto. These differences can be explained by the propagation of uncertainties during the recreation of the curves. In fact, the RMSE for the functionality at one, two, or three days shows that the median curve and the mean-field approximation provide very similar estimations. In contrast, the dynamic IIO model tends to produce less accurate estimations in the short term. The functionality at one day, two days, or three days from the data shows that the dynamic IIO yields consistently greater RMSE than the median curve and mean-field approximation of the proposed model, except for the 2004 Mid-Niigata earthquake, where it performs slightly better. In the task of estimating the whole curves, the IIO model is competitive with the proposed model, sometimes providing better RMSE. It is apparent that the dynamic IIO model it is better at fitting the global trends of the data, which may be explained due to its linear structure. Fig. 4 shows the performance of the dynamic IIO model in contrast to the model proposed in this work and the observed restoration curves. In general, the dynamic IIO model can recreate some infrastructure

To measure performance, the root mean squared errors (RMSE) between the observed data and the median of 1000 stochastic simulations, the mean-field approximation, and the dynamic IIO model were computed. Evaluating the RMSE of two or three days after day t (averaged over all days except for the last two or three in each dataset) allows evaluating how well the estimation methods perform after being updated with new data, and how much the estimation error propagates in each step. The RMSE is formally defined as

RMSE =

1 nI nT

it

(x i

x˜i (t ))2 .

(23)

Table 2 shows the RMSE over all infrastructures combined, disaggregated by event, estimation method, and estimation task. The first observation that can be extracted from this Table is that the median curve and the mean-field approximation yield similar results, except for the task of recreating the whole curve, where the RMSE differs in

(a) 2011 Tohoku, original

(b) 2011 Tohoku, cropped

Fig. 2. Infrastructure restoration after the Tohoku earthquake. Left: original data; right: after cropping the curves to day 27. 172

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) Maule region

(b) Bio-Bio region

Fig. 3. Infrastructure restoration data from the 2010 Maule earthquake, for two regions in Chile. Left: restoration in the Maule regions; right: restoration in the BioBio region.

observations near full restoration were weighted down; if x˜i (t ) 0.95, then a weight of 20(1 x˜i (t )) was used, just like with the proposed model. Fitting was performed by minimizing the least squares criterion against the differences x˜i (t + 1) x˜i (t ), i.e., the same as Eq. (14) but adapted to the IIO, and the inequality constraints were implemented using an active set system (pivots, like in SIMPLEX). In the estimation of functionality, it was also enforced that xi(t) ≤ 1, for all i, t. Furthermore, in the tests against the median curves and the mean-field approximations, the dynamic IIO model was not simulated; its linear structure ensures that estimation errors do not propagate non-linearly. This means that the mean-field approximation of the dynamic IIO yields the same results as its median and mean curves. Observe that this approach to fitting the dynamic IIO makes it fit within the recovery operator approach of González et al. [46], but considering constraints.

Table 2 Model performance characterized by root mean squared error, RMSE. In the column labels, Median stands for the median curve of the stochastic simulations, MFA stands for mean-field approximation, and DIIO stands for dynamic IIO. Whole curve Median

Earthquake 1996 2004 2010 2010 2011 2011 2016

Hanshin Mid-Niigata Maule, Maule Maule, Biobio Tohoku Christchurch Kumamoto

Earthquake 1996 Hanshin 2004 Mid-Niigata 2010 Maule, Maule 2010 Maule, Biobio 2011 Tohoku 2011 Christchurch 2016 Kumamoto

One day

0.058 0.009 0.067 0.074 0.089 0.029 0.045 Two days Median 0.012 0.008 0.067 0.057 0.059 0.017 0.017

MFA

DIIO

Median

MFA

DIIO

0.051 0.010 0.069 0.046 0.072 0.039 0.123

0.060 0.007 0.073 0.094 0.052 0.072 0.031

0.012 0.006 0.036 0.047 0.043 0.015 0.014

0.020 0.005 0.077 0.064 0.045 0.042 0.016

MFA 0.012 0.008 0.061 0.057 0.059 0.017 0.016

DIIO 0.028 0.006 0.102 0.084 0.060 0.057 0.024

0.012 0.006 0.040 0.047 0.043 0.015 0.014 Three days Median 0.016 0.009 0.072 0.070 0.062 0.020 0.022

MFA 0.016 0.009 0.066 0.069 0.063 0.020 0.021

DIIO 0.038 0.007 0.080 0.104 0.061 0.054 0.033

5.4. On the stochastic model The stochastic simulations of the model are shown in Fig. 5, while the mean-field simulations are shown in Fig. 6. It is apparent that the simulations of the fitted models, both stochastic and mean-field, resemble the original restoration data, which was already shown in Table 2. Insets (a), (b), (f), and (g) of Fig. 5 and Fig. 6 depict the simulated restoration after the four Japanese earthquakes from the suite of curves considered. The estimations deviate from the data in two cases. First, the restoration curves of the Gas service in the Hanshin earthquake are noticeably different from those of the simulations, being slower than predicted. Indeed, the restoration of the Gas service seems to have started one week after the earthquake, a situation that the proposed model cannot capture. Second, the simulated restoration of the Water distribution in the Tohoku earthquake does not match the maximum reached by the empirical data. The median simulated curves (stochastic and mean-field) reach about 90% in the window considered, while the empirical data reaches 95%. Apart from these two cases, the simulations match the original data closely. Insets (c) and (d) of Fig. 5 and Fig. 6 depict the simulated restoration after the Maule earthquake. In these cases, the start of the original restoration shows a delay, since the restoration of services begins only one or two days after the earthquake. Again, the model is not able to capture such behavior, so it cannot resemble the data with complete accuracy. However, as observed, the simulated curves show similar trends as the restoration curves, especially a few days after the event. Inset (e) of Fig. 5 and Fig. 6 show the simulation results on the Christchurch earthquake. As it can be seen, the simulated curves (stochastic and mean-field) closely follow the data. This might occur as a consequence of having more infrastructure restoration data in a single dataset, which allows a better calibration of the statistical model to reproduce the data.

restoration curves closely. However, it is unable to reproduce the more nonlinear curves that the proposed model approximate better. In some cases, the dynamic IIO model only predicts a straight line, which occurs when no depedendencies are incorporated into the model. Please note that the dynamic IIO model was implemented competitively, so that its comparison against the proposed model is fair. The dynamic IIO model for system restoration can be written as

x i (t + 1)

x i (t ) = a i +

bij xj (t ) + j

i

(t ),

(24)

where ai, bij, and σi are constants, and ϵ(t) is a standard Normal error term. In principle, it should hold that bij ≥ 0, since otherwise an infrastructure would contribute negatively to the restoration of another. However, if the current functionality of an infrastructure is a function of the rest and not of itself,

x i (t + 1) = ai +

bij x j (t ) + j i

i

(t ),

(25)

then convertion of Eq. (25) into Eq. (24)would yield bii = 1, because xi(t) needs to be subtracted from both sides of the equality. Therefore, in fitting (24) to data, it was required that

bii

1

bij

0 (i

j ),

(26)

instead of the more strict non-negativity, bij ≥ 0 for all i, j. No restrictions were imposed on ai to further boost the fitting ability of the IIO model. Also, to further improve the fitting of the IIO model to the data, 173

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin

(b) 2004 Mid-Niigata

(c) 2010 Maule, Maule

(d) 2010 Maule, Biobio

( e ) 2 0 1 1 C h r is t ch u r ch

( f ) 2 0 1 1 To h o k u

(g) 2016 Kumamoto

Fig. 4. Performance of the mean-field approximation (MFA) and the dynamic IIO model (DIIO) on estimating the whole curves.

5.5. Residuals and error model

increases, then it might be proportional to xi(t). Thus, if the error is proportional to the restoration rate, then (xi (t ) x˜i (t ))/ i (t ) should be approximately devoid of trends over time. The residuals are presented in Fig. 7 and do not exhibit a consistent increase over time, which provides evidence against the choice of an error term proportional to the current functionality. However, some of the curves show a pattern of decreasing magnitudes. For instance, the magnitudes of the residuals of the water supply in the Tohoku earthquake (Fig. 7(f)) are evidently decreasing. On the other hand, the magnitudes of the residuals of the gas and power infrastructures do not

The error considered in the proposed model is additive, although other modeling choices are possible. For instance, the error could have been modeled otherwise as multiplicative, in which case, the error could be proportional to the restoration rate, ϕi(t), or, say, to the current functionality, xi(t). An analysis of the residuals, x i (t ) x˜i (t ), can be used to assess which modeling choice is better. An additive error term would be the most adequate choice if the residuals were somewhat independent of ϕi(t) and xi(t). However, if it increases in magnitude as t 174

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin

(b) 2004 Mid-Niigata

(c) 2010 Maule, Maule

(d) 2010 Maule, Biobio

( e ) 2 0 1 1 C h r is t ch u r ch

( f ) 2 0 1 1 To h o k u

(g) 2016 Kumamoto

Fig. 5. Stochastic simulations of the proposed model on the task of estimating the whole curves.

exhibit such trends. The residuals in the Maule earthquake (Fig. 7(c) and (d)) more generally show this decrease in magnitude. Since the decreasing magnitudes of the residuals could be explained by an error proportional to the restoration rate, this hypothesis was also explored. Indeed, the plots depicted in Fig. 8 suggest that modeling an error proportional to the restoration rate might be adequate in some cases (Fig. 8(a), (b), (f)) and inadequate in others (Fig. 8(d), (e), (g)). If the magnitude of the relative residuals, |x i (t ) x˜i (t ))|/ i (t ), did not show an increasing or decreasing trend, then the random term should have been multiplicative instead of additive. However, given the mixed

evidence, it is not obvious that using a different choice of error would have improved the proposed model. Consecuently, the residual plots provide evidence justifying the use of an additive error term instead of a multiplicative one. This does not disregard, however, that the structure of the error could be more complicated. It must also be noted that the fluctuations embedded in the residuals are the results of the epistemic uncertainties associated with the choice of the model, the uncertainties associated with the data collected (number of infrastructures, sampling frequency, resolution), etc, all of which incide on the adequacy of the error term.

175

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin

(b) 2004 Mid-Niigata

(c) 2010 Maule, Maule

(d) 2010 Maule, Biobio

( e ) 2 0 1 1 C h r is t ch u r ch

( f ) 2 0 1 1 To h o k u

(g) 2016 Kumamoto

Fig. 6. Performance of the mean-field approximation (MFA) and the median curves of the stochastic simulations (Median) on estimating the whole curves.

5.6. Example application

x1 (t + 1) = x1 (t ) + 0.3745(1

(28)

To demonstrate a practical use of the model proposed, the model was calibrated for the Maule 2010, Chile earthquake on the Biobio region, and it was evaluated against a series of disruptions. For this region and event, the model is described by the following recurrence relations,

x 0 (t + 1) = x 0 (t ) + 1.5018(1 + 0.0511 0 (t ),

x1 (t ))1.3508x 0 (t )0.4965x2 (t )0.0522 + 0.0234 1 (t ),

x2 (t + 1) = x2 (t ) + 0.7118(1

x2 (t ))1.4985x1 (t )0.4394x2 (t )0.2491 + 0.0621 2 (t ), (29)

where x0(t) represents the functionality of the fixed telephony; x1(t) the functionality of the mobile telephony; and x2(t) the functionality of the power distribution. Eqs. (27)–(29) can be used to simulate the restoration of a system of infrastructures after an arbitrary shock. The following nine scenarios are

x 0 (t ))1.6301x 0 (t )0.4158x1 (t )0.1347x2 (t )0.4803 (27) 176

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin

(b) 2004 Mid-Niigata

(c) 2010 Maule, Maule

(d) 2010 Maule, Biobio

( e ) 2 0 1 1 C h r is t ch u r ch

( f ) 2 0 1 1 To h o k u

(g) 2016 Kumamoto

Fig. 7. Residuals of the predicted restoration of functionality over time.

considered here: (i) all infrastructures drop to 25% functionality; (ii) all infrastructures drop to at 50% functionality; (iii) all infrastructures drop to at 75% functionality; (iv) fixed telephony drops to 25% functionality while the rest to 75%; (v) mobile telephony drops to 25% functionality while the rest to 75%; and (vi) power distribution drops to 25% functionality while the rest to 75%. (vii) fixed telephony drops to 25% functionality while the rest to 100%; (viii) mobile telephony drops to 25% functionality while the rest to 100%; and (ix) power distribution drops to 25% functionality while the rest to 100%. These scenarios are hypothetical and chosen only for illustrative purposes.

Table 3 summarizes the results of the model when applied to the nine scenarios using the mean-field approximation. For each scenario, the following information is evaluated: the time to 95% recovery of infrastructure i, TTRi95% = min{t: x i (t ) 0.95}, for i ∈ {0, 1, 2}; and the sum of functionality losses until 95% recovery, TTR 95% 1

SLFi95% = t = 0 i (1 x i (t )), for i ∈ {0, 1, 2}. Both of these statistics characterize a different aspect of the resilience of the modeled system. The results show that mobile telephony takes longer to recover than power distribution, and that power distribution takes slightly longer to recover than fixed telephony. Averaging over all scenarios, fixed 177

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

(a) 1996 Hanshin

(b) 2004 Mid-Niigata

(c) 2010 Maule, Maule

(d) 2010 Maule, Biobio

( e ) 2 0 1 1 C h r is t ch u r ch

( f ) 2 0 1 1 To h o k u

(g) 2016 Kumamoto

Fig. 8. Relative residuals of the predicted restoration of functionality over time.

telephony takes 4.11 days to recover (TTR095% ), while power distribution takes 6.44 days, and mobile telephony 8.67 days. Furthermore, the average loss of functionality, SLFi95%, grows linearly with the time to recover, TTRi95% (Pearson correlation of 0.914, Spearman correlation of 0.915).

overcomes some issues of previous models by explicitly considering that the restoration rate of each infrastructure depends on the current level of functionality of the other dependent infrastructures rather than their restoration rates. Furthermore, fitting the model also provides an estimate for the coupling strength between the infrastructures modeled. The proposed model considerably extends and elaborates over a previous exploratory investigation on infrastructure restoration models [40]. The model presented is simple, easy to interpret, and easy to evaluate. Furthermore, the model performs well in the task of reproducing

6. Discussion In this work, a new statistical model for the simultaneous estimation of infrastructure restoration curves has been introduced. This model 178

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

incorporating different sources and types of data for calibration. Since these data are essential to any empirical research of infrastructure restoration, it is vitally important that the different systems measure and accurately report the recovery processes. Note also that statistical models cannot discover causality by themselves. Causality is determined using prior knowledge and is confirmed using controlled experiments. Hence, if the model states that, , one cannot conclude that say, ij = 0, for some infrastructures i , j infrastructure i does not depend on infrastructure j. It is true, however, that if the statistical model is right, it will most likely provide the best fit to the observations. Still, an alternative model might have an equally good or even better fit than the true model simply due to chance. Thus, when interpreting results, one must be aware of the shortcomings of statistical models, like the one presented herein.

Table 3 Results of the evaluation of the nine hypothetical scenarios. The formal definitions of TTRi95% (time to recovery) and SLFi95% (sum of functionality losses) are given in the text. Scenario

TTR 095%

TTR195%

TTR295%

SLF095%

SLF195%

SLF295%

(i) (ii) (iii) (iv) (v) (vi) (vii) (viii) (ix)

7 6 4 5 5 5 5 0 0

14 12 9 10 12 9 0 12 0

10 9 7 7 8 9 0 0 8

1.9893 1.1271 0.5321 1.3505 0.6208 0.6942 1.1969 0 0

3.5579 2.2772 1.1369 1.2860 2.9096 1.1593 0 2.7147 0

2.7207 1.6865 0.8570 0.8599 1.0368 2.1450 0 0 1.9365

7. Conclusion

infrastructure restoration processes, which is especially noteworthy considering that the model does not control the propagation of estimation errors in any way. In fact, the fitting algorithm only fits the functionality restoration of an infrastructure to the functionality of the other infrastructures in the same day, completely ignoring data from previous or following days. This suggests that the model proposed is adequate for the amount of information known about the infrastructures considered, and that its performance can serve as a baseline for more sophisticated models. The stochasticity of the model was encoded in the additive random term, σiϵi(t), of Eq. (9). The use of an additive random term supports the use of the Gauss–Markov theorem [54] for calibrating the model. However, since the functional form shown in Eq. (9) is essentially multiplicative, it makes sense to consider the use of a multiplicative random term. But, after inspecting the residuals of the model, no evidence supporting a multiplicative error was found. In fact, either modeling approach could have been used, but reliance on the Gauss–Markov theorem calls for the use of an additive error. Because of its simplicity, however, the model is a black box and thus neglects the internal structure of the infrastructures modeled. In particular, it ignores the network topology, spatial distribution, management, budgetary restrictions, and other features of the infrastructures at hand. The way these dynamics are encoded in the model is through parameters αi, βi, and γij, which are inferred from the functionality data used to train the model. Unfortunately, since these dynamics are not explicitly modeled, it is not possible to extract information about the physical characteristics of the involved infrastructures. To overcome this problem, the model should be extended so that it incorporates additional data. For example, the restoration data can be complemented with network topology data; this might be the subject of future work. Even though the model ignores the technical features of the infrastructures, its functional form encodes some assumptions regarding its management and budgetary restrictions during recovery. It was assumed that, during infrastructure restoration works, repairs that were the most effective, i.e., had the best restoration-effort relation, were prioritized over these that demanded more resources (personnel, materials, equipment) or provided less benefit to customers. This was the purpose of the (1 x i (t )) i term in Eq. (9); the term diminishes the amount of restored functionality as less functionality is there to restore, following the assumption that the least effective repairs are performed last. Another expense that affects restoration time is periodic maintenance. However, maintenance is not associated with restoration activities and is, therefore, ignored by the model, although it might condition parameters αi, βi, and γii for each infrastructure i. Application of the proposed model in actual decision-making requires collecting data on more infrastructures and with better granularity than the one found usually in the literature. This severely limits the applicability of the proposed model in research as well, except for the purposes of model development (rather than the study of the interdependent infrastructures) and improvements that consider

Urban and complex engineering systems, such as utilities, transportation, and healthcare, consist of social and technical subsystems embedded in intricate relationships that present interdependencies, which allow the propagation of disruptions, which, in turn, may result in large cascades of failures. This phenomenon of heightened fragility motivates the study of the extent of interdependence between infrastructures and how these interdependencies condition fragility and recovery. This work introduced a new statistical model aimed to quantify the interdependencies between different infrastructures and simulate the restoration process of a complete infrastructure system after a shock, such as a severe earthquake. A calibration algorithm that fits the model to observed time series of service or infrastructure restoration was introduced as well, and the model was tested against a selection of six iconic earthquakes whose restoration data was available in the literature. The model was capable of reproducing the concurrent recovery of different utility networks, which enables decision makers to use it to predict responses in different initial conditions. Estimated mean root square errors are usually less than 7% for the mean-field process. While the introduced model is simple, easy to interpret, and easy to evaluate, it is subject to several limitations. The model is subject to the typical limitations that statistical models have, especially with causality. The data the model requires to calibrate, which are the different service or infrastructure restoration time series, are difficult to find in practice and might suffer from data quality issues. The model also neglects the network topology of the infrastructures considered, as well as their spatial distribution. Overcoming some of these limitations will be the subject of future research work. Acknowledgments The authors kindly express their gratitude towards the following sponsorship from the government of Chile: CONICYT/FONDECYT/ 3170867 (FONDECYT Postdoctorado), CONICYT/FONDECYT/1170836 (FONDECYT Regular), and CONICYT/FONDAP/15110017 (Centro de Investigación para la Gestión Integrada de Desastres Naturales, CIGIDEN). References [1] Rinaldi SM, Peerenboom JP, Kelly TK. Identifying, understanding, and analyzing critical infrastructure interdependencies. IEEE Control Syst 2001;21(6):11–25. [2] Buldyrev SV, Parshani R, Paul G, Stanley HE, Havlin S. Catastrophic cascade of failures in interdependent networks. Nature 2010;464(7291):1025–8. [3] Brummitt CD, D’Souza RM, Leicht EA. Suppressing cascades of load in interdependent networks. Proc Natl Acad Sci 2012;109(12):E680–9. [4] Bashan A, Berezin Y, Buldyrev SV, Havlin S. The extreme vulnerability of interdependent spatially embedded networks. Nat Phys 2013;9(10):667–72. [5] Daqing L, Yinan J, Rui K, Havlin S. Spatial correlation analysis of cascading failures: congestions and blackouts. Sci Rep 2014;4:5381. [6] Winkler J, Dueñas-Osorio L, Stein R, Subramanian D. Interface network models for

179

Reliability Engineering and System Safety 181 (2019) 167–180

M. Monsalve, J.C. de la Llera

IEEE; 2010. p. 85–7. [32] Dueñas-Osorio L, Kwasinski A. Quantification of lifeline system interdependencies after the 27 February 2010 Mw 8.8 offshore Maule, Chile, earthquake. Earthquake Spectra 2012;28(S1):S581–603. [33] Cimellaro GP, Solari D, Bruneau M. Physical infrastructure interdependency and regional resilience index after the 2011 Tohoku earthquake in Japan. Earthquake Eng Struct Dyn 2014;43(12):1763–84. [34] Zorn CR, Shamseldin AY. Quantifying directional dependencies from infrastructure restoration data. Earthquake Spectra 2016;32(3):1363–81. [35] Krishnamurthy V, Kwasinski A, Dueñas-Osorio L. Comparison of power and telecommunications dependencies and interdependencies in the 2011 Tohoku and 2010 Maule earthquakes. J Infrastruct Syst 2016;22(3):04016013. [36] Nojima N. Restoration processes of utility lifelines in the great east Japan earthquake disaster, 2011. 15th world conference on earthquake engineering (15 WCEE). 2012. Lisbon [37] Kajitani Y, Sagai S. Modelling the interdependencies of critical infrastructures during natural disasters: a case of supply, communication and transportation infrastructures. Int J Crit Infrastruct 2009;5(1–2):38–50. [38] de la Llera JC, Rivera F, Mitrani-Reiser J, Jünemann R, Fortuño C, Ríos M, et al. Data collection after the 2010 Maule earthquake in Chile. Bull Earthquake Eng 2017;15(2):555–88. [39] Nobuoto Nojima YM. Comparison of functional damage and restoration processes of utility lifelines in the 2016 Kumamoto earthquake, Japan with two great earthquake disasters in 1995 and 2011. Tech. Rep. FS2016-L-0005. JSCE Disaster Reports; 2016. http://committees.jsce.or.jp/disaster/FS2016-L-0005 [40] Monsalve M, de la Llera JC. New models for estimating interdependencies and recovery of infrastructure systems. 16th European conference on earthquake engineering (16ECEE). 2018. Thessaloniki, Greece [41] Hernandez-Fajardo I, Dueñas-Osorio L. Probabilistic study of cascading failures in complex interdependent lifeline systems. Reliab Eng Syst Saf 2013;111:260–72. [42] Borgatti SP, Everett MG. A graph-theoretic perspective on centrality. Soc Netw 2006;28(4):466–84. [43] Ouyang M. Review on modeling and simulation of interdependent critical infrastructure systems. Reliab Eng Syst Saf 2014;121:43–60. [44] Isumi M, Nomura N, Shibuya T. Simulation of post-earthquake restoration for lifeline systems. Int J Mass Emergencies Disasters 1985. [45] Kozin F, Zhou H. System study of urban response and reconstruction due to earthquake. J Eng Mech 1990;116(9):1959–72. [46] es D. González A, Chapman A, nas Osorio LD, Mesbahi M, D’Souza RM. Efficient infrastructure restoration strategies using the recovery operator. Comput Aided Civil Infrastruct Eng 2017;32:991–1006. [47] Nguyen T-D, Cai X, Ouyang Y, Housh M. Modelling infrastructure interdependencies, resiliency and sustainability. Int J Crit Infrastruct 2016;12(1/2). [48] Lee EE, Mitchell JE, Wallace WA. Restoration of services in interdependent infrastructure systems: a network flows approach. IEEE Trans Man Syst Cybern Part C 2007;37(6):1303–17. [49] es D. González A, nas Osorio LD, Sánchez-Silva M, Medaglia AL. The interdependent network design problem for optimal infrastructure system restoration. Comput Aided Civil Infrastruct Eng 2016;32:334–50. [50] Liu H, Davidson RA, Apanasovich TV. Statistical forecasting of electric power restoration times in hurricanes and ice storms. IEEE Trans Power Syst 2007;22(4):2270–9. [51] Nateghi R, Guikema SD, Quiring SM. Comparison and validation of statistical methods for predicting power outage durations in the event of hurricanes. Risk Anal 2011;31(12):1897–906. [52] Barker K, Baroud H. Proportional hazards models of infrastructure system recovery. Reliab Eng Syst Saf 2014;124:201–6. [53] Zorn CR, Shamseldin AY. Post-disaster infrastructure restoration: a comparison of events for future planning. Int J Disaster Risk Reduct 2015;13:158–66. [54] Gauss–Markov theorem. The concise encyclopedia of statistics. New York, NY: Springer New York978-0-387-32833-1; 2008. p. 217–8. [55] Arlot S, Celisse A. A survey of cross-validation procedures for model selection. Stat Surv 2010;4:40–79.

complex urban infrastructure systems. J Infrastruct Syst 2011;17(4):138–50. [7] Borge-Holthoefer J, nos RAB, Gonzalez-Bailon S, Moreno Y. Cascading behavior in complex socio-technical networks. J Complex Netw 2013(1):3–24. [8] Kivelä M, Arenas A, Barthelemy M, Gleeson JP, Moreno Y, Porter MA. Multilayer networks. J Complex Netw 2014;2(3):203–71. [9] McDaniels T, Chang S, Peterson K, Mikawoz J, Reed D. Empirical framework for characterizing infrastructure failure interdependencies. J Infrastruct Syst 2007;13(3):175–84. [10] Chou C-C, Tseng S-M. Collection and analysis of critical infrastructure interdependency relationships. J Comput Civil Eng 2010;6(24):539–47. [11] Chang SE, McDaniels TL, Mikawoz J, Peterson K. Infrastructure failure interdependencies in extreme events: power outage consequences in the 1998 ice storm. Nat Hazards 2007;41(2):337–58. [12] Shoji G, Toyota A. Modeling of restoration process associated with critical infrastructure and its interdependency due to a seismic disaster. TCLEE 2009: lifeline earthquake engineering in a multihazard environment. 2009. p. 1–12. [13] Di Giorgio A, Liberati F. A Bayesian network-based approach to the critical infrastructure interdependencies analysis. IEEE Syst J 2012;6(3):510–9. [14] Sharkey TC, Nurre SG, Nguyen H, Chow JH, Mitchell JE, Wallace WA. Identification and classification of restoration interdependencies in the wake of hurricane sandy. J Infrastruct Syst 2016;1(22):04015007. [15] Cimellaro GP. A comprehensive methodology for the evaluation of infrastructure interdependencies. Urban resilience for emergency response and recovery. Springer; 2016. p. 139–223. [16] Holden R, Val DV, Burkhard R, Nodwell S. A network flow model for interdependent infrastructures at the local scale. Saf Sci 2013;53:51–60. [17] Nicholson CD, Barker K, Ramirez-Marquez JE. Flow-based vulnerability measures for network component importance: experimentation with preparedness planning. Reliab Eng Syst Saf 2016;145:62–73. [18] Gillespie DF, Robards KJ, Cho S. Designing safe systems: using system dynamics to understand complexity. Nat Hazard Rev 2004;5(2):82–8. [19] Andrijcic E, Haimes YY. Metamodeling of interdependent systems: application to bridge infrastructure management. J Infrastruct Syst 2016;23(2):04016028. [20] Pinnaka S, Yarlagadda R, Çetinkaya EK. Modelling robustness of critical infrastructure networks. 2015 11th international conference on the design of reliable communication networks (DRCN). IEEE; 2015. p. 95–8. [21] Tsuruta M, Goto Y, Shoji Y, Kataoka S. Damage propagation caused by interdependency among critical infrastructures. The 14th world conference on earthquake engineering (14WCEE). 2008. [22] Haimes YY, Jiang P. Leontief-based model of risk in complex interconnected infrastructures. J Infrastruct Syst 2001;7(1):1–12. [23] Haimes YY, Horowitz BM, Lambert JH, Santos JR, Lian C, Crowther KG. Inoperability input-output model for interdependent infrastructure sectors. i: theory and methodology. J Infrastruct Syst 2005;11(2):67–79. [24] Barker K, Haimes YY. Uncertainty analysis of interdependencies in dynamic infrastructure recovery: applications in risk-based decision making. J Infrastruct Syst 2009;15(4):394–405. [25] MacKenzie CA, Barker K. Empirical data and regression analysis for estimation of infrastructure resilience with application to electric power outages. J Infrastruct Syst 2012;19(1):25–35. [26] Liu R-R, Li M, Jia C-X. Cascading failures in coupled networks: the critical role of node-coupling strength across networks. Sci Rep 2016;6. [27] Correa-Henao GJ, Yusta JM, Lacal-Arántegui R. Using interconnected risk maps to assess the threats faced by electricity infrastructures. Int J Crit Infrastruct Prot 2013;6(3):197–216. [28] Muller G. Fuzzy architecture assessment for critical infrastructure resilience. Procedia Comput Sci 2012;12:367–72. [29] of the City LC, of San Francisco C. City and county of san francisco lifelines interdependency study. Tech. Rep.. 2014. [30] Barker K, Santos JR. A risk-based approach for identifying key economic and infrastructure systems. Risk Anal 2010;30(6):962–74. [31] Fioriti V, D’Agostino G, Bologna S. On modeling and measuring inter-dependencies among critical infrastructures. Complexity in engineering 2010 (COMPENG’10).

180