E Ferrario

N Pedroni

E Zio

Evaluation of the robustness of critical infrastructures by Hierarchical Graph representation, clustering and Monte Carlo simulation

Keywords: Interdependent critical infrastructures, robustness, Hierarchical Graph, hierarchical clustering, multi-state model, system of systems, Monte Carlo simulation, Markov and semi-Markov processes

In this paper, we present a methodological work that adopts a system-of-systems (SoS) viewpoint for the evaluation of the robustness of interdependent critical infrastructures (CIs).

We propose a Hierarchical Graph representation, where the product flow is dispatched to the demand nodes in consideration of different priorities. We use a multi-state model to describe different degrees of degradation of the individual components, where the transitions between the different states of degradation occur stochastically. The quantitative evaluation of the CIs robustness is performed by Monte Carlo simulation. The methodological approach proposed is illustrated by way of two case studies: the first one concerns small-sized gas and electricity networks and a supervisory control and data acquisition (SCADA) system; the second one considers a moderately large power distribution network, adapted from the IEEE 123 node test feeders. The large size of the second case study requires hierarchical clustering for performing the robustness analysis.

Introduction

Critical infrastructures (CIs) are complex systems essential for providing goods (e.g., energy) and services (e.g., transportation) across local, regional and national boundaries [START_REF] Kröger | Vulnerable Systems[END_REF]. Typically, they present both structural and dynamic complexities. The former derive from i) heterogeneity of components across different technological domains, due to the integration among different systems and ii) scale and dimensionality of connectivity through a large number of components (nodes), highly interconnected by dependences (unidirectional relationships) and interdependences (bidirectional relationships). Dynamic complexity manifests through the emergence of (unexpected) system behavior in response to changes in the environmental and operational conditions of its components. Furthermore uncertainties exist in the failure behavior of CI components, interconnections and interactions, so that the prediction of CI failure behavior is difficult [START_REF] Zio | Uncertainties in smart grids behavior and modeling: What are the risks and vulnerabilities? How to analyze them?[END_REF].

Engineered, physically networked CIs are considered in this paper. Examples are those providing: energy (electricity, oil and gas supply as subsectors); transportation (by rail, road, air, shipping); information and telecommunication (such as the internet); drinking water, including wastewater treatment, etc.

Due to the increasing complexity of CIs, random failures, natural events and malevolent attacks can have severe consequences on health, safety, security, economics and social well-being. In this respect, evaluating the robustness of CIs is fundamental to be able to improve their design and management so to reduce the impacts of disruptive events. There is no unique definition of robustness. Jensen defines it as the degree to which a system can function correctly in the presence of inputs different from those assumed [START_REF] Jensen | Improving Robustness and Flexibility of Tardiness and Total Flowtime Job Shops Using Robustness Measures[END_REF]; for [START_REF] Carlson | Complexity and robustness[END_REF] and [START_REF] Jen | Stable or robust? What's the difference?[END_REF], the robustness guarantees the maintenance of certain desired system characteristics, despite fluctuations in the behavior of its components or in its environment. [START_REF] Jen | Stable or robust? What's the difference?[END_REF] and [START_REF] Ali | Definition of a Robustness Metric for Resource Allocation[END_REF] specify that the concept of robustness should be defined for a given set of system features, under a given set of perturbations applied to the system. According to Foster, robustness is the ability of a system to react to noisy input parameters with little performance degradation [START_REF] Foster | Estimating Complex System Robustness from Dual System Architectures[END_REF]. A recent definition of robustness is given in the glossary of the specialty group on "Foundations of Risk Analysis" of the Society for Risk Analysis, as the antonym of vulnerability [START_REF] Sra | SRA Develops Glossary of Risk-Related Terms[END_REF]. In addition, a system is considered robust to uncertainty if specified goals are achieved, despite large information gaps (information gap is the disparity between what is known, and what needs to be known to ensure specified goals) [START_REF] Sra | SRA Develops Glossary of Risk-Related Terms[END_REF]. In this work, robustness is seen as the capability of the CIs to resist to failures or partial failures of the CIs components assuring the required level (or a high level) of supply of goods or services.

Traditionally, three steps are performed in the robustness assessment of CIs: (i) the system is represented to define the structural, logical and functional relations among its components; (ii) a mathematical model of the system is built to quantify its performance indicators; (iii) the model is solved, e.g., by simulating its behavior under different operational and accidental conditions. With respect to system representation (i), several types of approaches exist in literature, many of which rely on a hierarchy or graph structure. Hierarchical modeling has been often adopted to represent and model complex systems, since many organizational and technology-based systems are hierarchical in nature [START_REF] Courtois | On Time and Space Decomposition of Complex Structures[END_REF][START_REF] Haimes | Modeling complex systems of systems with Phantom System Models[END_REF]. Hierarchical functional models include Goal Tree Success Tree (GTST) [START_REF] Kim | Application of Goal Tree-Success Tree Model as the Knowledge-Base of Operator Advisory Systems[END_REF]]also combined with Master Logic Diagram (MLD) [START_REF] Hu | Logic-based hierarchies for modeling behavior of complex dynamic systems with applications[END_REF] and Multilevel Flow Models (MFM) [Lind, 2011a;[START_REF] Lind | Reasoning about causes and consequences in Multilevel Flow Models[END_REF]. In risk analysis, common representation techniques are hierarchical trees that are possibly used to identify i) the initiating causes of a pre-specified, undesired event or ii) the accident sequences that can generate from a single initiating event through the development of structured logic trees (i.e., fault and event trees, respectively) [START_REF] Zio | An Introduction to the Basics of Reliability and Risk Analysis[END_REF]. Recently, also networks have been represented by hierarchical modeling [START_REF] Gómez | Hierarchical infrastructure network representation methods for risk-based decision-making[END_REF][START_REF] Buriticá Cortés | A hierarchy-based approach to seismic vulnerability assessment of bulk power systems[END_REF].

In complex network theory approaches, instead, complex systems are represented by networks where the nodes stand for the components and the links describe the physical and relational connections among them [START_REF] Gheorghe | Critical Infrastructures at Risk. Securing the European Electric Power System[END_REF][START_REF] Dueñas-Osorio | Cascading failures in complex infrastructure systems[END_REF][START_REF] Hernandez-Fajardo | Probabilistic study of cascading failures in complex interdependent lifeline systems[END_REF][START_REF] Fang | Comparing Network-Centric and Power Flow Models for the Optimal Allocation of Link Capacities in a Cascade-Resilient Power Transmission Network[END_REF][START_REF] Praks | Monte-Carlo Based Reliability Modelling of a Gas Network Using Graph Theory Approach[END_REF]. Network-based approaches model interdependent CIs on the basis of their topologies or flow patterns [START_REF] Ouyang | Review on modeling and simulation of interdependent critical infrastructure systems[END_REF]. Physical and non-physical (heterogeneous) overlapping infrastructures have been represented as networks to identify risk scenarios and the impacts on connected networks in [START_REF] Lambert | Terrorism Scenario Identification by Superposition of Infrastructure Networks[END_REF]. Also probabilistic methods (e.g., Petri nets [START_REF] Laprie | Modelling interdependencies between the electricity and information infrastructures[END_REF], Bayesian networks [START_REF] Di Giorgio | A Bayesian Network-Based Approach to the Critical Infrastructure Interdependencies Analysis[END_REF] and flowgraphs [START_REF] Huzurbazar | Flowgraph models: a Bayesian case study in construction engineering[END_REF]) are based on graph representations.

In this paper, we present a methodological work that embraces a system-of-systems (SoS) framework of analysis [START_REF] Gheorghe | Mining intelligence data in the benefit of critical infrastructures security: vulnerability modelling, simulation and assessment, system of systems engineering[END_REF][START_REF] Eusgeld | System-of-systems" approach for interdependent critical infrastructures[END_REF][START_REF] Kröger | Vulnerable Systems[END_REF][START_REF] Tsilipanos | A System of Systems Framework for the Reliability Assessment of Telecommunications Networks[END_REF] and propose a Hierarchical Graph representation to evaluate the robustness of interdependent CIs, measured by its capability to deliver the required amount of product (e.g., energy, water, etc.) to the demand nodes of the infrastructure. In this respect, the system can be considered robust if it can maintain the required level (or a high level) of delivery when it is affected by perturbations (failures and partial failures). In doing so, we take into account the fact that the demand nodes may have different importances, which leads to possibly different priorities in the distribution of the product flow through the connections to the elements of the CI. For example, hospitals may be considered more important than residential buildings given their role in the health-care system; as a consequence, in the case of a reduction of electric power that can be delivered in the network, hospitals may receive the priority with respect to houses. This ranking of priority should be fixed by the analyst and, then, criteria (hereafter also referred to as "importance criteria") for the partition of the flow (e.g., electric power) in the network can be defined. In this work, we assume three different importance criteria that depend on the geographic position of the demand nodes, the quantity of product required by each of them and the assumption of equality of the demand nodes.

