Review on Energy Resilience

Reliability Engineering and System Safety 154 (2016) 106–116

Contents lists available at ScienceDirect

Reliability Engineering and System Safety

http://d 0951-83

E-m

journal homepage: www.elsevier.com/locate/ress

Critical location identification and vulnerability analysis of interdependent infrastructure systems under spatially localized attacks

Min Ouyang a,b

a Department of Systems Science and Engineering, Huazhong University of Science and Technology, 1037 Luoyu Road, Wuhan 430074, PR China b Key Lab. for Image Processing and Intelligent Control, Huazhong University of Science and Technology, Wuhan 430074, PR China

a r t i c l e i n f o

Article history: Received 3 October 2015 Received in revised form 5 May 2016 Accepted 13 May 2016 Available online 20 May 2016

Keywords: Infrastructure systems Interdependencies Vulnerability analysis Spatially localized attacks Critical location

x.doi.org/10.1016/j.ress.2016.05.007 20/& 2016 Elsevier Ltd. All rights reserved.

ail address: [email protected]

a b s t r a c t

Infrastructure systems are usually spatially distributed in a wide area and are subject to many types of hazards. For each type of hazards, modeling their direct impact on infrastructure components and analyzing their induced system-level vulnerability are important for identifying mitigation strategies. This paper mainly studies spatially localized attacks that a set of infrastructure components located within or crossing a circle shaped spatially localized area is subject to damage while other components do not directly fail. For this type of attacks, taking interdependent power and gas systems in Harris County, Texas, USA as an example, this paper proposes an approach to exactly identify critical locations in interdependent infrastructure systems and make pertinent vulnerability analysis. Results show that (a) infrastructure interdependencies and attack radius largely affect the position of critical locations; (b) spatially localized attacks cause less vulnerability than equivalent random failures; (c) in most values of attack radius critical locations identified by considering only node failures do not change when con- sidering both node and edge failures in the attack area; (d) for many values of attack radius critical locations identified by topology-based model are also critical from the flow-based perspective.

& 2016 Elsevier Ltd. All rights reserved.

1. Introduction

Infrastructure systems, such as electric power, water supply, telecommunication and emergence service systems, form the backbone for functioning of a community or nation and provide essential services to support the well-being of its citizens. How- ever, these systems are currently subject to many types of hazards, such as unavoidable natural hazards, component aging, sharp demand increase and climate change, which increase their failure probabilities and vulnerability. Also, infrastructure systems are not isolated but highly interconnected and mutually interdependent [31,33]. Interdependencies can improve infrastructure operational efficiency, but they may also increase system vulnerability, i.e., small failures in one system can result in cascading failures within it and across to other systems, causing large economic loss and affecting people’s life in the disaster area. Hence, to better protect interdependent infrastructures, for each type of hazards, it re- quires modeling their direct impact on infrastructure components and then analyzing their system-level vulnerability for identifying mitigation strategies, where the term “vulnerability” in this paper is interpreted to associate with a disruptive event and is quantified as the performance drop of a system under an event [44].

To model and analyze the vulnerability of interdependent in- frastructure systems, there exist many approaches in the literature, such as empirical approaches, agent based approaches, system dynamics based approaches, economic theory based approaches, network based approaches and others. A detailed review of these approaches is provided by the author in the reference [24]. This paper addresses the problem of identifying critical locations in infrastructure systems, which needs their topological and geo- graphical information, so a network-based approach is used in this paper for their vulnerability modeling. This paper next will only review some recent network-based vulnerability studies on in- terdependent infrastructure systems. This type of approaches models each system by a network and describes their inter- dependencies by inter-links, which enable to capture system to- pological features and flow characteristics for identifying critical system components [10,16,5].

Depending on the failure modes, network-based vulnerability studies on interdependent systems can be grouped into three types. The first is vulnerability studies under random failures, which are usually modeled by randomly removing a certain frac- tion or number of infrastructure components [15,29,3,34,36–38,6]. The second is vulnerability studies under natural hazards, such as earthquakes, hurricanes and lightening, and their impacts on in- frastructure components are usually modeled according to fragility

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116 107

curves, which provide the probability of exceeding a certain da- mage state threshold conditional to a selected hazard intensity measure, such as peak ground acceleration or peak ground velo- city or permanent ground deformation for seismic hazards [1,11,14,32,8], 3 second gust wind speed for hurricane hazards [13,22,26]. The third is vulnerability studies under intentional at- tacks, which can be further grouped into two types, depending on whether component geographical proximity is considered by the attack strategies.

The first type of intentional attacks is non-proximity-based attacks, where a set of important components (nodes or edges) is selected as the attack objective without considering their geo- graphical proximity. This type of attacks has been studied in two ways in the literature. One is random networks based studies where only node degree distributions of networks of concern are used for analysis and the largest component size is taken as the performance metric to quantify system vulnerability under the removal of a fraction of critical nodes selected according to their degree. This type of studies can be analytically solved by gen- erating function method when the network size tends to infinity, and results show that interdependencies increase system vulner- ability in contrast to independent scenarios, while interdependent networks are difficult to defend by strategies such as protecting the high degree nodes that are found useful to significantly reduce the vulnerability of single networks [18]. The other is specific networks based studies, where system topological information and sometimes component physical properties and flow-related parameters are used for analysis. This type of studies is usually analyzed by Monte-Carlo simulation technique or advanced algo- rithms, and usually takes more practical system operation models and performance metrics to quantify system vulnerability. In this case, the critical components can be selected not only by initial degree and recalculated degree [12,42], but also by component load levels [38], or by some optimization techniques or advanced algorithms to select components whose failures can cause the largest vulnerability [27,39,40,43].