The representation proposed consists of a graph structured in hierarchical levels that allows highlighting critical arcs and supporting the quantitative robustness evaluation by assigning different priorities to the demand nodes. Critical arcs are here defined as those links whose interruption or degradation affects several demand nodes. This concept of criticality can be related to that of "importance measure" used in reliability theory. Actually, importance measures quantify the contribution of a given component to a properly selected measure of system performance (e.g., robustness in this case): see, for example, the Birnbaum [START_REF] Birnbaum | On the importance of different components in a multicomponent system[END_REF], Fussell and Vesely [START_REF] Fussell | How to calculate system reliability and safety characteristics[END_REF], criticality importance measures [START_REF] Elsayed | Reliability Engineering[END_REF], etc.

More specifically, with respect to network system analysis, other importance measures have been defined to measure component criticality [START_REF] Borgatti | Centrality and network flow[END_REF][START_REF] Kröger | Vulnerable Systems[END_REF], like classical topological centrality measures including the degree of centrality [START_REF] Natvig | Multistate Systems Reliability Theory with Applications[END_REF][START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF], the closeness centrality [START_REF] Sabidussi | The centrality index of a graph[END_REF][START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF][START_REF] Wasserman | Social Network Analysis: Methods and Applications[END_REF], the betweeness centrality [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF] and the information centrality [START_REF] Latora | A measure of centrality based on network efficiency[END_REF].

For a more realistic representation, we adopt a multi-state model where different degrees of degradation of the individual components are contemplated [START_REF] Ferrario | Goal Tree Success Tree-Dynamic Master Logic Diagram and Monte Carlo simulation for the safety and resilience assessment of a multistate system of systems[END_REF]Zio, 2014] [Natvig, 2010]; the transitions between the different states of degradation occur stochastically and are modeled within Markov and semi-Markov processes.

For illustration purposes, we consider two case studies: the first one is characterized by smallsized interconnected gas and electricity networks and a supervisory control and data acquisition (SCADA) system [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF]; the second one is adapted from the IEEE 123 node test feeders [IEEE, 2000] and includes a large electricity distribution network. The first case study is chosen small enough to be able to clearly illustrate the Hierarchical Graph modeling of (three) connected systems, considering different priorities of the demand nodes. The second case study serves the purpose of showing how the approach can be extended when the size of the system increases.

As a measure of the robustness of the system, we evaluate the steady-state probability distributions of the product (e.g., gas and/or electricity) delivered to the demand nodes.

The quantitative evaluation of the system robustness is performed by Monte Carlo (MC) simulation [START_REF] Kalos | Monte Carlo Methods Volume 1: Basics[END_REF][START_REF] Zio | The Monte Carlo Simulation Method for System Reliability and Risk Analysis[END_REF]; in the second case study of larger dimension, an unsupervised spectral clustering algorithm is employed to make the size of the CI manageable and reduce the computational burden related to the analysis [START_REF] Fang | Unsupervised spectral clustering for hierarchical modelling and criticality analysis of complex networks[END_REF].

The remainder of the paper is organized as follows. In Section 2, the Hierarchical Graph representation is introduced and the importance criteria considered are illustrated; in Section 3, the procedural steps to evaluate the robustness of interconnected CIs by Hierarchical Graph and MC simulation are provided, and the combination of Hierarchical Graph and clustering analysis is given, then the advantages and limitations of the approach are discussed; Section 4 contains the description of the two case studies, the representation of the corresponding systems and the results obtained; in Section 5, some conclusions are provided. Finally, in the Appendix the data related to the second case study are illustrated.

Hierarchical Graph representation of systems of systems

The proposed representation technique can be applied to engineered, physically networked CIs (energy, transportation, information and telecommunication) characterized by a radial structure, i.e., by unidirectional flows of "products" (power, water, gas, data). Actually, the representation requires that the CI of interest be first modeled by a directed graph of nodes and arcs without loops (in this case, the arcs may represent elements of an infrastructure or the connections between different infrastructures). Typical radial systems are the distribution networks.

To build the Hierarchical Graph representation, we then, need to distinguish between input, demand (load) and transmission arcs: the "input arcs" connect the production sources to the network, the "demand arcs" terminate with nodes that require a given amount of product, whereas the "transmission arcs" transfer the product to other components in the network. Notice that the transmission and the demand arcs may coincide: for example, an arc may be needed to supply the connected node and in addition it may be required to transmit the product to other arcs/nodes.

In the Hierarchical Graph representation, the adjective "hierarchical" does not imply a "decomposition of the system into different levels of details", as in other hierarchical models (e.g., Goal Tree Success Tree -Dynamic Master Logic Diagram [START_REF] Hu | Logic-based hierarchies for modeling behavior of complex dynamic systems with applications[END_REF]] and hierarchical clustering [START_REF] Gómez | Hierarchical infrastructure network representation methods for risk-based decision-making[END_REF]), but it simply means that the graph of interconnected CIs is structured in hierarchical levels. In extreme synthesis, the representation is built as follows: at the bottom of the graph, the inputs (i.e., the arcs through which the product is injected into the networks) are represented; at the top, the goals (i.e., the demand nodes that have to be satisfied) are reported; in the middle, all the other arcs (transmission and/or load arcs) that provide product to the demand nodes are organized in hierarchical levels. These levels are numbered on the basis of the number of demand nodes that are served by the corresponding arcs: the higher the number of demands supplied by an arc, the higher the hierarchical level of that arc. For example, all the arcs that are required to supply LV demand nodes are "placed" at hierarchical level LV.

Formally, let us consider A interconnected infrastructure systems S (a) , a = 1, …, A, constituting the overall SoS, numbered in order in such a way that the first q exchange physical product (e.g., energy or material) and the last (Aq) exchange information, and are useful for the operation and control of the connected systems (e.g., a supervisory control and data acquisition -SCADA -system). The total number of components (arcs) transmitting physical flow is denoted as N.

For illustration purposes, refer to Figure 1 in which the graph of a SoS (top) and its corresponding Hierarchical Graph (bottom) are reported. The SoS in the example is composed by A = 4 systems, where the first two, i.e., S (1) and S (2) , exchange physical product (solid links in Figure 1, top) and the last two, i.e., S (3) and S (4) , support system S (1) (dotted links in Figure 1, top). The total number of components (arcs) is N = 8.

As described above, the Hierarchical Graph depicts the inputs at the bottom of the representation, i.e., in this case, arc S1 (1) _S2 (1) in Figure 1 (bottom); also, it shows the goals (i.e., the demand nodes) at the top: in this case, the demand nodes are represented by all the nodes of systems S (1) and S (2) , except S1 (1) , which is the source of product. Finally, it organizes the arcs in different hierarchical levels according to the number of demand nodes they supply: for example, in this case arc S1 (2) _S2 (2) is at hierarchical level 4 since it provides product to four demand nodes, i.e., S2 (2) , S3 (2) , S4 (2) and S5 (2) . The quantity of product required by the demand nodes is referred to as Ddem, where the subscript 'dem' is the indicator of the demand node among the N components.

Notice that the arcs referred to the (Aq) control and information systems (which are not contributing to the flow of product, but influence the state of the other arcs) do not appear in the hierarchical structure: instead, they are reported in a trapezoidal frame under the corresponding arc that they affect.

The squares located between the hierarchical levels mean that the product at that level has to be partitioned among the corresponding demand nodes. This representation allows highlighting all the paths going from the input sources to the end nodes: for example, in Figure 2, the path from input S1 (1) to node S3 (2) is highlighted. In addition, the representation is able to put in evidence the critical arcs as those located at higher hierarchical levels, since their interruption or degradation affects more demand nodes: for example, in Figure 2 arc S2 (1) _ S1 (2) is more critical than arc S2 (1) _ S3 (1) since the first one is required to supply five demand nodes (i.e., S1 (2) , S2 (2) , S3 (2) , S4 (2) and S5 (2) ), whereas the second one is necessary just for node S3 (1) . As illustrated in Section 1, in this view the concept of critical arcs here adopted can be related to that of "importance measure" employed in classical reliability theory. Figure 2: Hierarchical Graph of the system of systems in Figure 1, highlighting the path from the input to demand node S3 (2) ; LV: Level.

This representation is useful to analyze the robustness of interdependent CIs taking into account different priorities in the distribution of the product to the demand nodes according to their importance (relevance). In this context, Martinez et al. use path diagrams in representation of the relative importance of network elements across infrastructure systems accounting for the nonquantifiable (qualitative) features of the allocation of resources among the network elements [START_REF] Martinez | Diagramming qualitative goals for multiobjective project selection in large-scale systems[END_REF]. In this work, we neglect the qualitative aspects and we focus only on the quantitative objectives. Notice that the word "importance" is not here related to the concept of "importance measure" used in reliability theory, but it is employed to indicate the relevance of the demand nodes. For example, hospitals may be considered more important than residential buildings, given their role in the health-care system; as a consequence, in the case of a malfunctioning in the electrical transmission line, hospitals may be a priority with respect to common residential areas. Regarding the water distribution systems, a fire-fighting consumption should have the priority with respect to daily water consumption. These priorities of the demand nodes are defined by the analyst and drive the definition of importance criteria needed for the optimal flow partition in the network. However, these importance criteria are not explicitly shown in the representation, which instead is more focused in highlighting the hard, physical constraints that affect the product distribution. Three different importance criteria are considered in this work, namely, sequential, proportional and equal; such criteria are explained hereafter with respect to a simple example consisting of an input of 50 units and two demand nodes S1 and S2 that require 40 and 100 units, respectively (i.e., DS1 = 40 and DS2 = 100).

Sequential importance consists in ranking the demand nodes sequentially on the basis of a chosen "ranking criterion" (e.g., according to their distance from the source node: the closer the node to the source, the higher the importance). In this case, the nodes classified as more important are given higher priority; with respect to the example above, if S1 is more important than S2, the input product is given first to S1 until it is completely supplied and the rest (i.e., 50 -40 = 10 units) is fed to S2. Vice versa if S2 is more important than S1, the input is given entirely to S2 and there is no product left to supply S1.

Proportional importance orders the demand nodes on the basis of the quantity of product they need: the higher their demand, the higher their importance. Then, the product is partitioned into the network according to ratios of importance associated to each demand node, computed as the ratio between the quantity of product required by a node and the sum of all the demands of the entire system. With respect to the example above, S2 is more important than S1 because it requires a higher quantity of product. The ratios of importance are DS2/(DS1 + DS2) = 0.7 for S2

and DS1/(DS1 + DS2) = 0.3 for S1. Then, the input is partitioned as follows: 30% (i.e., 15 product units) is given to S1 and 70% (i.e., 35 product units) to S2.

Equal importance considers the demand nodes equal, even if their demands are different. Thus, the input is partitioned into equal parts on the basis of the number of the demand nodes. With respect to the example above, the input is divided in two equal parts: 50% is given to S1 and 50% to S2.

The detailed operative steps of the algorithm for partitioning the product according to these criteria are given in Section 3.1.

Monte Carlo simulation and Hierarchical Graphs for critical infrastructures evaluation

Within a multi-state system-of-systems (SoS) analysis framework, we wish to evaluate the performance of critical infrastructures (CIs) in terms of robustness, measured by the steadystate probability distribution of the product delivered at the demand nodes of the system. The quantitative evaluation is carried out by combining the Hierarchical Graph representation of Section 2 with Monte Carlo (MC) simulation.

In Section 3.1, the operative steps of the basic procedure are presented. Then, a modification of the basic procedure is proposed in Section 3.2 to deal with CIs of large size: in particular, a clustering algorithm is adopted to pre-process the CI in order to make its size manageable and reduce the computational burden associated to the analysis; in Section 3.3 the advantages and limitations of the approach are discussed.

Operative steps combining Monte Carlo simulation and Hierarchical Graph

We generically denote the state of a component of the CI (i.e., the capacity of the arcs) as ζcomp,i

, i  {1, 2, …, NScomp}, comp = 1, …, N, where N is the total number of components in the SoS, the subscript 'comp' indicates the component of interest, identified by its name or by an integer number from 1 to N, NScomp is the total number of states for component comp, and i is the state identification number (when i = 1, the component is in the worst state, whereas when i = NScomp, it is in the best state). For example, supposing that component S1 ( 2) _S2 ( 2) can enter three possible states, namely 0, 10 and 20, we denote the total number of states for the component as NSS1(2)_S2(2) = 3, and the corresponding states as ζS1(2)_S2( 2),1 = 0, ζS1(2)_S2(2),2 = 10, and ζS1(2)_S2( 2),3 = 20.

The quantity of product requested by the demand nodes is indicated by the vector {Ddem}, dem  {1, …, N}, where the subscript 'dem' indentifies the demand nodes.

In what follows, we describe an algorithm combining Monte Carlo simulation and Hierarchical

Graph representation for the evaluation of the robustness of CIs within a multi-state SoS framework; as mentioned before, in this work the robustness can be seen as the capability of the CIs to resist to failures or partial failures of the arcs assuring the required level (or a high level) of supply of product to the demand nodes and it is quantified in terms of the steady-state probability distribution of the product delivered to the demand nodes.

In extreme synthesis, the algorithm requires as inputs:

 the Hierarchical Graph that allows representing the origin-destination paths and the corresponding arcs in the hierarchical levels (see Section 2);

 the steady state probabilities of transition between the different arc states (i.e., capacities) ζcomp,i , i = {1, 2 ,…, NScomp}, comp = {1, 2 ,…, N};

 the vector {Ddem}, dem  {1, …, N}, of product required by demand nodes;

 the importance (relevance) of the demand nodes (see Section 2).

The output of the algorithm is represented by the steady state probability distributions of the product delivered to the demand nodes.

The assumptions made are the following:

 only unidirectional relationships are considered;

 the flow in the network is limited by the arcs capacities (i.e., capacitated networks are considered);

 the flow in the network satisfies the flow conservation principle, i.e., no flow in arcs will decrease or increase during flow transmission [START_REF] Ford | Flows in networks[END_REF];

 the capacities of the arcs are statistically independent, i.e., the state of an arc does not depend on the state of other arcs.

For clarification purposes, we describe the procedure with reference to the simple example of Figure 3, where two interconnected systems, S (1) and S (2) , are shown. The N = 5 components are: S1 (1) _S2 (1) , S2 (1) _S3 (1) , S2 (1) _S1 (2) , S1 ( 2) _S2 (2) and S2 ( 2) _S3 (2) . The input component is arc S1 (1) _S2 (1) that serves five demand nodes (i.e., the goals), S2 (1) , S3 (1) , S1 (2) , S2 (2) and S3 (2) , explicitly represented at the top of the diagram. The evaluation is carried out from the bottom to the top of the hierarchy and consists of the following steps:

1) Determine one possible system configuration by sampling the capacity of the arcs, ζcomp,i , i  {1, 2 ,…, NScomp}, comp = {1, …, N}, from the corresponding steady state probability distributions;

2) Identify the minimum arc capacity (mpathdem, dem  {1, …, N}) for each origindestination path: this capacity corresponds to the maximum product that can be delivered to the corresponding demand node dem, dem  {1, …, N}; for example, in Figure 4 the minimum arc capacity for the path from S1 (1) to S3 (1) is the minimum among the capacities of arcs S1 (1) _S2 (1) and S2 (1) _S3 (1) , connecting S1 (1) and S3 (1) ;

3) Set the input (inp) to the network equal to the capacity of the input arc, i.e., inp = ζcomp,i

, where i  {1, 2, …, NScomp} and comp is the index of the input arc (in the example of Figure 3, the input arc is S1 (1) _S2 (1) ); 4) If the input is zero (inp = 0), no product can be delivered to the demand nodes: EPdem = 0 for all dem  {1, …, N}; otherwise, estimate the optimal flows {EPdem} that can be delivered to the demand nodes by the following steps: a. Estimate the vector {EPdem} of optimal flows to the demand nodes taking into account i) the importance of the demand nodes and ii) the minimum capacity of each path (mpathdem, dem  {1, …, N}) that limits the quantity of product that can be delivered to the demand nodes (Figure 5, top). For example, referring to Figure 3, let us consider a proportional importance of the demand nodes (see Section 2) and assume that the ratio of importance of S2 (2) and S3 ( 2) is 0.2 and 0.3, respectively. According to the importance criterion considered and assuming a total input of 100 units, we assign 20 units to S2 (2) and 30 units to S3 (2) . On the contrary, if the minimum capacity of the path to S3 (2) were lower than 30 units (say, mpathS3(2) = 10 units), S3 (2) would receive at most a quantity equal to mpathS3(2) (i.e., 10 units) and the surplus quantity in exceedance (i.e., 20 units) would be distributed to other nodes (see the following steps 4 b. and 4 c.).

b. Initialize to zero an auxiliary variable surp (i.e., surp = 0). This variable is used to quantify the surplus, i.e., the amount of product that cannot be allocated in the network due to arc capacity constraints (i.e., due to the bottlenecks of the infrastructure).

c. Check if the capacities of the links, ζcomp,i, i  {1, 2, …, NScomp}, can support the sum of the estimated optimal products to the corresponding demand nodes, (dem) computed at the previous step 4 a. Such evaluation is performed from the bottom to the top of the diagram. If the sum of the estimated optimal product to the nodes served by a link is higher than its capacity, save the exceeding amount () in the auxiliary variable surp (i.e., surp = surp + ) and compute the optimal partition just for the nodes that are supplied by that link, considering as input the corresponding arc capacity inp = ζcomp,i , where i  {1, 2, …, NScomp} and comp is the link under analysis (Figure 5, middle).

d. Create a "new" graph, where the "new" capacities of all the arcs are updated on the basis of the quantity of product, {EPdem}, that has been effectively allocated at step 4 c. In particular, the arc capacities are reduced by the total quantity of product that they have already supplied to the corresponding demand nodes (Figure 5, bottom).

e. Compute again the minimum arc capacity for each path of the "new" graph (as in step 2) to evaluate the new maximum product that can reach the corresponding demand nodes (Figure 5, bottom).

f. Update the demands {Ddem}, dem  {1, …, N}, reducing them by the quantity {EPdem} that has been already allocated at step 4.c (Figure 5, bottom).