The second type of intentional attacks is proximity-based at- tacks that a set of components (nodes or edges or both) distributed in close proximity to each other in a localized area is selected as the attack objective. This type of attacks can simulate bomb at- tacks and military weapons with mass destruction. For inter- dependent randomly diluted square lattices with average node degree between 2.5 and 4.0, and the dependencies between nodes with and across networks are constrained less than a certain dis- tance, Berezin et al. [4] modeled localized attacks by removing all nodes within a certain distance from a random attack center, and found that localized attacks can be more harmful than random attacks with the same intensities. Note that this purely topology- based study only takes the largest component size as the perfor- mance metric to quantify system vulnerability, which has weak correlations to those results obtained when the infrastructure system flow properties are considered [17,23,25,9]. Also, this study only considers node failures within the attack influence area and the failures of edges that cross the area are not taken into account, but a recent work from the author shows that considering edge failures in the attack area will affect the position of critical loca- tions [28]. For specific infrastructure networks, some scholars have considered their vulnerability by removing both nodes and edges co-located in the same attack influence area, which enables to identify the critical area so that the failures of components inside such an area lead to the largest vulnerability. Patterson and Apostolakis [30] made such a vulnerability analysis on inter- dependent systems in a university campus area by dividing the infrastructure map into a generic hexagonal grid with a small ra- dius 7 m and each failure or attack influence area is modeled by a hexagon. Johansson and Hassel [19] have made a similar analysis

to identify critical locations in interdependent systems by con- sidering square attack influence areas with two different sizes 5 � 5 km2 and 2.5 � 2.5 km2. Note that the results from these studies depend on the shape of attack influence area (square, hexagon, etc.) and the attack influence area size as well as how the infrastructure map is partitioned into small attack influence areas.

This paper studies a type of spatially localized attacks that a set of infrastructure components located within or crossing a circle shaped spatially localized area is subject to damage while other components outside the area do not directly fail, which belongs to the second type of intentional attacks. For this type of attacks, this paper proposes an algorithm to identify critical locations in in- terdependent systems, where node and edge failures are both considered in the attack areas and infrastructure flow properties are also included for analysis. Note that different with previous research on identifying critical components or locations in infra- structure networks, which do not consider the spatial proximity of components in the critical component set [2,20,35], this paper identifies the critical components constrained to be all located within a limited circle shaped localized area. Also, different with the works from Patterson and Apostolakis [30] and Johansson and Hassel [19], which identify the critical components or locations based on the partition of infrastructure map, this paper does not partition the map and uses a proposed algorithm to exactly find out the critical locations. The rest of the paper is organized as follows: Section 2 provides a network-based vulnerability model of interdependent infrastructure systems. Section 3 introduces the model of spatially localized attacks and the approach to identify critical locations. Section 4 takes interdependent power and gas systems in Harris County, USA as an example to identify critical locations and make pertinent vulnerability analysis. Section 5 discusses the findings and extensions and Section 6 provides conclusions and future research.

2. Modeling vulnerability of interdependent infrastructure systems

This paper defines vulnerability of an infrastructure system under a specific disruptive event as its performance drop, then the critical location is an area that failures of components located within and crossing this area can make the largest vulnerability. Hence, to identify critical locations, it first needs a vulnerability model of interdependent systems, which can simulate the cas- cading failures within and between multiple systems and enable to compute system vulnerability under a particular disruption scenario. Define system performance in the normal state by PMn

and in the post-event state by PMd, then system vulnerability under a disruptive event is computed as follow:

= −

( )V PM PM

PM 1 n d

n

The disruptive event produces initial failure scenarios (failure states) of infrastructure systems components, and then to com- pute system vulnerability under this event, it needs to simulate the cascading failures within and between multiple systems. This paper next takes interdependent power and gas systems in Harris County, Texas as an example to illustrate how to describe them by networks and model their vulnerability under a specific disruption scenario.

For the gas system in Harris County, the gas compressors, the gas storage facilities, the gas delivery facilities, the gas receipt fa- cilities and the gas pipeline junctions are modeled as nodes while the gas pipeline segments are described as edges, with their re- presentation in the reference [22]. In addition, the gas storage

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116108

facilities and the gas receipt facilities are regarded as source nodes, the connection points and gas compressors are modeled as transmission nodes and the gas delivery facilities as load nodes. Finally, there are EG¼67 links and NG¼63 nodes in total, with NG,S¼14 source nodes, NG,T¼23 transmission nodes and NG,L¼26 load nodes. The capacity of each pipeline segment is set as c�Dr with r¼2.5 [7].

For the power system in Harris County, it includes the high voltage (35–345 kv) transmission network and the low voltage (0.12–35 kv) distribution network. For the distribution network, its information is not available for analysis due to security concerns and then only the transmission network is used for analysis. The transmission network has EP¼551 transmission lines and NP¼417 nodes, including NP,S¼23 power plants and 394 substations, with its geographical representation shown in the reference [22]. All transmission substations are initially regarded as load nodes, with their load levels estimated proportionally to the number of households in that tract.

Note that there exist four types of interdependencies among critical infrastructures: physical, cyber, geographic and logical [24,33]. Two critical infrastructures are regarded interdependent if it exist at least one type of these four interdependencies between them. The Harris power and gas systems of concern are located in the same region and the spatially localized attacks studied in this paper can cause the state change of both two systems simulta- neously, which means it exist geographical interdependencies between them and the two systems in the paper are inter- dependent. For the physical relationships between power and gas systems, for illustrative purpose, this paper only considers uni- directional interdependencies from the power grid to the gas system for simplification. This paper assumes that all types of gas nodes require electricity to keep their normal operation, and each gas node connects to its nearest transmission load substation. Fi- nally, it has 63 inter-links, which describe the physical relation- ships from the power transmission network to the gas system for power supply. Note that when considering the bi-directional physical interdependencies, the proposed approach is still valid.

Based on the above network-based description, this paper uses two types of models for vulnerability assessment. When only in- frastructure topological information is used for analysis, a topol- ogy-based performance model is used to simulate the cascading failures according to the following procedures: First, identify power system component failures directly from the event, and then check system topology to record all isolated sub-networks (a disruptive event may disconnect the original power network into several isolated sub-networks, where nodes in each sub-network are reachable from each other but are not reachable from nodes in other sub-networks). For each sub-network, if it includes no generator node, then all nodes in this sub-network are failed; if it includes at least one generator node, it assumes all load nodes in this sub-network can still supply power to the gas system, but the performance of each load node i is compromised and is quantified as the number of post-event connected generator nodes NP S i