g. Set the input inp equal to the auxiliary variable surp (inp = surp) and repeat step 4 until surp > 0 and the minimum ("new") arc capacity is not zero for at least one path. When one of these conditions is verified, the final vector {EPdem} of the optimal product that can be delivered to the demand nodes is determined (Figure 5, bottom).

The procedure above is repeated a large number of times (e.g., 10000) for many different MCsampled values of the arc capacities and the probability distribution of the product delivered at steady state to each demand node is obtained.

A final remark is in order with respect to the use of Monte Carlo simulation for the propagation of uncertainty of the arc capacities. Real-world complex infrastructure systems include a large number of components, many of which have uncertain (i.e., aleatory, stochastic) performance.

This high-dimensionality constitutes a major challenge for the extensive propagation of (aleatory) uncertainty through the corresponding mathematical models. Classical analytical or numerical schemes are not suitable. On the other hand, stochastic simulation methods (e.g.

Monte Carlo simulation), based on the repeated random sampling of possible component states and the evaluation of the system model for the different states sampled, offer a feasible means: such approaches are, therefore, adopted in the present work.

Concerning, instead, the issue of uncertainty representation, the classical way to address the uncertainty due to randomness is to collect data about the random phenomenon of interest and perform a statistical analysis to identify the probability distribution that best captures the variability of the available data. In a frequentist view, the available data may be interpreted as observable random realizations of an underlying, repeatable probabilistic model (e.g., a probability distribution) representing the stochastic (aleatory) phenomenon of interest which can be approximated with increasing precision by the analyst as the size of the available data set increases [START_REF] Apostolakis | The Concept of Probability in Safety Assessments of Technological Systems[END_REF]. It is worth noting that the procedure proposed is based on several iterative estimations of the vector {EPdem}, obtained by repeating steps 4 a.g. from the bottom to the top of the hierarchy:

in the very first iteration, the system configuration is the one sampled at step 1. and the input product corresponds to the capacity of the input arc; then, at each loop a "new" graph is considered where (i) the new product input value is represented by the surplus (surp), i.e., the amount of product that has not been allocated in the network at the previous iteration, (ii) the "new arc capacities are reduced by the total amount of product they have already supplied at the previous iteration and (iii) the "new" demands are scaled by the quantity already allocated in the previous iteration.

Finally, it is worth mentioning that a drawback of the Hierarchical Graph representation proposed may be represented by its difficult applicability to large networks since all the origindestination paths have to be identified and the bottlenecks of each path have to be spotted out.

To overcome this limitation, we propose to pre-process the infrastructure system by means of a clustering algorithm to reduce the systems dimension by "collapsing" many components in few representative clusters and, then, apply the Hierarchical Graph to the "clustered" infrastructure. The general concepts underlying the pre-processing phase based on clustering is discussed in the following Section 3.2; a discussion on the advantages and limitations of the approach is given in Section 3.3.

Combination of Hierarchical Graph representation and clustering for large-sized critical infrastructures

In order to manage large-sized CIs, it is useful to resort to clustering techniques to reduce the complexity and dimension of the system. Actually, complex systems are characterized by modularity that allows identifying groups of elements highly interconnected within them and sparsely linked to other dense groups in the network [START_REF] Sales-Pardo | Extracting the hierarchical organization of complex systems[END_REF][START_REF] Porter | Communities in Networks[END_REF]]. In addition, several studies show that networks often exhibit hierarchical organization [START_REF] Ravasz | Hierarchical organization in complex networks[END_REF][START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF]. These features lead to combine hierarchical modeling and clustering analysis to represent complex networks.

Cluster analysis aims at identifying groups of "similar behavior" in their data. For illustration purposes, refer to the simple example of Figure 6, left, where the original components of a network, namely S1 (1) , S2 (1) , …, S16 (1) , are shown. According to some features of interest (e.g., proximity), the components can be clustered in groups of "similar characteristics": in the example proposed, four clusters, C1, …, C4, are identified (dotted oval shape in Figure 6, left).

Then, a less complex analysis can be performed on the new fictitious, artificial (i.e., clustered) network, composed just by the identified clusters (Figure 6, right).

The cluster analysis can be carried out at different levels of details: an artificial network with a high number of clusters is closer to the original one and, thus, it is more detailed (i.e., it carries more information) than one with a small number of clusters. The system can be clustered at different levels of details, which allows building a hierarchical 1 clustering representation where the different hierarchical levels correspond to the different levels of detail of the analysis [START_REF] Fang | Unsupervised spectral clustering for hierarchical modelling and criticality analysis of complex networks[END_REF]. In order to reduce the size of the infrastructure under analysis, in this work we have applied the unsupervised spectral clustering algorithm (USCA) and the Fuzzy c-means (FCM) clustering as in [START_REF] Fang | Unsupervised spectral clustering for hierarchical modelling and criticality analysis of complex networks[END_REF], for its simplicity of implementation and its effectiveness in providing satisfactory results. The operative steps of the procedure are not reported here for brevity sake; the interested reader is referred to [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF][START_REF] Fang | Unsupervised spectral clustering for hierarchical modelling and criticality analysis of complex networks[END_REF]. The result of clustering is sensitive to the similarity function which defines the proximity of the nodes in the network; in this respect, the not trivial task is to ensure "meaningful" local neighborhoods [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF][START_REF] Fang | Unsupervised spectral clustering for hierarchical modelling and criticality analysis of complex networks[END_REF]. As similarity measure to perform the clustering, we have adopted the adjacency matrix, i.e., two nodes are similar if they are linked 1 Notice that in this case the term "hierarchical" refers to the level of detail of the clustering and not to the levels of the Hierarchical Graph representation that instead correspond to the number of demands served by a given arc of the network.

directly, otherwise they are not. Different similarity definitions can lead to different cluster partitioning of the network. In this case, the use of the adjacency affinity, which considers only local direct connection information, is justified by the fact that it allows obtaining radial networks for each (clustering) hierarchical level, a condition needed to employ the Hierarchical Graph (see Section 2).

By recursively operating the USCA and the FCM a hierarchical structure of the system can be obtained, where, at the top of the hierarchy, the system is represented by just one node and, at the bottom, by the whole original network. In the middle, each hierarchical level represents a different degree of resolution of the systems (from the top to the bottom it can be seen an increasing quantity of information about the local connectivity) and corresponds to artificial networks that include artificial nodes and links (these last ones are composed by those original network links connectingin parallelthe original nodes in the clusters forming the artificial nodes).

By applying the USCA and the FCM, an artificial (fictitious) network composed by kL clusters,

C1 (L) , …, Ck L (L)
, is produced at each (clustering) hierarchical level L. Notice that the last level of the clustering hierarchy coincides with the real SoS, i.e., the corresponding clusters coincide with the actual/original/real nodes of the SoS. The clustering is performed on the entire network except for the input nodes that are left out (only one generation node is considered in the application of the present work). For illustration purposes, Figure 7 depicts a sketch of the decomposition in five (clustering) hierarchical levels of a SoS with one input node, S1 (1) ; level 1 of the hierarchy is, then, composed by two nodes: the input, S1 (1) , and the rest of the system "condensed" in cluster C1 (1) . The clustering algorithm allows a new analysis at hierarchical level 2 and it decomposes cluster C1 (1) of hierarchical level 1 into two clusters C1 (2) and C2 (2) . At this point, if we want to increase the level of refinement of the analysis we can use the algorithm to further split clusters C1 (2) and C2 (2) . In the example of Figure 7, this results in the decomposition of cluster C1 (2) into three clusters (C1 (3) , C2 (3) and C3 (3) ) and cluster C2 (2) into five clusters (C4 (3) , C5 (3) , C6 (3) , C7 (3) and C8 (3) ). The Hierarchical Graph representation of the decomposed system at level 3 is also shown on the right.

A cluster k is characterized by its demand Dk that is the sum of the demands of the real nodes it contains: for example, cluster C1 (4) of Figure 7 has demand equal to the sum of the demands of nodes S3 (1) , S6 (1) and S7 (1) . For a given clustering hierarchical level L, the quantitative evaluation of the performance of the "artificial" clustered system is carried out as illustrated in Section 3.1 with reference to an indicator that represents the "global" state of the clusters C1 (L) , …, Ck L (L) of level L as a "synthesis" of the real capacity ζcomp,i , i  {1, 2, …, NScomp}, of the arcs, comp = 1, …, N, contained in the cluster itself. Actually, a measure of the cluster state is needed to approximately estimate the quantity of product that a cluster can receive and deliver to other clusters.

To represent the state of a cluster k (i.e., its performance), we consider an indicator idk based on the ratio of the expected capacity of cluster k at current and at nominal (optimal) conditions as follows:

𝑖𝑑 𝑘 = ∑ 𝑤 𝑐𝑜𝑚𝑝 * 𝜁 𝑐𝑜𝑚𝑝,𝑖 𝑛 𝑘 𝑐𝑜𝑚𝑝=1 ∑ 𝑤 𝑐𝑜𝑚𝑝 * 𝜁 𝑐𝑜𝑚𝑝,𝑁𝑆 𝑐𝑜𝑚𝑝 𝑛 𝑘 𝑐𝑜𝑚𝑝=1 , (1)
where comp indicates the component (arc) of the original network, nk is the number of arcs inside cluster k, 𝜁 𝑐𝑜𝑚𝑝,𝑖 , i  {1, 2, …, NScomp}, is the current (i.e., actual / sampled) state of the component comp, 𝜁 𝑐𝑜𝑚𝑝,𝑁𝑆 𝑐𝑜𝑚𝑝 is the maximum capacity of the arc comp and wi is the weight associated to the capacity of the arc comp. The weight 𝑤 𝑐𝑜𝑚𝑝 is computed as the ratio between the capacity of the arc comp and the sum of the maximum capacities of all the arcs of the network, i.e., 𝑤 𝑐𝑜𝑚𝑝 = 𝜁 𝑐𝑜𝑚𝑝,𝑖 ∑ 𝜁 𝑐𝑜𝑚𝑝,𝑁𝑆 𝑐𝑜𝑚𝑝 𝑁 𝑐𝑜𝑚𝑝=1 ⁄ , i  {1, 2, …, NScomp} and gives an idea of the weight of the arc in the entire network. The index of the cluster state (idk), takes value between 0 and 1.

Notice that the state of a cluster affects the cluster itself and the connected clusters, since the cluster is both a fictitious load node (which should provide itself with the required amount of product) and a fictitious transmission node (which should transmit the product to the other connected clusters). The top of Figure 8 shows two clusters, C1 and C2, supplied by the input source S1 (1) : cluster C2 is both a load and a transmission node, since on one side it contains five demand nodes (S2 (1) , S3 (1) , S4 (1) , S5 (1) and S6 (1) in Figure 8, bottom) and on the other side it is required to transmit the product to cluster C2. In particular, the product from input source S1 (1) has to pass through two arcs (S2 (1) _S5 (1) and S5 (1) _S6 (1) ) contained in C1 to reach cluster C2: if their capacities decrease, then the flow to nodes S5 (1) and S6 (1) (i.e., to the cluster C1) and to nodes S7 (1) , S8 (1) and S9 (1) (i.e., to the cluster C2) is reduced.

Figure 8: Top: artificial system composed by two clusters C1 and C2 supplied by one input node S1 (1) . Bottom: illustration of the real nodes inside the fictitious clusters: two arcs of C1 are needed to supply C2.

Thus, when the "capacity" of a cluster decreases, the consequence is twofold: the cluster cannot satisfy itself (i.e., the demands of the demand nodes it contains) and also the connected clusters.

In order to take into account the "twofold" reduction of performance, we "artificially reduce" the amount of product that can be given to the cluster itself and that can be delivered to the connected clusters, by multiplying by the indicator of the state of the cluster, idk i) the maximum demand Dk that it requires and ii) the maximum capacities, 𝜁 𝑐𝑜𝑚𝑝,𝑁𝑆 𝑐𝑜𝑚𝑝 , of the arcs (comp) that link the output clusters, respectively.

Discussion on the Monte Carlo simulation and Hierarchical Graph modeling for critical infrastructures analysis

The Hierarchical Graph representation supports both a qualitative and quantitative analysis of CIs. The former allows identifying critical arcs; the latter allows determining the partitioning of the flow in the network and computing a measure of system performance, e.g., the steady state probability distribution of the product delivered to demand nodes. However, the representation does not make evident how the partition of the flow is carried out, since the algorithm used for its computation is recursive.

The Hierarchical Graph representation can aid capturing the structural and dynamic complexities of CIs as follows. It can represent components of different technological domains (e.g., physical, hard componentslike pipelinesand soft componentslike information and communication systems) and their interconnections by dependences, distributed over large geographic areas. Different types of service and product demands can be included, e.g., supply of power to critical buildings (hospitals, industrial plant) and residential areas, supply of water for daily use or for emergency situations (fire-fighting), etc. However, a limitation is that the interdependences (bidirectional relationships) cannot be represented, since the representation is based only on unidirectional flow. Also, the Hierarchical Graph becomes large when the scale and dimensionality of connectivity is large, and the identification of all the origin-destination paths may become computationally prohibitive. To overcome this limitation, the clustering technique is employed to reduce the size of the system, obtaining a smaller radial artificial (clustered) network, i.e., without loops. This allows a preliminary analysis of system robustness to be performed in short time, followed by a more detailed analysis on the most relevant parts of the system.

With respect to the dynamic complexity, the quantitative analysis of the system robustness supported by the Hierarchical Graph gives indications on how much the system is able to adapt to failure scenarios for continuing to efficiently deliver the products and services through the network, according to different criteria of importance. It is to be noted that the Hierarchical Graph is not well suited to represent growing systems, since the addition of a new component to the system methods may change its entire structure (i.e., it lacks of flexibility).

The overall approach is based on the assumption that capacitated networks (i.e., networks where the flows of product in the links are limited by arcs capacities) provide an appropriate representation of infrastructure performance, e.g., robustness. This assumption does not represent a limitation for the analysis and, indeed, several works in the scientific literature are based on the same hypothesis. Actually, in many real flow networks such as electric power networks, gas distribution networks, computer networks, urban traffic networks, logistics networks, each arc has several possible capacities or states due to failures, partial failures, maintenance, etc. Works related to capacitated networks can be found in [START_REF] Lin | On reliability evaluation of a capacitated-flow network in terms of minimal pathsets[END_REF][START_REF] Ramirez-Marquez | Multi-state component criticality analysis for reliability improvement in multi-state systems[END_REF][START_REF] Yeh | A simple minimal path method for estimating the weighted multicommodity multistate unreliable networks reliability[END_REF][START_REF] Ramirez-Marquez | Stochastic network interdiction optimization via capacitated network reliability modeling and probabilistic solution discovery[END_REF][START_REF] Lin | On transmission time through k minimal paths of a capacitated-flow network[END_REF][START_REF] Lin | On transmission time through k minimal paths of a capacitated-flow network[END_REF]Yeh, 2011] [Nozick et al., 2005].

Table 1 summarizes the advantages and limitations of the approach."

Table 1: Advantages and limitations of the Hierarchical Graph representation

Advantages

 Both qualitative analysis (identification of critical elements) and quantitative analysis (computation of a measure of system performance) can be performed.

 The complexity of CIs can be handled.

 Different types of demands (daily water consumption, fire-fighting consumption, power consumption, etc.), with different importance (priority), can be taken into account.  A clustering algorithm can be employed to reduce the model dimension.

 A multi-state model can be adopted for a more realistic representation of the system states. Limitations  Interdependences (bidirectional relationships) cannot be included in the representation.  Difficult applicability to large networks (but clustering can reduce this problem).

 The representation itself shows the physical constraints that affect the product and service distribution, but it does not illustrate how the partition of the products and service is carried out, since the algorithm used for its computation is recursive.  The representation is not flexible, i.e., the addition of a new component may change the entire structure.

Applications

In this Section, we apply the proposed Hierarchical Graph representation to evaluate the robustness of two case studies (hereafter referred to as "A" and "B"): case study A (Section 4.1) consists of two interdependent infrastructures (gas and electric power networks) and a supervisory control and data acquisition (SCADA) system connected to the gas network; case study B (Section 4.2) considers an electric power distribution network adapted from the IEEE 123 node test feeders. As mentioned before, the robustness addressed in this paper refers to the capability of the CIs to resist to failures or partial failures of the CIs components (arcs) assuring the required level (or a high level) of supply of goods or services. It is measured as the steadystate probability distributions of the product (e.g., gas and/or electricity) delivered to the demand nodes.

In both cases studies, we adopt a multi-state model to account for the different degrees of degradation of the components and we describe the state transition processes by Markov and semi-Markov modeling.

Case study A

The case study is taken from [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF] and deals with two interconnected infrastructures, i.e., a natural gas distribution network and an electricity generation/distribution network (Figure 9, solid and dashed lines, respectively). The gas distribution network is supported by a SCADA system (Figure 9, dotted lines). The objective of this interconnected system of systems (SoS) is to provide the necessary amount of gas and electricity (hereafter also called "product") to four demand nodes (end-nodes), namely D1 and D2 (gas), and L1 and L2 (electricity).

Figure 9: Interdependent gas (solid lines) and electric (dashed lines) infrastructures and SCADA system (dotted lines) [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF]]. The possible states (i.e., capacities) of the arcs are given in square brackets; the quantities of product demanded by end-nodes D1, D2, L1, L2 are reported in bold.

The gas distribution network, supplied by two sources of gas (namely, S1 and S2, connected to the network by arcs S1_DS1 and S2_DS2, respectively), provides gas to the end-nodes D1 and D2 and to two nodes of the electricity network (E1 and E2). Once the gas enters into nodes E1 and E2, it is transformed into electrical energy that flows through arcs E1_G1 and E2_G2

(representing the electric power generation stations) to supply the end-nodes of electricity (L1 and L2); notice that demand L2 can be supplied by both electrical generations E1_G1 and

E2_G2.

The assumption is made that the gas-electricity transformation occurs with a constant coefficient, i.e., 100 cu. ft. of natural gas produces 1 MWh of electricity [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF].