d , , di-

vided by the total number of connected generator nodes NP, S, i in the normal state, i.e., N N/P S i

d P S, , , , then the topology-based post-

event power system performance PMT P d , that is the total sum of the

performance of each power load node i is computed as follow:

∑= ( )=

PM N

N 2 T P d

i

N P S i d

P S ,

1

, ,

,

P L,

Second, identify gas system component failures directly from the disruptive event or indirectly due to the failures of their de- pended power load nodes or the failures of inter-links, then em- ploy the similar connectivity model for the power system to es- timate topology-based post-event gas system performance PMT G

d ,

that is the total sum of the performance of each gas load node as follow:

∑= ( )=

PM N

N 3 T G d

i

N G S i d

G S i ,

1

, ,

, ,

G L,

Where NG S i d

, , is the number of post-event connected source nodes

NG S i d

, , for gas load node i, NG,S, i is the total number of connected

source nodes in the normal state, N N/G S i d

G S i, , , , is the post-event performance of gas load node i. Note that if there have bi-direc- tional interdependencies between two systems, repeating the above two steps can compute the post-event steady-state system performance.

When infrastructure flow properties are further considered for analysis, a flow-based performance model is used to simulate the cascading failures according to the following procedures: First, identify power system component failures directly from the event, and then check the system topology to record all unconnected sub-networks. For each sub-network, its response is modeled ac- cording to the following rules: (1) If it does not contain any power plant, then all load nodes are assumed failed (in terms of services); (2) If the sum of all power plant capacities is larger than the sum of the demand, then check the line flow constraints after running a DC power flow equation; if there are violations, cut the load to the load node with the smallest load level until the line constraints are satisfied; (3) If the sum of all power plant capacities is smaller than the sum of the demand, cut the load to the load node with the smallest load level first to balance the supply and the demand, and then run the step (2). The DC power flow equation F¼AP provides a relationship between power flow vector F consisting of the flow through each line and the power injection vector P at the nodes, with A as a constant matrix. When line flow constraints are sa- tisfied, power grid component states and performance level can then be recorded and computed. Define the load level of a node i in normal state by Pi

n and in post-event steady state by Pi d, then the

flow-based post-event power system performance PMF P d , that is

the ratio of post-event total power supply ∑ = Pi N

i d

1 P L, to pre-event

total power supply ∑ = Pi N

i n

1 P L, is computed as follow:

∑ ∑= ( )= =

PM P P/ 4

F P d

i

N

i d

i

N

i n

, 1 1

P L P L, ,

Second, identify gas system component failures directly from the disruptive event or indirectly due to power system outage or the failures of inter-links, then employ the network maximum flow model [7] to estimate gas system performance. Define the maximum flow from the source nodes to the load nodes in the normal state by MFn and in the post-event state by MFd, then the flow-based post-event gas system performance PMF G

d , that is the

ratio of post-event maximum flow MFd to pre-event maximum flow MFn is computed as follow:

= ( )PM MF MF/ 5F G d ,

d n

Similarly, if there exist bi-directional physical inter- dependencies, repeating the above two steps can compute their post-event system performance.

Based on the above performance metrics in Eqs. (2)–(5) and the vulnerability metric in Eq. (1), when system topology-based or flow-based performance model is used, each system vulnerability can be respectively computed, denoted by VT P, , VT G, from the to- pology-based model and VF P, , VF G, from the flow-based model. In addition, for interdependent systems, this paper uses the weighted sum of two system vulnerability to reflect their overall vulner- ability, with equations as follow:

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116 109

= * + ( − )*

= * + ( − )* ( )

V w V w V

V w V w V

1

1 6

T PG T P T G

F PG F P F G

, , ,

, , ,

To estimate overall vulnerability of multiple systems and identify their critical locations, this paper simply sets w¼0.5 for illustrative purpose, its setting in practice can be made by using the Analytic Hierarchy Process (AHP) method [41].

3. Modeling spatially localized attacks and identifying critical locations

To model spatially localized attacks, it needs to first define spatially localized areas and then describe how components within an area would fail. As a starting point, this paper considers circle shaped spatially localized areas, where a localized area has a circle shape determined by a center and a radius r, and assume all components located within or crossing the area are all failed. These assumptions can be relaxed and will be further discussed in the Section 5. Some other types of spatially localized attacks have been studied by the author in another work [28]. This model can simulate bomb and military attacks with mass destruction. Note that the attack center in this model can be in any position of the infrastructure map, and then the number of candidate attack centers is infinite. Because the number of infrastructure compo- nents is limited, as an alternative, the problem of identifying cri- tical location can be converted to first search a set of infrastructure

Fig. 1. A simple network and several spatially localized attack scenarios. The attack ra components n1, e1, e2, e6; (b) an attack scenario which covers a maximal component set set and has its circle boundary pass through two nodes n1 and n2; (d) an attack scenario boundary tangent to an edge e6 and pass through one node n4.

components (attack objectives) that can be covered by a circle with radius r and can cause the largest vulnerability than all other component sets (each of them can be covered by a circle r).

To better introduce the algorithm, this paper first takes a simple example shown in Fig. 1 for demonstration. The network in Fig. 1 has five nodes and seven edges, and the four circle shaped attack scenarios have the same radius. In Fig. 1(a), components n1, e1, e2, e6 is a set of components which can be covered by a circle (r), but this set of components is not the optimum attack objective, because if moving the attack center into some position like in Fig. 1 (b), three more components n2, e3, e4 can be covered by a circle (r). This indicates the attack scenario in Fig. 1(b) can cause larger vulnerability than the scenario in Fig. 1(a), and then the attack location in Fig. 1(b) is more critical than that in Fig. 1(a). The next question is whether the attack center in Fig. 1(b) can further move to some position to cover more components besides n1, n2, e1, e2, e3, e4, e6. The answer is negative, then these seven components {n1,n2,e1,e2,e3,e4,e6} can be called a maximal component set, which is defined as a set of components that can be covered by a circle(r) and any new component added into will make the com- ponents not able to be covered by a circle(r). Given a network position and an attack radius r, if it is able to identify all its maximal component sets MCSs, then among them, the MCS whose contained component failures can cause the largest vulnerability is the optimum attack objective. Hence, the problem is further con- verted into how to identify all maximal component sets.