A SCADA system controls the gas flow through arcs a_b, b_c, c_d and d_e. It is assumed that: i) the SCADA has two core subsystems controlling different sets of arcs (in particular, the first one -SUB1refers to links a_b and b_c, whereas the second one -SUB2controls arcs c_d and d_e); ii) the SCADA is always provided with electric power [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF].

The capacities of the arcs of the gas and electricity networks (determining the maximum flows of gas or electricity supported by the arcs) can be deterministic (i.e., fixed constant values) or stochastic (i.e., randomly evolving in time) (Figure 9, values in the square brackets). The stochastic capacities give rise to a multi-state model that reflects the possibly different degrees of degradation of the arcs. On the contrary, the SCADA system state is defined by a binary random variable, whose values 1 and 0 represent its complete and partial functioning, respectively. For example, when the state of the SCADA subsystem SUB1 (controlling arcs a_b and b_c) is 0, the capacity of these arcs decreases because of the incorrect information provided by the SCADA subsystem (even if the arcs are not subject to a direct damage). On the basis of the two states of the SCADA subsystems, two different vectors of capacities are identified for each arc a_b, b_c, c_d and d_e: as illustrated in Figure 9, the first vector is used when the corresponding SCADA subsystem is in state 0, whereas the second one is employed when the SCADA subsystem is in state 1.

Changes in the arc capacities are due to random degradations or recovery actions. The state transitions over time are modeled by Markov and semi-Markov processes as in [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF]. Semi-Markov processes are adopted to model the evolution of the capacities of the gas supply links (S1_DS1 and S2_DS2), whereas Markov processes are used for all the others arcs.

Both Markov and semi-Markov processes for a generic component 'comp' are defined by a transition probability matrix

} ..., , 2 , 1 , : { P comp ij S j i p comp  
, where pij is the one-step probability of transition from state i to state j. In addition, the semi-Markov processes are characterized by continuous probability distributions for the holding time T ij comp, i.e., for the time of residence of a component in state i before performing a transition to state j.

The steady-state probability vectors for a generic component 'comp' described by a Markov process is computed as [START_REF] Zio | Computational methods for reliability and risk analysis[END_REF]. For a semi-Markov process, this equation is weighted by the expected time of residence τ i in a given state i before performing a transition, i.e., as [START_REF] Barry | Stochastic modeling: analysis and simulation[END_REF].

comp comp comp P    
i i comp i comp      , ,     comp S j j j comp 1 , /  , i = 1, …, NScomp
Figure 10 reports the steady-state probability vectors of the arcs of the system of Figure 9 assuming the state transition probabilities given in [START_REF] Nozick | Assessing the performance of interdependent infrastructures and optimising investments[END_REF].

Hierarchical Graph representation for the system of case study A

In Figure 11, the Hierarchical Graph of the system of case study A (Section 4.1) is illustrated.

The injection of product (i.e., gas) in the SoS is made through arcs S1_DS1 and S2_DS2, located at the bottom of the diagram (Section 2). Since both arcs enter in node a, also the following links DS1_a and DS2_a are considered part of the inputs and reported at the bottom of the hierarchy. Four demand nodes, i.e., D1, D2 (gas) and L1, L2 (electricity), represent the goals of the analysis and they are explicitly located at the top of the diagram. The graph presents four hierarchical levels: in level 4, arc a_b, is reported since it supplies all the four demand nodes; in level 3, arc b_c is depicted, since it serves three demand nodes (i.e., L1, L2 and D2); in level 2, arcs c_E1, E1_G1, c_d and d_e are considered, since they supply two demand nodes: in particular, arcs c_E1 and E1_G1 supply L1 and L2, whereas arcs c_d and d_e serve L2 and D2;

in level 1, there are the remaining arcs that are related just to one demand node: for example, e_E2 serves only node L2. The influence of the SCADA subsystem SUB1 on the arcs a_b and b_c and of the SCADA subsystem SUB2 on the arcs c_d and d_e is illustrated in the trapezoidal frames under the corresponding arcs. For illustration purposes, three different importance criteria for the demand nodes are considered (see Section 2):

 sequential importance: the demand nodes are ranked on the basis of the proximity to the sources: 1) D1 (the most important), 2) L1, 3) L2 and 4) D2 (the least important).

 proportional importance: the demand nodes are satisfied on the basis of their demands (the nodes that require more product are given higher priority). Since D1 and D2 require 100 and 80 cu. ft. of gas and L1 and L2 need 500 and 400 MWh (equivalent to 50 and 40 cu. ft.), the resulting importance ranking of the demand nodes is: 1) D1 (the most important), 2) D2, 3) L1 and 4) L2 (the least important).

 equal importance: the product is divided equally among four demand nodes.

Results of case study A

The robustness of the case study A is given by steady state probability distributions of the product delivered to the four demand nodes D1, L1, L2 and D2, by applying the operative steps illustrated in Section 3.1. From these distributions two quantities can be analyzed as illustrated in Table 2: (i) Steady state probabilities of delivering the (maximum, optimal) required product to the demand nodes (top);

(ii) Steady state probability of delivering a quantity of product exceeding 90% of the corresponding demands (bottom). As expected, in the case of sequential importance of the demand nodes, D1 is the demand node most satisfied, whereas D2 is the least served: the corresponding probabilities of being completely satisfied are 0.9927 and 0.7526, respectively (Table 2, top). Differently, in the ranking produced by the proportional importance criterion, D2 is more important than L1 and L2. Thus, the probability of delivering the required product to D2 should increase, whereas the probability of satisfying the last two demand nodes should decrease. In fact, the steady state probabilities of delivering the required maximum product to L1 and L2 decrease (Table 2, top).

On the contrary, the probability P(D2 = 80 [1000 cu. ft.]) of delivering the maximum product to D2 remains almost the same, since (i) the path needed to reach D2 is affected by the uncertainty on the capacity of many arcs and (ii) D2 is not the most important demand node.

However, an effect of the proportional importance criterion on D2 can be seen by analyzing the steady state probability P(D2 > 72 [1000 cu. ft.]) that the product given to D2 exceeds 90% of its demand: actually, this value increases considerably from 0.7526 of the previous case (sequential importance, Table 2, top) to 0.9160 (proportional importance, Table 2, bottom).

Finally, the criterion of equal importance, turns out to give preference to the nodes that have the lowest demands (this is expected since a lower demand has higher probability to be satisfied by an equal partition of the product). Actually, in this case the steady state probabilities of supplying L1, L2 and D2 at the optimum level increase with respect to the case of proportional importance (from 0.7563 to 0.9306 for L1, from 0.7563 to 0.9678 for L2 and from 0.7568 to 0.9063 for D2, see Changing the importance criterion induces changes in the ranking of the demand nodes that receive less product. However, it can be noticed that the demand node D2 is the most difficult to supply to, even in the case of proportional importance where D2 is the second most important demand node. The highest probability of supply of D2 is 0.9160 while the others demand nodes can reach a probability of supply higher than 0.98. Then, to improve the robustness of the system the path from the sources to the demand node D2 should be strengthened.

Case study B

Figure 12 shows the electric power distribution network here considered, as adapted from the IEEE 123 nodes test feeder [IEEE, 2000] in the sense that regulators, capacitors, switches and feeders with length equals to zero are neglected. With these simplifications, the network is composed of 114 nodes: 1 generation point (node 115) and 113 load/transmission nodes. Node 61 of the original IEEE 123 node test feeders is missing here, since after the removal of switches and transformers it turns out to be an end node with load equal to zero. The arcs (i.e., the feeders) connect different nodes and distribute the power through the network. In the analysis by Hierarchical Graph, we focus on the arcs (and not on the nodes), hereafter also called "components"; thus, the total number N of components is 113.

Results of case study B

In order to evaluate the robustness of the electric power distribution network of Section 4.2, an unsupervised spectral clustering algorithm (Section 3.2) is first applied to reduce the size of the system and, then, Monte Carlo simulation is performed at a given clustering level in combination with the Hierarchical Graph representation (Section 3.1). The outcome of the procedure is represented by the steady state probability distribution of the electricity delivered at the demand nodes (or clusters). In this case, we consider proportional importance of the clusters/demand nodes (see Section 2).

As illustrated in Section 4.2.1, five hierarchical clustering levels are identified by the unsupervised spectral clustering algorithm. In the following, for the sake of brevity the results obtained from the analysis of the "fictitious" clustered networks are mainly given for the levels 3 and 4 of the clustering hierarchy and are compared to those obtained from the analysis of the original (i.e., not clustered) system.

In order to perform a fair comparison, the values of the product delivered to the nodes obtained analyzing the real network are "grouped" on the basis of the clustered structure of the corresponding "fictitious" network: for example, at hierarchical clustering level 4 the products delivered to nodes 9, 10, 11 and 14 have to be summed to allow the comparison with the corresponding cluster 6 of the "fictitious" network (Figure 15). The comparison is, then, made by means of the expected product not supplied idNSk to cluster k, computed as:

𝑖𝑑𝑁𝑆 𝑘 = 𝐷 𝑘 -𝜇 𝑘 𝐷 𝑘 , ( 2
)
where Dk is the demand of cluster k (equal to the sum of the demands of the nodes of the original network contained in cluster k) and µk is the mean value of the product delivered to the cluster k at steady state.

We have performed the comparison between the performances of the clustered and real networks in three cases:

 Case 1: all the arcs are characterized by the steady state probability vector Π 1 of Figure 16 left; in this case, the system is expected to have a satisfactory performance as all the arcs remain in their best state (i.e., state 4) with very high probability (i.e., 0.9982).

 Case 2: all the arcs are characterized by the steady state probability vector Π 2 of Figure 16, middle; in this case the arcs turn into state 3 (i.e., a state of partially reduced functionality) with high probability (i.e., around 0.5) so that the global performance of the system is not expected to be very high;

 Case 3: same as case 2, except for the fact that the arcs contained in clusters 8 and 12 at hierarchical clustering level 4 are assumed to be more prone to degradation than the other arcs and then characterized by the steady state probability vector Π 3 of Figure 16, right; in particular, they degrade to their worse states (e.g., state 2) with high probability. 12 for cases 1 (left), 2 (middle), 3 (right).

In the following, comments on the results obtained for the three cases above are given with reference to Tables 34where the clusters at hierarchical levels 3 and 4, respectively, are ranked in ascending order based on the product not supplied (idNSk) (i.e., the mean value of the product supplied to clusters at the top of the ranking is closer to their demands than that of the clusters at the bottom). Tables 34, right, illustrate also the ranking of the clusters based on their distance (measured in terms of number of arcs) from the source point: for example, at level 3 cluster 8 is connected to the input source by four arcs (i.e., 115_1, 53_54, 60_67 and 101_105) (see Figure 14). The reference to the distance of a demand node from the source is motivated by the fact that in general the higher the distance from the source, the higher the probability of reduced product supply.

Case 1

The results obtained from simulations on the clustered and real networks are very similar at all hierarchical clustering levels. In particular, the higher is the hierarchical level (i.e., the finer is the clustering), the closer are the results of the fictitious network to those of the real network.

The probability distribution functions (pdfs) of the product delivered to the clusters of the fictitious and real networks at steady state present the same modal values (i.e., the same peaks) at hierarchical levels 2, 3 and 4 and mean values very close to each other at all hierarchical clustering levels. For illustration purposes, the mean and modal values of the pdf of the product delivered at steady state to cluster 1 at hierarchical levels 2, 3 and 4 obtained from the analysis of the real and clustered networks is given in Table 5, together with their percentage difference. by a reduction of performance and those connected to these last ones), followed by clusters 14, 15, 22, 28, 29 (i.e., those that are farther from the generation point). Thus, it can be concluded that in those cases where the performance of some arcs is biased towards very low values, a finer level of analysis is needed to capture the global behavior of the network.

Finally, the computational cost of the simulations depends on the configuration of the system: the higher the capacities of the arcs, the faster the process to allocate the product by Hierarchical

Graph. However, in all the three cases explored, it is evident a considerable reduction in the computational time when the level of detail of the analysis is reduced (i.e., when the network is clustered). In Table 8, the computational times needed to perform the simulations on the clustered networks at hierarchical levels 2, 3 and 4 and on the original network are given for cases 1, 2 and 3 detailed above. In conclusion, a robustness analysis has been carried out at different (hierarchical) clustering levels by Hierarchical Graph and Monte Carlo simulation for an electric power distribution network adapted from the IEEE 123 node test feeders. The results have shown that it is possible to apply the approach resorting to clustering techniques to reduce the computational time of the simulations since:

-The modal and mean values of the distribution of the pdfs of the product delivered at steady state to the clusters are similar for the analysis of the real and clustered networks;

-The simulations of the fictitious and original networks produce a similar rankings of the product-not-supplied index: it can be noticed a strong correspondence of "groups" of clusters at the tails of the ranking.

However, attention should be paid when some arcs inside the clusters are more prone to degradation than other arcs. In this case, the simulation carried out on a clustered network characterized by a low level of detail (i.e., by a low number of clusters) may not always capture this reduction of performance: a finer clustering level is, then, necessary.

These analyses carried out on clustered networks provide approximate results that are useful in a decision making process for a preliminary design step, reducing considerably the computational cost of the simulation.

Conclusions

In this paper, we have proposed a new representation technique, i.e., the Hierarchical Graph, to analyze the performance of interconnected critical infrastructures (CIs) under a multi-state system-of-systems (SoS) framework. In particular, the robustness of the SoS, i.e., the capability of the CIs to resist to failures or partial failures assuring the required level (or a high level) of supply of goods or services, has been evaluated in terms of the product delivered at steady state to the demand nodes of the CIs.

The Hierarchical Graph representation supports both a qualitative and quantitative analysis to evaluate the robustness of CIs. Actually, critical arcs, which interruption or degradation affects more demand nodes, can be visualized and the product delivered at steady state to the demand nodes of the CIs can be quantified, while accounting for:

-the importance (relevance, priority) of the demand nodes;

-the quantity of product required by the demand nodes;

-the state of the arcs of the networks (possibly affected by failures and partial failures).

The quantification algorithm based on the Hierarchical Graph follows a logical analysis, while not resorting to a flow model.

In the paper, we have applied the approach on two case studies. First, we have analyzed a smallsized SoS composed by two interconnected CIs (gas and electricity networks) and a supervisory control and data acquisition (SCADA) system, and we have evaluated its robustness by Monte Carlo (MC) simulation considering different importance criteria (sequential, proportional, and equal) for the demand nodes. We have shown that the Hierarchical Graph representation can support this kind of analyses that are useful for decision makers to understand margins of improvement of the SoS to optimize the delivery of product to the demand nodes by changing their importance. The small size of this case study has allowed understanding the details of Hierarchical Graph modeling.

Then, we have evaluated a moderately large-sized power distribution network by adopting a hierarchical clustering algorithm to analyze the SoS at different levels of detail and simplify the Hierarchical Graph representation. In this case, only a proportional importance of the demand nodes has been considered. The results have shown that the Hierarchical Graph can be adopted together with hierarchical clustering algorithms to provide approximate results by analyzing clustered networks instead of the entire large-sized, real network. This can be useful in a first preliminary phase of design of SoS, in order to have satisfactory, physically coherent results with relatively low computational cost.

The work presented is methodological in nature and intended to explore the feasibility of the proposed approach, while understanding its functioning and potential benefits. To this aim, simplifying assumptions have been made in the numerical case studies analyzed. In spite of some limitations due to the assumptions made (e.g., the use of exclusively unidirectional relationships, the prediction of possibly imprecise results when the clustering algorithm is employed), the proposed approach has shown some interesting benefits that make it potentially attractive for the analysis of CIs robustness, e.g., the fact that qualitative and quantitative analyses can be combined and multistate modeling can be adopted for more refined insights on system behavior.

Figure 1 :

 1 Figure 1: Top: graph of the components of the system of systems; the links represent the exchange of physical product (solid lines) and influence/support relationships (dotted lined). Bottom: corresponding Hierarchical Graph; LV: Level.

Figure 3 :

 3 Figure 3: Hierarchical Graph of a generic example taken as reference to illustrate the algorithm; LV: Level.

Figure 4 :

 4 Figure 4: Exemplification of step 2 of the algorithm with respect to the example proposed in Figure 3.

Figure 5 :

 5 Figure 5: Exemplification of step 4 of the algorithm with respect to the example proposed in Figure 3.

Figure 6 :

 6 Figure 6: Exemplification of the clustering procedure.

Figure 7 :

 7 Figure 7: Left: sketch of the decomposition of a system in five hierarchical levels (L), where the last one (L = 5) coincides with the actual nodes of the system; right: Hierarchical Graph of the corresponding hierarchical level 3; LV: Level of the Hierarchical Graph.

Figure 10 :

 10 Figure 10: Steady state probability vectors for the states of the components of the case study A (comp = S1_DS1, S2_DS2, a_b, b_c, c_d, d_e, E2_G2, E1_G1, SCADA).

Figure 11 :

 11 Figure 11: Hierarchical Graph of the SoS depicted in Figure 9 (case study A); LV: Level.

Figure 12 :

 12 Figure 12: IEEE 123 node test feeders adapted to the purposes of the present analysis.

  Four states, i = 1, …, 4, characterize the capacity (ζcomp,i, comp = 1, …, N) of the arcs: the first one (i = 1) represents the worst state and corresponds to the complete failure of the link, i.e., no

Figure 13 :

 13 Figure 13: IEEE network clustered at hierarchical level 2 (left) together with the corresponding Hierarchical Graph representation (right); LV: Level.

Figure 14 :

 14 Figure 14: IEEE network clustered at hierarchical level 3 (left) together with the corresponding Hierarchical Graph representation (right); LV: Level.

Figure 15 :

 15 Figure 15: IEEE network clustered at hierarchical level 4.