Note that each maximal component set MCSi can be covered by

dius is the same in four cases. (a) an attack scenario which causes the failures of {n1, n2, e1, e2, e3, e4, e6}; (c) an attack scenario which covers a maximal component which covers a maximal component set {n2,n4,e1,e3,e4,e5,e6,e7} and has its circle

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116110

a circle(r) with its center at any position in an area described by a series of inequalities (the spherical distance between the attack center and each component of the attack objectives is less than or equal to r), and then the circle center can move into some position so that it can pass through or be tangent to at least two compo- nents (two nodes or two edges or one node and one edge) in MCSi. For example, the attack center in Fig. 1(b) can move to the position shown in Fig. 1(c) that can make the seven components still cov- ered and make its boundary pass through two nodes n1 and n2. Fig. 1(d) shows a maximal component set {n2, n4, e1, e3, e4, e5, e6, e7} covered by a circle, with its boundary tangent to the edge e6 and passing through the node n4. Based on this fact, to identify all maximal component sets, it can be realized according to the fol- lowing procedures:

Step One: for any two infrastructure components, find out the attack center with radius r and its circle boundary passing through or tangent to these two components. Depending on the two components’ type, there exist three cases: (a) two components are both nodes (including the case that the two nodes are the same node); (b) one component is a node and the other is an edge; (c) two components are both edges. For each of these cases, this paper next introduces the equations to compute the attack center. The infrastructure component location information is usually mea- sured by latitude and longitude, which actually corresponds to a three-dimensional coordinate system. To make the attack center related equations easier to be solved, the manuscript uses a three- dimensional coordinate system for calculation. For any node with longitude g and latitude t, its location can be converted into a point (x, y, z), with x¼cos(t)*sin(g), y¼cos(t)*cos(g), z¼sin(t). When the two components in step one are both nodes, if these two nodes are the same node, the attack center is just the position of the node; if these two nodes are different nodes and if their spherical distance is not larger than 2r (as the infrastructure components are on the earthquake or ball surface, the shortest path between two components is measured in terms of spherical distance), denote their points in three dimensional space by (x1, y1, z1) and (x2, y2, z2), respectively and the earth radius by RE6400 km, then the attack center (x0, y0, z0) with its circle boundary passing through these two nodes satisfies the following equations:

+ + =

( − ) + ( − ) + ( − ) =

( − ) + ( − ) + ( − ) = ( )

⎜ ⎟

⎜ ⎟

⎧

⎨

⎪ ⎪⎪

⎩

⎪ ⎪ ⎪

⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥

⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥

x y z R

x x y y z z R r R

x x y y z z R r R

2 sin 2

2 sin 2 7

0 2

0 2

0 2 2

0 1 2

0 1 2

0 1 2

2

0 2 2

0 2 2

0 2 2

2

This equation can be converted into the following equations

+ + = + + = + + = ( )

⎧ ⎨ ⎪

⎩ ⎪

x y z R

a x a y a z a

b x b y b z b 8

0 2

0 2

0 2 2

1 0 2 0 3 0 4

1 0 2 0 3 0 4

With the parameters a1, a2, a3, a4, b1, b2, b3, b4 satisfying the following equations:

= = = = −

= = = = − ( )

⎜ ⎟

⎜ ⎟

⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥

⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥

a x a y a z a R R r R

b x b y b z b R R r R

2 , 2 , 2 , 2 2 sin 2

2 , 2 , 2 , 2 2 sin 2 9

1 1 2 1 3 1 4 2

2

1 2 2 2 3 2 4 2

2

By solving the above equations, the circle center can be de- termined as:

= −( + ) ± ( + ) − ( + + )( + − )

+ +

= +

= +

= − −

= − −

= − −

= − +

( )

x c d e f c d e f c e d f R

c e

y c x d

z e x f

c a b a b a b a b

d a b a b a b a b

e a a c

a

f a d a

a

1

1

10

0 1 1 1 1 1 1 1 1

2 1 2

1 2

1 2

1 2 2

1 2

1 2

0 1 0 1

0 1 0 1

1 1 3 3 1

3 2 2 3

1 3 4 4 3

3 2 2 3

1 1 2 1

3

1 2 1 4

3

When the two components in step one are one node C (x3, y3, z3) and one edge connecting two nodes A (x1, y1, z1) and B (x2, y2, z2), then the circle center (x0, y0, z0) with its circle boundary passing through the node C and tangent to the edge AB also sa- tisfies the Eq. (8), with the parameters a1, a2, a3, a4, b1, b2, b3, b4 provided by the following equations:

= − = − = −

= ± ( )

+ + ( + + ) = =

= = − ( )

⎜ ⎟ ⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥

a y z y z a x z x z a x y x y a

R r R a a a

a a a b x b y b

z b R R r R

, , ,

sin / 2 , 2 ,

2 , 2 2 sin 2 11

1 1 2 2 1 2 2 1 1 2 3 1 2 2 1 4

2 2

1 2

2 2

3 2 1 2 3 1 3 2 3 3

3 4 2

2

When the two components in step one are both edges, one connecting nodes A(x1, y1, z1) and B(x2, y2, z2) and the other connecting nodes C(x3, y3, z3) and D(x4, y4, z4), then the circle center (x0, y0, z0) with its circle boundary tangent to both two edges AB and CD also satisfies the Eq. (8), with the parameters a1, a2, a3, a4, b1, b2, b3, b4 provided by the following equations:

= − = − = −

= ± ( )

+ + ( + + ) = −

= − = −

= ± ( )

+ + ( + + )

( )

a y z y z a x z x z a x y x y a

R r R a a a

a a a b y z y z b

x z x z b x y x y b

R r R b b b

b b b

, , ,

sin / ,

, ,

sin /

12

1 1 2 2 1 2 2 1 1 2 3 1 2 2 1 4

2 2

1 2

2 2

3 2 1 2 3 1 3 4 4 3 2

4 3 3 4 3 3 4 4 3 4

2 2

1 2

2 2

3 2 1 2 3

Step Two: for each of the attack centers computed in step one, compute the spherical distance between the attack center and all infrastructure components to identify the components that can be covered by a circle r (with the distance less than or equal to r), and then these components are a candidate maximal component set. If the infrastructure component is a node (xi, yi, zi), whether it can be covered by a circle can be judged by the following inequality:

( − ) + ( − ) + ( − ) ≤ ( )

⎜ ⎟ ⎡ ⎣⎢

⎛ ⎝

⎞ ⎠

⎤ ⎦⎥x x y y z z R

r R

2 sin 2 13

i i i0 2

0 2

0 2

2

If the infrastructure component is an edge connecting two nodes A(x1, y1, z1) and B(x2, y2, z2), whether this edge can intersect with the circle can be checked by the shortest spherical distance between the circle center and the edge, as the edge is a part of the following curve

+ + = ( − ) + ( − ) + ( − ) = ( )

⎪

⎪⎧⎨ ⎩

x y z R

y z y z x x z x z y x y x y z 0 14

2 2 2 2

1 2 2 1 2 1 1 2 1 2 2 1

The point C (x, y, z) on this curve that has the shortest spherical distance to the circle center will satisfy the following equation:

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116 111

[( − ) − ( − ) ]

+ [( − ) − ( − ) ]

+ [( − ) − ( − ) ] = ( )

x z x z z x y x y y x

x y x y x y z y z z y

y z y z y x z x z x z 0 15

2 1 1 2 0 1 2 2 1 0

1 2 2 1 0 1 2 2 1 0

1 2 2 1 0 2 1 1 2 0

Solving the above Eqs. (14) and (15) can get the location of point C, then according to the locations of nodes A, B and C, if the spherical distance between A and B is the sum of spherical dis- tance between A and C and the spherical distance between C and B, then the shortest spherical distance from the circle center to the edge AB is the spherical distance from the circle center to C; if not, it is the minimum value of the spherical distance between A and the attack center and the spherical distance between B and the attack center.

Step three: refine all candidate maximal component sets so that each component set is not the proper subset of any other component set, then the refined candidate maximal component sets are just all maximal component sets. For each maximal component set, its contained components’ failures induced vul- nerability can be computed for both single and multiple systems. The optimum maximal component set that causes the largest

Fig. 2. Critical locations in single systems and interdependent systems in Harris County analysis: (a) power system; (b) gas system; (c) interdependent power and gas systems.

vulnerability can be identified, and the critical location is an area that can cover all components in the optimum maximal compo- nent set. Note that the critical location can be characterized by critical location center (x, y, z) determined by a series of inequal- ities, which describe the shortest spherical distance between (x, y, z) and each component in the optimum maximal component set is less than or equal to r.

4. Applications and results

This section will take interdependent power and gas systems in Harris County, Texas, USA as an example to identify critical loca- tions in these systems and make pertinent vulnerability analysis from four aspects: effects of attack radius on critical locations and their associated vulnerability, comparison of vulnerability under spatially localized attacks and equivalent random failures, com- parison of vulnerability by considering purely node failures and by considering both node and edge failures in the attack area, com- parison of critical locations identified from topology-based and flow-based models.

, Texas, USA when the attack radius is 2.5 km and the flow-based model is used for

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116112

First, given an attack radius r¼2.5 km (this value is selected for illustrative purpose), critical locations in single systems and in- terdependent systems can be identified, as shown in Fig. 2. When only considering single power system, the critical location is at right middle side of the map and it causes the largest flow-based power system vulnerability VF,P¼0.2288. If attacking this critical location and further considering gas system failures as well as interdependencies, it causes a flow-based overall vulnerability VF,PG¼0.1654. When only considering single gas system, the cri- tical location is at right bottom side of the map and it causes the largest flow-based gas system vulnerability VF,G¼0.2653. If at- tacking this critical location and further considering power system failures as well as interdependencies, it causes a flow-based overall vulnerability VF,PG¼0.1509. When considering both two systems, the critical location is at left upper side of the map and it causes the largest flow-based overall vulnerability VF,PG¼0.1736, which is 5% (15%) larger than attacking the critical location in single power (gas) system. These results indicate that infra- structure interdependencies largely affect the position of critical locations.

Second, with the increase of attack radius r, how do the critical location and its associated vulnerability change? To address this problem, this paper identifies the critical locations and records their vulnerability when r ranges from 0.5 km to 25 km with a step of 0.5 km. Fig. 3(a) shows three typical critical locations for the studied interdependent systems when the flow-based vulner- ability model is used for analysis. When rr2.5 km, the critical locations are at the left upper side of the map with attack centers very close to the center for r¼2.5 km; when 2.5 kmorr10 km, the critical locations are at the right side of the map with attack centers around the center for r¼10 km; when r becomes larger and larger, the critical location center moves to the middle of the map, as shown in the figure for r¼25 km. For single systems, si- milar trends are also found. Fig. 3(b) shows the flow-based vul- nerability curves when attacking the critical locations in single and multiple systems. For comparison purpose, for critical locations identified in single systems, their induced overall flow-based vulnerability are also computed and shown in the figure. From the figure, for single systems, the critical locations induced power or gas system vulnerability is non-decreasing function, and at some range of r, the largest vulnerability keeps constant, which indicates the critical locations are still in the same area and the optimum

Fig. 3. (a) Critical locations in interdependent power and gas systems when the attack r locations in single systems and multiple systems are attacked. For critical locations identi shown in the figure.