Table 2 :

 2 Steady state probabilities of (i) delivering the (maximum, optimal) required product to the demand nodes (top) and (ii) delivering a quantity of product exceeding 90% of the corresponding demands (bottom).

  Table 2 top); the probabilities to serve D1 remain almost the same (around 0.82, see Table 2 top).

	Sequential	0.9927	0.9867	0.9723	0.7526
	Proportional	0.8195	0.7563	0.7563	0.7568
	Equal	0.8205	0.9306	0.9678	0.9063
				P(L2 > 36	P(D2 > 72
				[10 MWh])	[1000 cu.
	Sequential	0.9927	0.9867	0.9723	0.7526
	Proportional	0.9778	0.9155	0.9507	0.9160
	Equal	0.9717	0.9482	0.9816	0.9063

  Figure 16: Steady state probabilities of the four levels of capacity of the arcs of the network in Figure

		Steady state		Steady state		Steady state
		probability		probability		probability
	1	0.0006	1	0.0005	1	0.0291
	Π 1 = 2	0.0003	Π 2 = 2	0.0010	Π 3 = 2	0.0971
	3	0.0009	3	0.4995	3	0.8651
	4	0.9982	4	0.4990	4	0.0087

Table 3 :

 3 Product not supplied, idNSk, to the cluster k for the real and clustered networks at hierarchical level 3, with respect to cases 1, 2 and 3. Right: ranking (in ascending order) of the clusters based on their distance from the input source.

	Level 3
			Case 1				Case 2				Case 3		Distance -
	REAL		CLUSTERED	REAL		CLUSTERED	REAL		CLUSTERED	ranking
	idNS k	k	idNS k	k	idNS k	k	idNS k	k	idNS k	k	idNS k	k	dist k	k
	0.005	1	0.002	1	0.298	1	0.250	1	0.264	1	0.245	1	1	1
	0.011	2	0.004	2	0.338	5	0.298	5	0.338	5	0.296	5	2	2
	0.011	3	0.004	3	0.338	6	0.304	2	0.340	6	0.324	7	3	2
	0.013	5	0.004	5	0.339	2	0.305	3	0.341	7	0.325	2	5	2
	0.013	7	0.004	7	0.339	3	0.322	7	0.345	8	0.332	3	7	3
	0.014	6	0.008	4	0.341	7	0.332	6	0.346	4	0.332	6	6	4
	0.018	8	0.008	6	0.342	4	0.332	8	0.396	2	0.332	8	8	4
	0.019	4	0.009	8	0.345	8	0.338	4	0.412	3	0.335	4	4	5

Table 4 :

 4 Product not supplied, idNSk, to the cluster k for the real and clustered networks at hierarchical level 4, with respect to cases 1, 2 and 3. Right: ranking (in ascending order) of the clusters based on their distance from the input source.

Table 8 :

 8 Computational time [s] of 10000 Monte Carlo simulations for the analysis of the real and clustered networks at hierarchical levels 2, 3, 4 for cases 1, 2 and 3 considered.

			Case 1		Case 2		Case 3
	Level	REAL	CLUSTERED	REAL	CLUSTERED	REAL	CLUSTERED
	2	454	11	1304	11	1392	50
	3	454	33	1304	40	1392	78
	4	454	110	1304	201	1392	256

product can flow through it (ζcomp,1 = 0, comp = 1, …, N); the fourth one (i = 4) corresponds to the best state, i.e., to its nominal designed capacity. To obtain a reasonable value for the nominal capacity of the arcs useful for the purposes of the present analysis, we have solved the DC power flow equations (DCPF) that provide the physical flows of electricity through the lines given i) the total power injected in the network (assumed equal to the sum of all the powers required by the demand nodes), ii) the network topology (i.e., the adjacency matrix) and iii) the reactance of all the arcs [START_REF] Mccalley | The DC Power Flow Equations 1.0 Introduction[END_REF]; the nominal capacity of a link is set equal to the optimal power flow through that link.

In the Appendix, the demands of the network nodes provided by the IEEE database and the maximal capacity ζcomp,4 of the arcs (comp = 1, …, N) obtained by the DC power flow equations, are given.

Changes in the arc capacities are due to random failures (as in the previous case study A): the state transitions over time are modeled as Markov processes.

Hierarchical Graph representation and clustering for case study B

An unsupervised spectral clustering algorithm (Section 3.2) is applied to the IEEE electric power distribution network of interest and identifies five hierarchical clustering levels: level 1 is composed by the generation node and one cluster representative of the whole network, whereas level 5 coincides with the original network. In Figures 131415, the clusters generated at levels 2, 3 and 4 are reported, respectively: the corresponding Hierarchical Graph representation to levels 2 and 3 is also given for illustration purposes. Notice that level 2 is characterized by two clusters, level 3 by 8 clusters and level 4 by 29 clusters (besides the generation point): these clustered representations correspond to different possible "levels of detail" that the analyst may choose to study the network. The values of the product not supplied idNSk are in general very low (i.e., lower than 0.025), so that the presence of less supplied cluster(s) is not evident. In general, the ranking produced by the clustered and real networks are similar (Tables 34, case 1). Differences in the ranking can be found at clustering level 4, but they can be neglected given the very low values of the corresponding product not supplied index. Also a correspondence between the rankings of "groups" of clusters that have similar product not supplied index can be found: for example, the clusters closer to the generation source, e.g., 4, 5, 13, 18, 19, are slightly more supplied than those farther, e.g., 14,15,16,22,27,28,29 for both the fictitious and real systems, confirming the physical coherence of the approach.

Case 2

In general, the results obtained from simulations on the fictitious and real networks are less similar to each other than in case 1, where the arc state capacities present less variability.

However, the trends in the network behavior that have been pointed out in the previous case can still be identified.

The pdfs of the product delivered to the clusters of the fictitious and real networks at steady state present (almost) the same peak values at hierarchical levels 2, 3 and 4. Also the mean values are very close to each other at hierarchical levels 3 and 4 (actually, the maximum percentage difference is 6.9 % for cluster 1 at level 3 and 3.9 % for cluster 5 at level 4); instead, at level 2 the difference is larger, i.e., 10.7 % for cluster 1. Thus, as expected, the higher the hierarchical level (i.e., the finer the clustering), the closer the results of the artificial network to those of the real network. For illustration purposes, the mean and modal values of the pdf of the product delivered at steady state to cluster 1 at hierarchical levels 2, 3 and 4 obtained from the analysis of the real and clustered networks is given in Table 6, together with their percentage difference. With respect to the previous case 1, the probability distributions present higher variability; in particular, this effect is more evident in the results obtained on the artificial networks. For example, referring to cluster 1, the variances of the distributions of the product delivered are 3514, 725 and 25 at levels 2, 3, 4, respectively, for the original network, whereas they are 25283, 3199 and 35, respectively, for the clustered network. This can be due to the presence of the clusters and to the definition of their performance (Section 3.2). Actually, the state of a cluster is represented by indicator idk (eq. 1) that tries to capture and synthetize the main features of the nodes inside the cluster itself. By so doing, the (detrimental) effect of the degradation of an arc is "spread" through the entire cluster instead of having a contained local impact only on the physically connected nodes. This leads to an increase in the variability of the performance of the network: as expected, such variability decreases with the level of detail of the analysis, i.e., with the reduction of the cluster size.

The ranking based on the product not supplied (Tables 34, case 2) is similar for the groups of clusters at the top and at the bottom (i.e., at the tails) of the ranking, so that it is possible to highlight those that are more supplied (e.g., clusters 1 and 5 at level 3, and clusters 1, 2, 3, 4, 5, 6 at level 4) with respect to those who are less supplied (e.g., clusters 4 and 8 at level 3, and 14, 28, 29 at level 4, see Table 4, case 2). As for case 1, the tails of the ranking reflect the distance of the clusters from the generation point: at the top we find the elements closer to the input source, whereas at the bottom those that are farther (see Tables 34, right).

Case 3

The pdfs of the product delivered to the clusters of the artificial and real networks at steady state present almost the same modal (peak) values at hierarchical levels 2, 3 and 4. Instead, the mean values are farther from each other than in the previous case 2: for example, at hierarchical level 3 the maximum percentage difference is 14.5 % for cluster 3, whereas at level 4 it is 7.9 % for cluster 10. At level 2 the maximal percentage difference is larger, i.e., 16 % for cluster 1.

Again, the higher the hierarchical level, the higher the similarity of the results produced by the artificial network and by the real network. For illustration purposes, the mean and modal values of the pdf of the product delivered at steady state to cluster 1 at hierarchical levels 2, 3 and 4 obtained from the analysis of the real and clustered networks is given in Table 7, together with their percentage difference. In this case, the reduction of performance of clusters 8 and 12 at level 4 strongly influences the ranking performed according to the product not supplied (idNSk) indicator.

The simulations carried out on the artificial networks at hierarchical level 3 cannot capture the reduction in the performances of the arcs inside clusters 8 and 12 at level 4. Actually, according to the simulation performed on the original network, at level 3 clusters 2 and 3 should be the most impacted, since they are those that contain the arcs with the worst capacities; those are followed by clusters 4 and 8 that are the farthest from the input source. On the contrary, the simulations of the fictitious networks at level 3 identifies clusters 4 and 8 as the least supplied (