maximal component sets include almost the same components. If attacking critical locations in single systems and further con- sidering another system failures as well as interdependencies, it is interesting to find that the overall vulnerability curve is not a non- decreasing function. In Fig. 3(b), when r increases from 18.5 km to 19 km, attacking the critical location for single power system will cause a sharply decrease of the overall vulnerability, this is be- cause attacking the power critical location with r¼18.5 km (it causes direct failures of 99 gas components) will directly affect much more gas components than that with r¼19.0 km (it causes direct failures of 20 gas components), which makes gas system vulnerability decrease from 0.8367 to 0.2041, then the overall vulnerability decreases from 0.7418 to 0.4260. It indicates that critical location for one system is not necessary critical for another system, and critical location in interdependent systems depends on how multiple systems are geographically interdependent. When considering both systems, the overall vulnerability curve is also a non-decreasing function with VF,PG increasing from 0.1736 at r¼0.5 km to 0.7988 at r¼25 km.

Third, as studied by Berezin et al. [4], for interdependent spa- tially embedded networks, in the case of only considering system topologies and taking the largest component size as the perfor- mance metric, spatially localized attacks with attack centers ran- domly selected in the network may be more harmful than equivalent random failures. For the studied systems, when their flow properties and flow-based performance model are con- sidered, whether the above results hold true will be studied next. This paper uses the ARCGIS technique to generate 100,000 ran- domly points as the attack centers in the Harris County territory, then the averaged flow-based vulnerability under all these sce- narios is estimated. For each r, the average fraction of node or edges covered in all 100,000 attack area scenarios is also com- puted and taken as equivalent intensity for random node or edge failures. Based on this equivalent random failure intensity f, 10,000 failure scenarios are generated through comparing a uniformly distributed random number with f for each node or edge, and these failure scenarios induced vulnerability are estimated. Fig. 4 shows the averaged flow-based vulnerability results under ran- dom location based localized attacks, and random node or edge failures with equivalent intensities. From the figure, for single systems, the random location based attacks are less detrimental than random node failures, but are more harmful than random

adius is 2.5 km, 10 km and 25 km. (b) Flow-based vulnerability curves when critical fied in single systems, their induced overall vulnerability for multiple systems is also

Fig. 4. Averaged flow-based vulnerability under random location based spatially localized attacks, random node and edge failures with equivalent intensities. (a) power system; (b) gas system; (c) interdependent power and gas system.

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116 113

edge failures when r is less than 22.5 km for power system and 17 km for gas system. For interdependent systems, random loca- tion based attacks are less harmful than both random node and edge failures, which are opposite with the results found by Berezin et al. [4]. To check the reason, this paper further uses the purely topology-based model and the largest component size in Berezin et al. [4] to compare system vulnerability under spatially localized attacks and equivalent random failures, and it gets similar results with those in this paper. Hence, the opposite results could be mainly because the two systems studied in this paper do not be- long to the type of networks that can make spatially localized attacks more harmful.

Fourth, for previous research on spatially localized attacks, some map partition based modeling approaches [19,30] consider the failures of some special links crossing the area but with their endpoints outside the attack areas, but some theoretical based approaches only consider node failures in the attack areas, the failures of those special links are not considered [4]. If considering the failures of those special links, whether the critical locations and their associated vulnerability are affected will be discussed

next. Given an attack radius, for each of its associated maximal component set, a candidate maximal node set that only considers node failures can be obtained from the maximal component set by deleting edge components; then for all derived candidate maximal node sets, refining them to make each of them non-void proper subset of any other will produce all maximal node sets. Fig. 5(a) shows the critical locations for only considering node failures and for considering both node and edge failures in the attack areas with r¼3.5 km when the flow-based vulnerability model is used for analysis. Despite the two critical locations are close to each other, they actually include many different components and cause different vulnerability. The critical location for considering only node failures causes an overall flow-based vulnerability 0.2121, while the critical location for considering both node and edge failures causes an overall flow-based vulnerability 0.2445. Similar results can be also found for single systems under some values of r, but at most of r, the critical locations and their associated vul- nerability do not change, as shown in Fig. 5(b) that depicts the differences of the largest flow-based vulnerability under the two cases. In most of r, the critical location does not change, but in a

Fig. 5. (a) Critical locations in interdependent power and gas systems when considering only node failures and considering both node and edge failures in the attack area; (b) the largest vulnerability difference between considering only node failures and considering both node and edge failures.

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116114

few of r, the vulnerability difference is not zero, which means the critical locations in the two cases are in different positions.

Fifth, based on topology-based and flow-based vulnerability models, it enables to compare the critical locations identified from these two models and analyze whether the topology-based critical locations are still critical from the flow-based perspective. Fig. 6(a) shows critical locations in interdependent systems when different vulnerability models are used at r¼2.5 km. It can be found that these two locations are far from each other. The topology-based critical location causes the largest topology-based overall vulner- ability VT,PG¼0.4219 but a flow-based vulnerability VF,PG¼0.1407, while the flow-based critical location causes the largest flow- based overall vulnerability VF,PG¼0.1736 but a topology-based vulnerability VT,PG¼0.2238. Fig. 6(b) shows the flow-based vul- nerability curves by attacking both topology-based and flow-based critical locations. It can be found that under many attack radius r, attacking the topology-based critical location can also cause the flow-based vulnerability close to or equal to the largest flow-based vulnerability, which indicates topology-based critical locations are also critical from the flow-based perspective in these cases.

Fig. 6. (a) Critical locations in interdependent systems when topology-based and flow-ba based vulnerability curve when attacking topology-based and flow-based critical locatio

5. Discussions

The results in Section 4 show that (1) Infrastructure inter- dependencies and attack radius largely affect the position of cri- tical locations; (2) Spatially localized attacks on interdependent systems are less harmful than equivalent random failures; (3) In most values of attack radius, critical locations by considering only node failures do not change when considering both node and edge failures in the attack area; (d) for many values of attack radius critical locations identified by topology-based model are also cri- tical from the flow-based perspective. These results are based on several assumptions, and this section will further discuss these assumptions and consider possible extensions.

(1) For the infrastructure system vulnerability model, this paper only considers the fraction of power supply and the fraction of gas supply as the performance metrics to quantify system flow-based vulnerability. However, infrastructure systems usually have many stakeholders with different preferences, and one performance metric alone is not sufficient to reflect

sed performance models are used respectively at attack radius r¼2.5 km; (b) Flow- ns.

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116 115

their overall expectation. Hence, considering other types of performance metrics and combining them into one composite vulnerability metric as the attack objective function can help identify critical locations that can better reflect the concerns from all stakeholders. Note that in this case, searching all maximal component sets by the proposed approach is still the first step for identifying critical locations.

(2) This paper assumes that the attack areas are circle shaped and all components located within or crossing the area are all failed. The first assumption can be relaxed and the proposed algorithm can be extended to other shapes of attack areas, such as triangle, square, hexagon. The solution is similar and can be realized by first searching all maximal component sets and then identifying the optimum maximal component set that can cause the largest vulnerability, but in these cases the equations required to search the maximal component set ac- cording to two component locations are more complicated than Eqs. (7)–(12). The second assumption can be also relaxed. If each component in the attack area has the same failure probability, identifying critical location can still start from searching the maximal component sets. If components have different failure probabilities depending on their locations to the attack center, the proposed method is no longer valid and then the infrastructure map partition method can be used as an alternative [19,30]. Also, this paper does not consider ground variation, such as hills, valleys, towers, which may affect the failure probabilities of components in the attack area. When data related to these factors is available, combing them into analysis will be more realistic, but note that the proposed method can at least provide an upper-bound esti- mation of the vulnerability.

(3) In the case study, this paper considers an extensive range of attack radius from 0.5 km to 25 km for analysis. Large attack radius may be not realistic as they could be only caused by nuclear weapons. But note that the systems studied in this paper are both transmission systems and they cover a large- scale area, to make these systems significantly affected by spatially localized attacks, it needs to assume a relatively large attack radius. Hence, the extensive range of attack radius is considered mainly for illustrative purpose. If the data for a distribution system in a relatively small area is available for analysis, like infrastructure systems in a university campus [30], using the proposed approach under small attack radius for analysis will be more realistic.

(4) This paper considers not only node failures, but also edge failures within the attack area. However, the proposed algo- rithm is based on the assumption that each edge between two nodes is deployed along the shortest path in terms of spherical distance. This assumption makes sense in most cases due to the constraint of construction cost, but in some case, the edge may be not deployed according to the shortest path due to geologic, political or some other reasons. In this case, it can divide the edge into a series of short segments, each of which are shortest path between its endpoint nodes, and then by using the proposed approach in Section 2, the critical locations can be also identified.

6. Conclusions

This paper proposes an approach to identify critical locations in interdependent infrastructure systems, and the critical location is characterized by critical location center determined by a series of inequalities, which describe the shortest spherical distance be- tween the center and each component in the optimum maximal component set is less than or equal to the attack radius. To better

protect infrastructure systems, increasing the security inspection in the critical location center area or hardening components in the optimum attack objectives are possible mitigation strategies. De- spite this paper only takes the interdependent power and gas systems in Harris County, Texas, USA for analysis, the proposed approach can be easily adapted to analyze other types of inter- dependent systems. Also, this paper only considers uni-directional physical interdependencies between multiple systems for analysis, the proposed approach can be also used for the case of bi-direc- tional physical interdependencies. Moreover, the proposed algo- rithm can be also used to model the interactions between at- tackers and defenders [21] for identifying the optimal mitigation strategies.

Acknowledgments

This material is based upon work supported in part by the National Natural Science Foundation of China under Grant 51208223, 61572212 and the Fundamental Research Funds for the Central Universities under Grant 2014QN166. Any opinions, find- ings, and conclusions or recommendations expressed in this ma- terial are those of the author and do not necessarily reflect the views of the sponsors. The author also acknowledges the data shared by Dr. Dueñas-Osorio through his Structural and Infra- structure Reliability and Risk Assessment (SISRRA) research group at Rice University.

References

[1] Adachi T, Ellingwood BR. Serviceability of earthquake-damaged water sys- tems: effects of electrical power availability and power backup systems on system vulnerability. Reliab Eng Syst Saf 2008;93:78–88.

[2] Alderson DL, Brown GG, Carlyle WM. Operational models of infrastructure resilience. Risk Anal 2015;35(4):562–86.

[3] Baxter GJ, Dorogovtsev SN, Goltsev AV, Mendes JFF. Avalanche collapse of interdependent networks. Phys Rev Lett 2012;109:248701.

[4] Berezin Y, Bashan A, Danziger MM, Li D, Havlin S. Localized attacks on spatially embedded networks with dependencies. Sci Rep 2015;5:8934.

[5] Buldyrev SV, Parshani R, Paul G, Stanley HE, Havlin S. Catastrophic cascade of failures in interdependent networks. Nature 2010;464:1025–6.

[6] Buldyrev SV, Shere N, Cwilich GA. Interdependent networks with identical degrees of mutually dependent nodes. Phys Rev E 2011;83:016112.

[7] Carvalho R, Buzna L, Bono F, Gutierrez E, Just W, Arrowsmith D. Robustness of trans-European gas networks. Phys Rev E 2009;80:016106.

[8] Cavallaro M, Asprone D, Latora V, Manfredi G, Nicosia V. Assessment of urban ecosystem resilience through hybrid social–physical complex networks. Comput-Aided Civ Infrastruct Eng 2014;29(8):608–25.

[9] Cavalieri F, Franchin P. Models for seismic vulnerability analysis of power networks: comparative assessment. Comput-Aided Civ Infrastruct Eng 2014;29:590–607.

[10] Cavdaroglu B, Mitchell JE, Sharkey TC, Wallace WA. Integrating restoration and scheduling decisions for disrupted interdependent infrastructure systems. Ann Oper Res 2013;203:279–94.

[11] Dueñas-Osorio L, James IC, Barry JG. Seismic response of critical inter- dependent networks. Earthq Eng Struct Dyn 2007;36:285–306.

[12] Dueñas-Osorio L, James IC, Barry JG, Ann B. Interdependent response of net- worked systems. J Infrastruct Syst 2007;13(3):185–94.

[13] Federal Emergency Management Agency (FEMA), Hazards U.S. Multi-Hazard (HAZUS- MH) Assessment Tool v1.4, ⟨www.fema.gov/plan/prevent/hazus/in dex.shtm⟩; 2015.

[14] Franchin P, Cavalieri F. Probabilistic assessment of civil infrastructure resi- lience to earthquakes. Comput-Aided Civ Infrastruct Eng 2015;30:583–600.

[15] Gao J, Buldyrev SV, Havlin S, Stanley HE. Robustness of a network of networks. Phys Rev Lett 2011;107(19):195701.

[16] Hernandez I, Dueñas-Osorio L. Sequential propagation of seismic fragility across interdependent lifeline systems. Earthq Spectra 2011;27(1):23–43.

[17] Hines P, Cotilla-Sanchez E, Blumsack S. Do topological models provide good information about electricity infrastructure vulnerability? Chaos 2010;20:033122.

[18] Huang X, Gao J, Buldyrev SV, Havlin S, Stanley HE. Robustness of inter- dependent networks under targeted attack. Phys Rev E 2011;83:065101(R).

[19] Johansson J, Hassel H. An approach for modeling interdependent infra- structures in the context of vulnerability analysis. Reliab Eng Syst Saf 2010;95 (12):1335–44.

M. Ouyang / Reliability Engineering and System Safety 154 (2016) 106–116116

[20] Jönsson H, Johansson J, Johansson H. Identifying critical components in technical infrastructure networks. J Risk Reliab 2008;222(2):235–43.

[21] Murray-Tuite PM. A methodology for assessing transportation network ter- rorism risk with attacker and defender interactions. Comput-Aided Civ In- frastruct Eng 2010;25:396–410.

[22] Ouyang M, Dueñas-Osorio L. An approach to design interface topologies across interdependent urban infrastructure systems. Reliab Eng Syst Saf 2011;96 (11):1462–73.

[23] Ouyang M. Comparisons of purely topological model, betweenness based model and direct current power flow model to analyze power grid vulner- ability. Chaos 2013;23:023114.

[24] Ouyang M. Review on modeling and simulation of interdependent critical infrastructure systems. Reliab Eng Syst Saf 2014;121:43–60.

[25] Ouyang M, Zhao LJ, Hong L, Pan ZZ. Comparisons of complex network based models and real train flow model to analyze Chinese railway vulnerability. Reliab Eng Syst Saf 2014;123:38–46.

[26] Ouyang M, Wang ZH. Resilience assessment of interdependent infrastructure systems: With a focus on joint restoration modeling and analysis. Reliab Eng Syst Saf Spec Issue Resil Eng 2015;141:74–82.

[27] Ouyang M, Pan ZZ, Hong L, He Y. Vulnerability analysis of complementary transportation systems with applications to railway and airline systems in China. Reliab Eng Syst Saf 2015;142:248–57.

[28] Ouyang M, Tian H, Hong L, Wang ZH, Mao ZJ. Critical Infrastructure Vulner- ability to Spatially Localized Failures with Applications to Chinese Railway System. Risk Analysis 2016 (under review).

[29] Parshani R, Buldyrev SV, Havlin S. Interdependent networks: reducing the coupling strength leads to a change from a first to second order percolation transition. Phys Rev Lett 2010;105:048701.

[30] Patterson SA, Apostolakis GE. Identification of critical locations across multiple infrastructures for terrorist actions. Reliab Eng Syst Saf 2007;92:1183–203.

[31] Peerenboom JP, Fisher RE, Rinaldi SM, Kelly TK. Studying the chain reaction. Electr Perspect 2002;27(1):22–35.

[32] Poljansek K, Bono F, Gutierrez E. Seismic risk assessment of interdependent critical infrastructure systems: The case of European gas and electricity

networks. Earthq Eng Struct Dyn 2012;41:61–79. [33] Rinalidi SM, Peerenboom JP, Kelly T. Identifying, understanding and analyzing

critical infrastructure interdependencies. IEEE Control Syst Mag 2001;12:11– 25.

[34] Rosato V, Issacharoff L, Tiriticco F, Meloni S. Modeling interdependent infra- structures using interacting dynamical models. Int J Crit Infrastruct ((IJCI)) 2008;4(1-2):63–79.

[35] Schrijver A. On the History of the Transportation and Maximum Flow Pro- blems. Math Program Ser B 2002;91(3):437–45.

[36] Shao J, Buldyrev SV, Havlin S, Stanley HE. Cascade of failures in coupled net- work systems with multiple support-dependent relations. Phys Rev E 2011;83 (3):036116.

[37] Svendsen NK, Wolthusen SD. Connectivity models of interdependency in mixed-type critical infrastructure networks. Inf Sec Tech Rep 2007;12(1):44– 55.

[38] Zio E, Sansavini G. Modeling interdependent network systems for identifying cascade-safe operating margins. IEEE Trans Reliab 2011;60(1):94–101.

[39] Zio E, Golea LR, Rocco CM. Identifying groups of critical edges in a realistic electrical network by multi-objective genetic algorithms. Reliab Eng Syst Saf 2012;99:172–7.

[40] Zio E, Golea LR, Sansavini G. Optimizing protections against cascades in net- work systems: a modified binary differential evolution algorithm. Reliab Eng Syst Saf 2012;103:72–83.

[41] Saaty TL. The analytic hierarchy process: planning, priority setting, resource allocation.New York: McGraw-Hill; 1980.

[42] Nan C, Eusgeld I, Kröger W. Analyzing vulnerabilities between SCADA system and SUC due to interdependencies. Reliab Eng Syst Saf 2013;113:76–93.

[43] Pinar A, Meza J, Donde V, Lesieutre B. Optimization strategies for the vul- nerability analysis of the electric power grid. SIAM J Optim 2010;20(4):1786– 810.

[44] Hong L, Ouyang M, Peeta S, He XZ, Yan YZ. Vulnerability assessment and mitigation for the Chinese railway system under floods. Reliab Eng Syst Saf 2015;137:58–68.