¹

1 Introduction↩︎

Networks offer a powerful tool for modeling and analyzing diverse systems with complex interactions, such as ecological, social, biological, and technological systems [1]–[4]. Examples include modeling interdependent critical infrastructure (ICI) systems (e.g., coupled water and power systems) as multi-layer networks to analyze their vulnerability and resilience to disruptive events [5]. By definition, interdependency refers to bidirectional dependencies between infrastructure networks through which individual networks interact, influencing overall system operations [6], [7]. There are primarily four categories of interdependencies [7], [8]: (i) physical, (ii) cyber, (iii) geographic, and (iv) logical. This study examines physical interdependencies. Additionally, multi-layer networks in this context refer to different layers inside each critical infrastructure network, including supply, transmission, and demand layers. Ideally, performing such network-level assessments requires data from real infrastructure networks with complete information on the network topology, flow, and interdependencies. However, data on the topology of real-world infrastructure networks is often unavailable either due to security concerns or decentralized operations of different infrastructure sectors [9]. As a result, infrastructure network performance has been evaluated using model-driven techniques relying on assumptions of the existence of interdependencies. Recent research advances have explored synthetic infrastructure data generation to overcome real network data challenges [6], [10], [11]. In this study, we propose a data-driven technique founded in network reconstruction methods to infer the structure of complex interdependent infrastructure networks based on cascading failure data containing independent cascading failure scenarios (See Definition 1 in Section 2 for mathematical details). A cascading failure occurs when an initial failure in an interconnected system triggers successive dependent failures through positive feedback loops. This process typically continues until no further failures occur or the system stabilizes. A failure of a small fraction of nodes in one network can initiate an iterative cascade of failures in several interdependent networks [7]. This work constitutes the first attempt at developing network reconstruction methods for solving data and uncertainty challenges in complex infrastructure systems.

Network reconstruction from observations of the dynamics on the target network is a fundamental yet challenging inverse problem [12]–[14]. The research advancements centering around this important problem can be categorized into (i) capturing network structure using partial information [15], [16], and (ii) anticipating dynamic changes [17], [18]. While applications such as social, biological, and ecological networks have obvious dynamic patterns, the structure of critical infrastructure networks has been assumed to be static and deterministic unless the network is subject to uncertain disruptions. Network reconstruction has been applied to an individual critical infrastructure network scenario, such as road networks [19]; however, reconstructing interdependencies of multiple infrastructure sectors is still lacking. To achieve more representative network interdependencies, we aim to propose a data-driven network reconstruction approach to infer the structure of ICI networks from (observed) cascading failures.

1.1 Background and Related Work↩︎

Network reconstruction aims to infer the network topology from direct or indirect data that are inherently connected to the underlying network, such as observations of dynamic processes on the network [2], [3]. Instances include the diffusion of news among social media [20], the propagation of cascading failures in infrastructure networks [21], and the spreading of influence among political parties [22], among others. Although a single sequence of cascading failures may not provide sufficient information on the underlying network, combining sufficiently rich sequences of cascading failures can enable robust network reconstruction to obtain valuable insights into the functionality of the underlying network [23]. Data-driven network reconstruction approaches can be divided into: (i) graph embedding-based approaches [24]–[27], (ii) optimization-based approaches [28], [29], and (iii) the Bayesian approaches [12]–[14]. The first two approaches require information that is typically unavailable for critical infrastructure, whereas Bayesian approaches are commonly used to uncover the complete network structure from partial observations while accounting for uncertainty. Bayesian approaches for network reconstruction are either (i) parametric or (ii) nonparametric. Parametric approaches require the network type to be known, whereas many real-world systems are comprised of a diverse set of interconnected networks that do not fall into a well-defined network type. Nonparametric approaches employ a probabilistic approach via a sufficient number of samples of possible topologies to identify the most probable network. Prior studies equip Bayesian approaches with the epidemic spreading models and the information diffusion cascade models to infer structure from network dynamics [13], [14]. This study develops a nonparametric Bayesian approach to infer the discrete posterior distribution for the topology of the target infrastructure network from observed cascading failures.

1.2 Paper Contributions↩︎

In addition to the novel application of network reconstruction approaches in ICIs, this work advances the methodology in two ways. First, we devise an infrastructure-dependent proposal (IP) that incorporates topological and physical constraints of critical infrastructure and significantly reduces the sampling space in the application of the Markov Chain Monte Carlo (MCMC) algorithm for network reconstruction. Using the block structure of ICI networks, we leverage the hierarchical stochastic block model (HSBM) to calculate the prior probability of the potential graphs. Moreover, the convergence of the MCMC algorithm is ensured and numerically demonstrated. The second methodological contribution lies in the design of three computational optimization techniques to improve the efficiency of the MCMC algorithm in reconstructing large-scale networks where the number of possible network topologies increases exponentially over the number of nodes. The performance of these three techniques is numerically demonstrated via experiments on reconstructing synthetic ICI networks.

The remainder of this paper is structured as follows: Section 2 introduces the general Bayesian model for reconstructing the network topology based on cascading failures data. The MCMC algorithm with the infrastructure-dependent proposal (IP) to infer the Bayesian model is described in Section 3. In Section 4, we present three computational optimization techniques to improve the computational efficiency of the proposed network reconstruction approach. Extensive numerical experiments for validating the proposed approach are shown in Section 5, followed by concluding remarks in Section 6.

2 Bayesian Approach for Network Reconstruction↩︎

In the Bayesian approach for network reconstruction, we estimate the discrete posterior distribution of the network topology given the observations, i.e., a set of initial node failure scenarios and the following cascading failure sequences. Specifically, conditioned on the cascading failure data \(\boldsymbol{C}\), we estimate the discrete posterior distribution \(P(\boldsymbol{A}|\boldsymbol{C})\) for the adjacency matrix \(\boldsymbol{A}\) of the underlying network via Bayes’ rule [13], [14] as \[\begin{align} \label{eq:bayesianupdate} & P(\boldsymbol{A}|\boldsymbol{C}) = {P(\boldsymbol{C}|\boldsymbol{A})P(\boldsymbol{A})}/{P(\boldsymbol{C})} \propto P(\boldsymbol{C}|\boldsymbol{A})P(\boldsymbol{A}), \end{align}\tag{1}\] where \(P(\boldsymbol{C}|\boldsymbol{A})\) is the likelihood that the cascading failures \(\boldsymbol{C}\) occur in a network with topology \(\boldsymbol{A}\), \(P(\boldsymbol{A})\) encodes the prior information on the network topology, and \(P(\boldsymbol{C})\) is the normalization constant which represents the total evidence for the cascading failure data \(\boldsymbol{C}\) [13]. Details on calculating the graph prior \(P(\boldsymbol{A})\) and the likelihood \(P(\boldsymbol{C}|\boldsymbol{A})\) are presented in the following two sections.

2.1 Hierarchical Stochastic Block Model (HSBM)↩︎

Calculating the \(P(\boldsymbol{A})\) requires realistic possible samples of the target network to be reconstructed. Various types of networks having various prior probabilities \(P(\boldsymbol{A})\) can be generated by various graph models. The adopted graph model should represent realistic behaviors of the target network. Generative graph models for real networks, such as the classical random, scale-free, and small-world networks, cannot represent real-world ICIs as nodes in infrastructure networks are grouped by blocks and hierarchy [8]. Moreover, random graphs will generate unrealistic networks that do not comply with the physical constraints and properties of real infrastructure networks. As such, a suitable model for characterizing ICI networks would be the HSBM [30].

The HSBM is built upon the Stochastic Block Model (SBM), which divides different nodes into different blocks based on the node memberships. The HSBM adds hierarchical structures in each block according to node functionality levels. The probability of an edge between two nodes depends on the blocks to which the nodes belong and the hierarchical levels at which the nodes are positioned. Considering a multi-layer network with a set of blocks \(\mathcal{M}\) and the number of nodes in each block is \(n_{M}, M\in\mathcal{M}\), each single block \(M \in \mathcal{M}\) has a hierarchical structure \(\mathcal{L}\) where \(n_{M}\) nodes are further divided into \(|\mathcal{L}|\) levels. Considering the set of nodes and edges in the multi-layer network as \(\mathcal{V}^{\mathcal{M}} = \bigcup\limits_{M\in\mathcal{M}}{\mathcal{V}^M}\) and \(\mathcal{E}^{\mathcal{M}} = \bigcup\limits_{M\in\mathcal{M}}{\mathcal{E}^M}\), and denoting the block and the level labels of the node \(k\) as \(m_k\in\mathcal{M}\) and \(l_k \in \mathcal{L}\), the probability of an edge \((k,j)\) is an independent Bernoulli random variable \(p_{kj}\) conditioned on the block and the level labels \(m_k, m_j, l_k, l_j\) of nodes \(k\), \(j\), respectively. Then, the prior probability conditioned on the block and the layer assignment is calculated as \[\begin{align} & P(\boldsymbol{A}|\mathcal{M}, \mathcal{L}) = \prod\limits_{\forall k, j\in \mathcal{V}^{\mathcal{M}},k \ne j}{p_{kj}^{\boldsymbol{A}_{kj}}(1 - p_{kj})^{1 - \boldsymbol{A}_{kj}}}\\ & p_{kj} = g(m_{k}, m_{j}, l_{k}, l_{j}) \end{align}\]where \(g()\) is an HSBM function calculating the edge probability based on the block and the layer assignments of nodes. Conditioning the graph prior probability on the HSBM model, we introduce prior knowledge on the network topology \(P(\boldsymbol{A}|\mathcal{M}, \mathcal{L})\) from the function \(g()\). In ICIs, \(\mathcal{M}\) corresponds to individual infrastructure sectors (e.g., power grid and gas distribution network) and \(\mathcal{L}\) represents different types of facilities with different functionalities (e.g., power supply and distribution nodes).

2.2 Cascading Failure Model↩︎

After employing the HSBM to model the prior information of the network topology, we use the cascading failure model to calculate the likelihood in 1 . Denoting the adjacency matrix of the network to be reconstructed as \(\boldsymbol{A}^{\ast}\), and the prior knowledge about the network topology is updated with the likelihood \(P(\boldsymbol{C}|\boldsymbol{A})\) towards the target topology \(\boldsymbol{A}^{\ast}\) determining the diffusion of the cascading failure encoded in \(\boldsymbol{C}\). Disasters cause initial failures of infrastructure components, which trigger flow redistribution of resources and energy within and across the infrastructure networks [31], causing overload and damage to additional components, leading to cascading failures until the system reaches a new steady-state with no newly-failed components. We leverage the susceptible-infected (SI) epidemic model to simulate the cascading failure across infrastructure sectors, wherein the facility failure is modeled by the node infections and the flow redistribution is modeled by the probabilistic propagation of the infection in the network.

Due to the independence between any two distinct cascading failure scenarios, the likelihood of the cascading failure data \(\boldsymbol{C}\) encodes the failure propagation across all cascading failure scenarios \(i \in \{1,\dots,C\}\), over all time-steps for the entire disruption duration \(\{1,\dots,T(\boldsymbol{c}^i)\}\), and among all nodes in the multi-layer network \(j\in\mathcal{V}^{\mathcal{M}}\). Such a likelihood is given as \[\begin{align} \label{eq:likelihood} & P(\boldsymbol{C}|\boldsymbol{A}) = \prod\limits^{C}_{i=1}{\prod\limits_{t = 1}^{T(\boldsymbol{c}^i)}{\prod\limits_{j=1}^{|\mathcal{V}^{\mathcal{M}}|}{P(\boldsymbol{c}^i_{t +1,j}|\boldsymbol{c}^i_{t})}}} \end{align}\tag{2}\] where \(T(\boldsymbol{c}^i)\) is the number of time-steps that the cascading failure scenario \(\boldsymbol{c}^i\) lasts and \(P(\boldsymbol{c}^i_{t+1, j}|\boldsymbol{c}^i_t)\) is the probability of the status of operation of node \(j\) at time-step \(t + 1\), which is calculated by incorporating operable node status as \[\begin{align} & P(\boldsymbol{c}^i_{t+1,j}|\boldsymbol{c}^i_{t}) = \nonumber\\& {\Bigg( {1 - {\prod _{k \in {{\mathcal{V}}^{\mathcal{M}}}\backslash \{j\} }}\left[ {1 - {q_{kj}}{\boldsymbol{A}_{kj}}\boldsymbol{c}_{t,k}^i (1 - {\rm{ }}\boldsymbol{c}_{t - 1,k}^i)} \right]} \Bigg) ^{\boldsymbol{c}_{t + 1,j}^i(1 - \boldsymbol{c}_{t,j}^i)}}\nonumber\\ & \times {\Bigg( { {\prod _{k \in {{\mathcal{V}}^{\mathcal{M}}}\backslash \{j\} }}\left[ {1 - {q_{kj}}{\boldsymbol{A}_{kj}}\boldsymbol{c}_{t,k}^i(1 - {\rm{ }}\boldsymbol{c}_{t - 1,k}^i)} \right]} \Bigg) ^{(1-\boldsymbol{c}_{t + 1,j}^i)(1 - \boldsymbol{c}_{t,j}^i)}} \end{align}\]where \(q_{kj}\) denotes the probability of failure propagating from node \(k\) to node \(j\). The term \(\boldsymbol{c}^i_{t, k}(1 - \boldsymbol{c}^i_{t-1, k})\) ensures that the failure propagation is Markovian, i.e., the failure of nodes at the current time-step is only impacted by nodes that failed at the last time-step, matching the failure process of infrastructure networks, such as the progressive collapse of structures and gradual outages of the power stations [32].

3 MCMC for network reconstruction↩︎

By generating possible network topologies with the graph prior model from the HSBM and incorporating the cascading failure model, 1 becomes \[\begin{align} \label{eq:bayesianfull} & P(\boldsymbol{A}|\boldsymbol{C}, \mathcal{M}, \mathcal{L}) \propto \underbrace{P(\boldsymbol{C}|\boldsymbol{A}, \mathcal{M}, \mathcal{L})}_{\text{likelihood}}\underbrace{P(\boldsymbol{A} | \mathcal{M}, \mathcal{L}) P(\mathcal{M}, \mathcal{L})}_{\text{prior probability}}. \end{align}\tag{3}\] Since the cascading failure scenarios \(\boldsymbol{C}\) are independent of the blocks and the layers assignments, the likelihood \(P(\boldsymbol{C}|\boldsymbol{A}, \mathcal{M}, \mathcal{L})\) reduces to \(P(\boldsymbol{C}|\boldsymbol{A})\). Then, since \(P(\mathcal{M}, \mathcal{L})\) does not contain \(\boldsymbol{A}\), 3 reduces to \[\begin{align} \label{eq:finalform} & P(\boldsymbol{A}|\boldsymbol{C}, \mathcal{M}, \mathcal{L}) \propto P(\boldsymbol{C}|\boldsymbol{A}) P(\boldsymbol{A}|\mathcal{M}, \mathcal{L}). \end{align}\tag{4}\]

We clarify that while our approach utilizes a Hierarchical Stochastic Block Model (HSBM) and an SI likelihood—both of which involve parameter estimation—we characterize the method as nonparametric in the sense that it avoids assuming a rigid, predefined network topology, allowing the latent structure to be inferred directly from the data.

3.1 Metropolis-Hastings (M-H) Algorithm↩︎

Since 4 generally does not admit a closed-form solution, simulation techniques are usually leveraged to generate samples of the discrete posterior distributions of possible network topologies. MCMC methods such as the Metropolis-Hastings (M-H) algorithm and Gibbs sampling algorithm are commonly used to perform Bayesian inference via sampling [33], [34]. We choose the M-H algorithm over the Gibbs sampling algorithm as the latter requires an analytical solution to the conditional distributions of each model parameter. At each iteration, a new multi-layer network \(\mathcal{M}'\) with topology \(\boldsymbol{A}'\) is proposed based on the current network \(\mathcal{M}\) with topology \(\boldsymbol{A}\) via the proposal distribution \(Q(\boldsymbol{A}'|\boldsymbol{A})\). The algorithm probabilistically accepts the new proposal with the following acceptance ratio: \[\begin{align} \label{eq:mcmcupdate} \gamma = \min \bigg\{ 1, \frac{P(\boldsymbol{C}|\boldsymbol{A}')P(\boldsymbol{A}'|\mathcal{M},\mathcal{L})}{P(\boldsymbol{C}|\boldsymbol{A})P(\boldsymbol{A}|\mathcal{M},\mathcal{L})}\frac{Q(\boldsymbol{A}|\boldsymbol{A}')}{Q(\boldsymbol{A}'|\boldsymbol{A})} \bigg \}. \end{align}\tag{5}\]

The performance of the algorithm relies on the choice of the \(Q(\boldsymbol{A}'|\boldsymbol{A})\). The random-walk graph proposal is a common choice where a random node pair is selected either by adding new edges or removing the existing edges between them. However, this approach can generate invalid or unrealistic proposals for infrastructure networks. For instance, proposals of possible water-power networks, a likely edge is from pumping stations to storage tanks, representing water extraction from nearby rivers at pumping stations and its transfer to storage tanks to be stored for future use. Nevertheless, random graph generation may add edges going from storage tanks to pumping stations, which is unrealistic. We thus devise an infrastructure-dependent proposal (IP) method that imposes additional topology constraints on the proposed networks.

3.2 Infrastructure-dependent Proposal (IP)↩︎

The infrastructure-dependent proposal (IP) ensures that only networks with a topology conforming to ICIs are admissible. This is achieved by first analyzing the standard topology of ICIs, from which we abstract the topological constraints for designing the IP. Each network in ICIs has a hierarchical structure with three levels (i.e., \(|\mathcal{L}| = 3\)) corresponding to supply, transmission, and demand facilities. Given a set of infrastructure networks \(\mathcal{M}\) where supply, transmission, and demand nodes are denoted as \(\textit{s}\), \(\textit{t}\), and \(\textit{d}\), and a set of interdependent links \(\mathcal{I}\) across the networks (each interdependent link \(I \in \mathcal{I}\) is associated with an ordered pair of infrastructure networks), we can express the topological constraints as follows:

1. \(\forall M\in \mathcal{M}, \forall k\in \mathcal{V}_{\text{s}}^M, \exists j\in \mathcal{V}_{\text{d}}^M, k \mathop{\rightarrow}\limits^{\text{path}} j\).

2. \(\forall M\in \mathcal{M}, \forall k \in \mathcal{V}_{\text{d}}^M, \exists j\in \mathcal{V}_{\text{s}}^M, j\mathop{\rightarrow}\limits^{\text{path}} k\).

3. \(\forall M\in \mathcal{M}, \forall k\in \mathcal{V}_{\text{t}}^M, \exists j\in \mathcal{V}_{\text{d}}^M, k\mathop{\rightarrow}\limits^{\text{edge}} j\).

4. \(\forall M\in \mathcal{M}, \forall k\in \mathcal{V}_{\text{t}}^M, \exists j\in \mathcal{V}_{\text{s}}^M, j\mathop{\rightarrow}\limits^{\text{edge}} k\).

5. \(\forall I\in \mathcal{I}, \forall k \in \mathcal{V}_{\text{d}}^I, \exists j\in \mathcal{V}_{\text{s}}^I, k\mathop{\rightarrow}\limits^{\text{edge}} j\).

6. \(\forall I\in \mathcal{I}, \forall k \in \mathcal{V}_{\text{s}}^I, \exists j\in \mathcal{V}_{\text{d}}^I, j\mathop{\rightarrow}\limits^{\text{edge}} k\).

7. \(\forall M\in \mathcal{M}, \forall (k, j)\in \mathcal{E}^M, k\in \mathcal{V}_1, j\in\mathcal{V}_2, (\mathcal{V}_1, \mathcal{V}_2) \in \{(\mathcal{V}^M_{\text{s}}, \mathcal{V}^M_{\text{t}}), (\mathcal{V}^M_{\text{t}}, \mathcal{V}^M_{\text{d}}), (\mathcal{V}^M_{\text{s}}, \mathcal{V}^M_{\text{d}})\}\).

8. \(\forall I\in \mathcal{I}, \forall (k, j)\in \mathcal{E}^I, k\in \mathcal{V}^I_d, j\in\mathcal{V}^I_s\).

9. \(\forall M\in \mathcal{M}\), no cycles in \(M\).

10. \(\forall I\in \mathcal{I}\), no cycles in \(I\).

Customarily, resources are generated at the supply nodes and transported via the transmission nodes to the demand nodes, where the resources are further distributed to the residents. Thus, in every infrastructure network \(M\in\mathcal{M}\), every supply node \(k \in \mathcal{V}^M_{\text{s}}\) is connected by a path to at least one demand node \(j \in \mathcal{V}^M_{\text{d}}\) and vice versa, corresponding to constraints [constraint1]–[constraint2]. Every transmission node \(k \in \mathcal{V}^M_{\text{t}}\) is connected by an edge to at least one demand node \(j \in \mathcal{V}^M_{\text{d}}\), corresponding to the constraint [constraint3]. Also, every transmission node \(k \in \mathcal{V}^M_{\text{t}}\) is connected by an edge from at least one supply node \(j \in \mathcal{V}^M_{\text{d}}\), corresponding to the constraint [constraint4]. Similarly, for every interdependent link \(I \in \mathcal{I}\), we have the same constraints for connectivity from the demand nodes to the supply nodes, according to constraints [constraint5]–[constraint6]. The resources move from the supply nodes (high-level) to the demand nodes (low-level) with no resources flowing back (either directly or indirectly via the transmission nodes), which results in only forward edges from the high-level nodes to the low-level nodes according to constraint [constraint7], which automatically leads to no cycles in each infrastructure network as shown by constraint [constraint9]. For every interdependent link, the directions of the edges are from the demand nodes to the supply nodes according to constraint [constraint8], which automatically leads to no cycles in each interdependent link, as shown by constraint [constraint10]. The rationale behind the existence of the transmission nodes is driven by the physical nature of infrastructure networks. In most spatial network models, transmission nodes are employed to represent the aggregation points. While constraint [constraint7] allows for direct supply-to-demand edges, the intermediate transmission nodes are necessary for two main reasons: (i) dimensionality reduction (avoiding highly dense network connections via a limited set of transmission nodes), and (ii) economies of scale (enabling unit cost minimization via bulk handling at aggregated transmission nodes). From a sampling perspective, multi-hop paths (i.e., having sequentially connected transmission nodes) significantly increase the diameter of the graph. This often slows down the convergence of MCMC-based sampling schemes because the proposal must now account for a chain of dependencies rather than a simple bipartite or tripartite structure. Therefore, we have not included \((\mathcal{V}_t^M,\mathcal{V}_t^M)\) in constraint [constraint7].

The M-H algorithm with the IP and considerations of the \(10\) constraints is outlined in Alg. \(1\). First, a basic adjacency matrix \(\boldsymbol{A}^0\) is generated following lines \(1\)–\(6\) to initiate the algorithm. While we construct \(\boldsymbol{A}^0\) by connecting failed nodes in sequential time-steps, any method that generates graphs with \(P(\boldsymbol{C}|\boldsymbol{A}) >0\) can be used. In lines \(10\)–\(11\), we randomly select a node pair \((k,j)\) from the subset of the feasible set \(\mathcal{S}=\{(\mathcal{V}^{M}_{\text{s}}, \mathcal{V}^{M}_{\text{t}}), (\mathcal{V}^{M}_{\text{s}}, \mathcal{V}^{M}_{\text{d}}), (\mathcal{V}^{M}_{\text{t}}, \mathcal{V}^{M}_{\text{d}}), (\mathcal{V}^{I}_{\text{d}}, \mathcal{V}^{I}_{\text{s}})\}, \forall M\in\mathcal{M}, \forall I\in\mathcal{I}\}\), which guarantees the satisfaction of constraints [constraint7]–[constraint10]. We then add or remove the edge between that node pair \((k,j)\) and generate a new candidate network following lines \(12\)–\(16\). The candidate network is accepted as a sample of the discrete posterior distribution of \(\mathcal{A}\), if the addition or removal of this edge does not violate constraints [constraint1]–[constraint6]. Otherwise, we reject this candidate, propose another node pair, and check the feasibility of the proposed topology, as shown in lines \(17\)–\(25\).

The advantage of this IP over the traditional (original) one lies in the improvement in the sampling efficiency. By imposing the additional topological constraints, the invalid topology proposals are excluded, thereby increasing the accuracy of sampling. Also, such exclusion of the invalid topology proposals shrinks the space of the candidate topology (See Appendix 7), which significantly decreases the mixing time of the M-H algorithm. We also prove the convergence of the Markov chain formed by the IP (See Appendix 8).

Alg. 1: M-H algorithm with the IP

4 Sampling Efficiency Improvements↩︎

Although the IP significantly shrinks the space of the candidate network topology and thus improves the sampling efficiency, the computational expense can still be very high due to the computation of proposing and updating the topology of the network at each iteration. At each iteration, lines \(8\)–\(26\) in Alg. \(1\) take \(\mathcal{O}(1)\) operations to randomly choose node pairs from the feasible set and add or remove the associated edges. Evaluating the likelihood \(P(\boldsymbol{C}|\boldsymbol{A})\) in calculating the acceptance ratio \(\gamma\) takes \(\mathcal{O}(\sum_{i=1}^{C}{T(\boldsymbol{c}^i)|\mathcal{V}^\mathcal{M}|^2})\) time to double-loop over all nodes in \(\mathcal{V}^\mathcal{M}\) at time-step \(t\) to calculate the failure probability of every node at \(t+1\) and iterate over all time-steps in \(\{1,\dots,T(\boldsymbol{c}^i) - 1\}\) through all cascading failure scenarios \(\boldsymbol{c}^i, i\in\{1,\dots,C\}\). Validating whether the generated network \(\boldsymbol{A}\) violates constraints [constraint1]–[constraint6] requires \(\mathcal{O}(|\mathcal{V}^\mathcal{M}|(|\mathcal{V}^\mathcal{M}| + |\mathcal{E}^\mathcal{M}|) + |\mathcal{V}^\mathcal{I}|(|\mathcal{V}^\mathcal{I}| + |\mathcal{E}^\mathcal{I}|))\)—\(\mathcal{V}^\mathcal{I} := \bigcup_{I \in \mathcal{I}} \mathcal{V}^{I}\) and \(\mathcal{E}^{\mathcal{I}} = \bigcup_{I \in \mathcal{I}} \mathcal{E}^{I}\)—as we check the existence of paths by tree-search algorithms such as Depth First Search (DFS) or Breadth First Search (BFS) on every node in every infrastructure network and every interdependent link. Each iteration in Alg. \(1\) takes \[\begin{align} & {\mathcal{O}}(\underbrace {\sum\limits_{i = 1}^C {T({\boldsymbol{c}^i})|{{\mathcal{V}}^{\mathcal{M}}}} {|^2}}_{{\text{Likelihood calculation}}} + \underbrace {|\mathcal{V}^\mathcal{M}|(|\mathcal{V}^\mathcal{M}| + |\mathcal{E}^\mathcal{M}|) + |\mathcal{V}^\mathcal{I}|(|\mathcal{V}^\mathcal{I}| + |\mathcal{E}^\mathcal{I}|)}_{{\text{Graph validation}}}) \end{align}\]which is nonlinear due to \(|\mathcal{V}^\mathcal{M}|^2\) in the likelihood calculation and the \(|\mathcal{V}^{\mathcal{M}}|^2 + |\mathcal{V}^{\mathcal{I}}|^2\) in the graph validation. Since the interdependent infrastructure network is a sparse and tripartite graph (s-t-d), we devise the following three optimization techniques to speed up the computation:

4.1 Likelihood Calculation↩︎

The quadratic term \(|\mathcal{V}^\mathcal{M}|^2\) is due to evaluating the failure probability of each node \(j\) by considering failure propagation from all other nodes, \(1 - \prod\limits_{k\in \mathcal{V}^\mathcal{M}\backslash \{j\}}(1 - q_{kj}\boldsymbol{A}_{kj}\boldsymbol{c}^i_{t, k}(1 - \boldsymbol{c}^i_{t-1, k}))\). This computation can be simplified by only considering the failed neighborhoods of node \(j\), which is \(1 - {\prod _{k \in {{\mathcal{N}}^{\mathcal{M}}}(j) }}(1 - {q_{kj}}{\boldsymbol{c}}_{t,k}^i(1 - {\boldsymbol{c}}_{t - 1,k}^i))\). Instead of double-looping over all nodes, we traverse the edge list and integrate the failure propagation probability over edges with failed heading nodes and the same tailing nodes. Since the whole process requires iterating over the adjacency list (an array of linked lists describing the neighbors of a particular node), the time complexity of the likelihood calculation reduces to \(\mathcal{O}(\sum_{i=1}^{C}{T(\boldsymbol{c}^i)(|\mathcal{V}^\mathcal{M}| + |\mathcal{E}^\mathcal{M}|)})\).

4.2 Graph Validation↩︎

The second-order term \(|\mathcal{V}^\mathcal{M}|(|\mathcal{V}^\mathcal{M}| + |\mathcal{E}^\mathcal{M}|) + |\mathcal{V}^\mathcal{I}|(|\mathcal{V}^\mathcal{I}| + |\mathcal{E}^\mathcal{I}|)\) is due to checking the corresponding paths on every node in the network. However, considering the tripartite structure of ICI networks and since every iteration of the M-H algorithm only adds or removes one edge, Theorem 1 guarantees that constraints [constraint1]–[constraint6] are validated in linear time.

The proof of Theorem 1 is presented in Appendix 9. Theorem 1 indicates that we can validate the proposed topology by only checking the out-degree of the heading node \(k\) and the in-degree of the tailing node \(j\), which only takes linear time \({\mathcal{O}(|{\mathcal{V}^{{\mathcal{M}}}}| + |{\mathcal{E}^{{\mathcal{M}}}}| + |{\mathcal{V}^{{\mathcal{I}}}}| + |{\mathcal{E}^{{\mathcal{I}}}}|)}\).

4.3 Tie-No-Tie Sampler↩︎

Another technique to improve the sampling efficiency is to replace the random sampler at line \(10\) in Alg. \(1\) with the ‘tie-no-tie’ (TNT) sampler [14]. The random sampler selects node pairs with equal probability, resulting in a higher frequency of proposing node pairs with non-edges (edges not present in the network) and adding the corresponding edges, then proposing node pairs with edges and removing the corresponding edges in sparse graphs. However, the newly added edges are likely to be rejected due to the violation of topological constraints [constraint1]–[constraint10] or the generation of a network topology discouraged by the cascading failure data. Conversely, the TNT sampler first selects the set of edges or the set of non-edges with equal probability, and then changes a random node pair in that set, which increases the frequency of removing edges in the proposal. Since the removed edges are more likely to be accepted due to graph sparsity, the sampling efficiency is improved.

5 Numerical Experiments↩︎

Here, we apply the proposed Bayesian approach to one synthetic ICI. We design five methods from different combinations of the IP and different computational optimization techniques (Table 1). By comparing the performance of these methods on reconstructing the synthetic system of ICI networks, we can identify the best method.

5.1 Network Reconstruction of the Synthetic Networks↩︎

To assess the effectiveness of the proposed network reconstruction Bayesian approach, we primarily consider the class of synthetic networks. The synthetic network contains \(3\) blocks corresponding to water (W), power (P), and gas (G) networks, each of which has \(3\) hierarchical levels corresponding to supply s, transmission t, and demand d facilities (Fig. 1 (a)). All \(3\) networks have \(2\) supply nodes, \(3\) transmission nodes, and \(5\) demand nodes. The interdependency across networks is simulated using a previously developed approach for generating Synthetic ICI Networks (SICINs) [6], which requires as input the degree distribution extracted from real-world infrastructures of the same type. This prior work [6] introduces a method for generating SICINs to overcome the research bottleneck caused by a lack of real-world infrastructure data. The SICIN framework utilizes simulation and nonlinear optimization. It employs a modified simulated annealing (MSA) algorithm to determine optimal asset placement within individual networks. Moreover, it uses a novel pseudo-tripartite graph algorithm to generate the complex links that connect different infrastructure systems. SICIN has been previously validated through a comparison with real-world infrastructure networks that led to the smallest difference in the adjacency matrices (DA). Particularly, for the baseline case of the Water-Power-Gas (WPG) ICI network in Shelby County, TN, USA, almost all properties (i.e., topological characteristics) of SICIN-based synthetic networks are the closest to the real ICI network of Shelby County, affirming the trustworthiness of the SICIN framework developed by [6].

We consider \(4\) types of interdependencies: \(1\)) pumping stations require electricity from \(12\)kV substations for pumping water from nearby rivers \({P \rightarrow W}\), \(2\)) power gate stations rely on water from water delivery stations for cooling purposes \({W \rightarrow P}\), \(3\)) gas gate stations require electricity from \(12\)kV substations for extracting natural gas from underground \({P \rightarrow G}\), and \(4\)) power gate stations depend on natural gas from gas delivery stations for generating electricity \({G \rightarrow P}\). We denote the water-power-gas (WPG) SICIN to be reconstructed as \(\boldsymbol{A}^{\ast}\) hereafter. In summary, for the WPG SICIN, we have \(\mathcal{M} = \{W,P,G\}\) and \(\mathcal{I} = \{P \rightarrow W, W \rightarrow P, P \rightarrow G, G \rightarrow P\}\).

The cascading failure data \(\boldsymbol{C}\) is simulated using the SI epidemic model, the ratio of the initial failed nodes \(f\) is set to \(0.2\), and the failure propagation probability \(q\) is set to \(0.4\) to get sufficient cascading failure data, aligned with prior work in the literature [14]. We investigate the performance of the proposed approach considering different amounts of cascading failure data \(\boldsymbol{C}\). We design \(3\) experiments with at least \(5\) time-steps each: (i) \(E_5^5\): \(5\) cascading failure scenarios, (ii) \(E_5^{15}\): \(15\) cascading failure scenarios, and (iii) \(E_5^{40}\): \(40\) cascading failure scenarios. All three experiments are performed with and without the IPs.

We evaluate the accuracy of network reconstruction by considering it as a classification problem. Suppose that the set of networks in the discrete posterior distribution is \(\mathcal{A}\), the edge probability \(\hat{p}_{kj}\) for each node pair \((k,j)\) is defined as the ratio of the number of graphs with the edge between \(k\) and \(j\) in the posterior and the total number of graphs in the posterior as follows: \[\begin{align} \label{eq:margin} & {\hat{p}_{kj}} = {|\{{\mathbf{A}\in \mathcal{A}}: \mathbf{A}_{kj} = 1\}|}/{|\mathcal{A}|}. \end{align}\tag{6}\] The edge probability for each node pair in the network is displayed on the heatmap of the adjacency matrix in Fig. 2. In Figs. 2 (a)–2 (c), some false edges that are not present in the real networks are proposed (i.e., false positives) and shown in blank areas. In contrast, when using the IP (shown in Figs. 2 (d)–2 (f)), the heatmap of the adjacency matrix conforms to the block structure in the original network (shown in Fig. 1 (b)). Such a difference in reconstruction accuracy highlights the advantage of topological constraints [constraint1]–[constraint10] incorporation via the IPs. The heatmap of the adjacency matrix is darker in \(E_5^{40}\) than in the other two experiments, indicating higher edge probability. Since \(E_5^{40}\) uses more cascading failure data to update the network topology, we have higher confidence in the predicted edges.

5.2 Performance of the Network Reconstruction Approach↩︎

5.2.1 Accuracy of the Network Reconstruction↩︎

We now examine the accuracy of the proposed network reconstruction approach. We set a probability threshold \(\rho\) and classify node pairs into two categories: node pairs connected by an edge if \(\hat{p}_{kj} \ge \rho\) and node pairs that are not connected by an edge if \(\hat{p}_{kj} < \rho\). We further count the number of edges for each node pair that are correctly classified and calculate the F1-score, which is the harmonic mean of precision and recall. The precision-recall curves are depicted in Fig. 3. Precisely, for \(10,001\) equidistant values of the probability threshold \(\rho\) ranging from \(0\) to \(1\), we calculate the F1-scores and identify the maximum F1-score among those \(10,001\) values. For the three experiments \(E_5^{5}\), \(E_5^{15}\), and \(E_5^{40}\), the maximum F1-scores are \((0.383,0.565,0.710)\) for experiments without the IP (Configuration: TNT + Edge List), and \((0.750,0.865,0.964)\) for experiments with the IP (Configuration: IP + TNT + Edge List), respectively. The accuracy of the reconstructed network improves with additional cascading failure data. Also, the IP significantly improves the performance of reconstructing networks using the same amount of cascading failure data.

To further shed light on the network reconstruction accuracy of the IP-based method, we conduct comparisons with a graph-embedding-based network reconstruction baseline [24], [25]. Prior work proposes the auto-encoder framework, namely variational graph auto-encoder (VGAE), consisting of a graph convolutional network (GCN) encoder and a simple inner product decoder used to reconstruct the adjacency matrix [25]. Precisely, using node features (e.g., node types or historical failure patterns), VGAE learns a latent representation for each node and predicts edges based on the similarity between these representations. Feeding the same cascading failure data \(\boldsymbol{C}\) to the VGAE, we obtain the following maximum F1-scores for the three experiments \(E_5^5\), \(E_5^{15}\), and \(E_5^{40}\): \((0.711,0.714,0.712)\). Comparing these values for the VGAE with those of the non-IP-based \((0.383,0.565,0.710)\) and the IP-based \((0.750,0.865,0.964)\) methods, we realize that the accuracy-based superiority is as follows: IP-based \(>\) VGAE \(>\) non-IP-based.

Among the non-IP-based, the IP-based, and the VGAE, Fig. 4 corroborates that the IP-based method is the only threshold \(\rho\)-based reconstructed network satisfying the topological constraints [constraint1]–[constraint10]. As mentioned earlier in Section 1, the graph-embedding-based approaches like the VGAE [25] require a massive amount of information to attain an acceptable reconstruction performance, while the Bayesian approaches like our IP-based method achieve an acceptable reconstruction performance given a relatively small amount of cascading failure data, which is driven by the specific nature of limited data availability for the ICI networks.

To further verify the applicability of the IP-based method to larger ICI networks, we generate a \(120 \times 120\) random WPG ICI network by setting edge probability \(p\) equal to \(0.5\) for feasible edges, ensuring the satisfaction of topological constraints [constraint1]–[constraint10]. All \(3\) networks have \(8\) supply nodes, \(12\) transmission nodes, and \(20\) demand nodes. Similarly, the cascading failure data \(\boldsymbol{C}\) is simulated using the SI epidemic model, the ratio of the initial failed nodes \(f\) is set to \(0.2\), and the failure propagation probability \(q\) is set to \(0.4\) to get sufficient cascading failure data. It is worth noting that generating the cascading failure data \(\boldsymbol{C}\) for larger ICI networks takes a longer time. Feeding the generated cascading failure data \(\boldsymbol{C}\) to Alg. \(1\) for an experiment with \(40\) cascading failure scenarios and with \(4\) time-steps, namely \(E_{4}^{40}\), we obtain the network reconstruction results in Fig. 5. Accordingly, Fig. 6 depicts the precision-recall curve associated with the network reconstruction in Fig. 5.

Repeating the process for \(10\) replicates, the average maximum F1-score and the average running time (s) are \(0.664\) and \(223.004\), respectively. The corresponding \(95\%\) confidence intervals (CIs) are \([0.660, 0.668]\) and \([211.969, 234.040]\), respectively, indicating the high confidence for the calculated average values. For relatively dense (due to \(p = 0.5\)) and large-scale ICI networks, increasing the number of cascading failure scenarios and the failure length at the expense of higher computational cost could lead to an improvement in the achieved fair/moderate network reconstruction accuracy (i.e., F1-score = \(0.664\)).

5.2.2 Effect of the Ratio of the Initial Failed Nodes↩︎

To investigate the effect of the ratio of the initial failed nodes \(f\) on the network reconstruction accuracy (quantified by the maximum F1-score), we repeat the simulations with distinctly generated cascading failure data \(\boldsymbol{C}\) by varying \(f\) from \(0.1\) to \(0.9\) with a step size of \(0.1\). The observations have been reflected in Fig. 7 for the three experiments \(E_5^{5}\), \(E_5^{15}\), and \(E_5^{40}\) without the IP (Configuration: TNT + Edge List) and with the IP (Configuration: IP + TNT + Edge List) for the WPG SICIN. As Fig. 7 depicts, the moderate-value choices of \(f \in \{0.2,0.3,0.4\}\) seem reasonable as such choices lead to better maximum F1-scores. It is particularly consistent with the choice of \(f = 0.2\) in prior work in the literature [14]. The observations in Fig. 7 further corroborate the fact that moderate values of \(f\) ensure the high inferential precision by triggering a controlled, multi-branch cascade that propagates organically through the network’s true paths.

On one hand, the low values of \(f\) usually lead to high ambiguity due to the lack of data to differentiate between many candidate network topologies and failing to reach deep, isolated, or sparsely connected sections of the network. On the other hand, the high values of \(f\) usually lead to significant loss of information due to overloading and suffering an immediate or near-instantaneous global collapse, weakening the effectiveness of the Bayesian inference. Moreover, comparing Fig. 7 (b) with Fig. 7 (a), we realize that the magnitudes of relative sensitivity (i.e., elasticity) values of the reconstructions with the IP, \((0.036,0.130,0.098)\), are relatively less than those of the counterparts without the IP, \((0.349,0.304,0.123)\), respectively. In other words, the IP-based network reconstruction approach outperforms the non-IP-based counterpart by achieving a less sensitive network reconstruction (in terms of accuracy) against the ratio of the initial failed nodes.

5.2.3 Effect of the Noise Level in the Cascading Failure Data↩︎

To explore the effect of the noise level \(\eta\) in the cascading failure data \(\boldsymbol{C}\) on the network reconstruction accuracy (quantified by the maximum F1-score), for a fixed value of \(f = 0.2\), we repeat the simulations with distinctly generated noisy cascading failure data \(\boldsymbol{C}\) with different equidistant noise levels \(\eta\)’s ranging from \(0\) (i.e., the noiseless baseline) to \(0.18\) with a step size of \(0.01\). By trial and error (e.g., a bisection-based search), we find that \(0.18\) is the maximum tolerable noise level \(\eta\) for which the network can be reconstructed. We have incorporated two classes of noises: (i) missing data (i.e., dropouts) simulating a sensor failing to record a node’s collapse, and (ii) timing jitter (i.e., permutations) simulating logging delays where node \(k\) is recorded as failing after node \(j\), even if \(k\) actually failed first. Through the noise level \(\eta\), we will parametrically simulate observational noise (missing nodes or reporting delays) by altering the failure tracks immediately after they are generated but before they are used to initialize the prior. The observations have been visualized in Fig. 8 for the three experiments \(E_5^{5}\), \(E_5^{15}\), and \(E_5^{40}\) without the IP (Configuration: TNT + Edge List) and with the IP (Configuration: IP + TNT + Edge List) for the WPG SICIN.

As Fig. 8 demonstrates, an overall descending trend (i.e., trade-off) exists between the noise level \(\eta\) and the maximum F1-score. Quantitatively, the Kendall tau values for the three experiments \(E_5^{5}\), \(E_5^{15}\), and \(E_5^{40}\) are: \((-0.194,-0.193,-0.453)\) for the case without the IP in Fig. 8 (a) and \((-0.342,-0.152,-0.324)\) for the case with the IP in Fig. 8 (b), respectively. A negative value of Kendall tau indicates the existence of an overall descending trend. Such a trade-off between the noise level \(\eta\) and the maximum F1-score is intuitively aligned with our expectation regarding an adverse effect of the noise on the network reconstruction accuracy. Furthermore, comparing Fig. 8 (b) with Fig. 8 (a), we realize that the magnitudes of relative sensitivity (i.e., elasticity) values of the reconstructions with the IP, \((0.017,0.014,0.014)\), are less than those of the counterparts without the IP, \((0.086,0.033,0.084)\), respectively. In other words, the IP-based network reconstruction approach outperforms the non-IP-based counterpart by achieving a less sensitive network reconstruction (in terms of accuracy) against the noise level.

5.3 Computational Performance of the Three Incorporated Optimization Techniques↩︎

To demonstrate the computational performance of the three incorporated optimization techniques, we devise five methods by combining different techniques and the IP. The details of the five methods are outlined in Table 1, where the \(✔\) denotes that the corresponding method uses the corresponding technique and Validation refers to validating the proposed graph via Theorem 1. We implement each method \(10\) times (i.e., \(10\) replicates) to reconstruct the synthetic networks and report their average (reconstruction) running time and their average maximum F1-score in Fig. 9. Given that the \(95\%\) confidence interval (CI)-based error bars in Fig. 9 are small, we are confident about the average values. In each run, we propose \(3,000\) samples and use the first \(2,000\) samples as a warm-up. The numerical simulations in Fig. 9 have been performed on a Windows 10 laptop with a 1.9 GHz Intel Core i7-8650 U CPU and 16 GB RAM.

In Fig. 9 (a), M-1 and M-2 have the shortest running time and the difference in running time between M-1, M-2, and the rest becomes larger as the amount of data increases from \(E_5^{5}\) to \(E_5^{40}\), indicating that traversing edge lists rather than adjacency matrix to calculate the likelihood significantly reduces the computational load and this reduction is significant when more failure data is used. The comparable running time between M-1 and M-2 indicates that applying graph validation according to Theorem 1 brings almost no improvement in computational efficiency. This is because proposing networks by removing edges in the TNT sampler can easily generate disconnected node pairs, and once finding such node pairs, constraints [constraint1]–[constraint6] are broken, leading to saving the time that would be otherwise needed to check other node pairs. Comparing M-3 and M-4, replacing the random sampler with the TNT sampler reduces the computational load because the TNT sampler proposes more valid topological variation by removing edges rather than adding invalid edges that would be later rejected. The running time of M-4 is greater than that of M-5 in \(E_5^{5}\) while it is vice versa in \(E_5^{15}, E_5^{40}\). Given the smaller amount of cascading failure data in \(E_5^{5}\), the proposed networks can easily be accepted without applying the IP to check any topological constraint in advance. Thus, the sampling process is completed in a shorter time in M-5, where no extra operations on checking the network topology are required. However, for larger cascading failure data in \(E_5^{15}\) and \(E_5^{40}\), and without verifying the network topology for topological constraints [constraint1]–[constraint10] in advance, the topology candidates proposed by the random sampler are more likely to be rejected as they are not supported by the cascading failure data (low likelihood). Then, the time spent on proposing such invalid samples is compensated by the time for checking the topological constraints, leading to shorter times in M-4 than M-5 for \(E_5^{15}\) and \(E_5^{40}\).

In Fig. 9 (b), the first four methods lead to roughly the same and better average maximum F1-score than M-5, which indicates the power of the IPs’ incorporation. Nevertheless, M-4, which applies the random sampler, is slightly better than the previous three methods that use the TNT sampler, as some edges that are present in the target network have already been removed by the TNT sampler in M-1, M-2, and M-3, and the random sampler has a higher chance of adding edges than the TNT sampler, thus the former is more likely to recover those edges that should exist in the target network but have already been removed. The average running times and the average maximum F1-scores of different experiments are provided in Table 2. Considering the trade-off between the running time and the reconstruction accuracy, M-1 is identified as the best method.

Table 2 indicates that M-4 achieves the best average maximum F1-score, which stems from the choice of the random sampler. This motivates us to further investigate the effect of sampler choice on the average running time and the average maximum F1-score. To that end, we conduct similar simulations for three methods (M-1, M-1-r, and M-2-r) along with \(2\) experiments (\(E_5^5\) and \(E_5^{15}\)) where the modified methods M-1-r and M-2-r are simply obtained from the methods M-1 and M-2 by only replacing the TNT sampler with the random sampler, respectively (i.e., keeping the rest of the techniques unchanged in the configuration). The average running times and the average maximum F1-scores of different experiments are provided in Table 3. Considering the trade-off between the running time and the reconstruction accuracy, M-2-r (i.e., Configuration: IP + Random + Edge List) and M-1 (i.e., Configuration: IP + TNT + Edge List + Validation) are suggested as the best method among the three methods for the experiments \(E_5^5\) and \(E_5^{15}\), respectively. It is worth noting that in the case of \(5\) cascading failure scenarios, the best method is no longer identified as M-1. However, increasing the amount of data to the \(15\) cascading failure scenarios, the best method is still identified as M-1. Such an observation highlights the positive impact of collecting more cascading failure scenarios on the reconstruction accuracy. Furthermore, it is remarkable that depending on the considered weighted sum of effects of the average (reconstruction) running time and the average maximum F1-score, e.g., \(w/T_{\text{avg}} + (1-w) F_1\) (weight: \(w \in [0,1]\)) and the choice of weight \(w\), the best method could be identified differently. Therefore, the human operator can play a significant role in terms of striking a balance between the average (reconstruction) running time and the average maximum F1-score via properly selecting the weighted sum coefficient \(w\) based on the physics-informed priorities.

6 Conclusion↩︎

This paper presents a Bayesian approach for topology reconstruction of the interdependent critical infrastructure (ICI) networks based on observed cascading failure data, which is the first application of topology reconstruction methods for ICI networks. The hierarchical stochastic block model (HSBM) is used to capture the clustering and hierarchical structure of ICI networks. The computational complexity resulting from the exponential growth of the potential topology is addressed with a new infrastructure-dependent proposal (IP), taking into account the topological constraints of ICI networks. We show the superiority of the IP-based method over a graph-embedding-based approach, namely variational graph auto-encoder (VGAE), in terms of the accuracy of the network reconstruction. Furthermore, we propose three optimization techniques to significantly improve the sampling efficiency. Results of a synthetic network case study demonstrate the effectiveness and transferability of the proposed approach, making it an appropriate choice for the topology reconstruction of the real-world ICI networks, wherein data on the topology is often unavailable either due to security concerns or decentralized operations of different infrastructure sectors. The more accurately we reconstruct the topology of real-world ICI networks, the greater the potential for improved risk mitigation and economic benefits, as it provides policymakers and practitioners with a clearer basis for informed decision-making. We further verify the applicability of the IP-based method to larger ICI networks using a \(120 \times 120\) random WPG ICI network where the edge probability \(p\) is considered to be equal to \(0.5\) for feasible edges, ensuring the satisfaction of topological constraints [constraint1]–[constraint10] and highlighting the fact that for relatively dense (due to \(p = 0.5\)) and large-scale ICI networks, increasing the number of cascading failure scenarios and the failure length at the expense of higher computational cost could lead to an improvement in the achieved fair/moderate network reconstruction accuracy (i.e., F1-score = \(0.664\)). Finally, we conduct sensitivity analyses of the network reconstruction accuracy with respect to the two factors: (i) the ratio of the initial failed nodes, and (ii) the noise level in the cascading failure data. The former sensitivity analysis suggests that the moderate-value choices of the initial failed nodes lead to better network reconstruction accuracy. The latter suggests that the noise in the cascading failure data has an adverse effect on the network reconstruction accuracy. Furthermore, both sensitivity analyses corroborate that the IP-based network reconstruction approach outperforms the non-IP-based counterpart by achieving a less sensitive network reconstruction (in terms of accuracy) against the two factors.

Future work can be carried out to thoroughly investigate the influence of cascading failure data, network structure, and network size on the reconstruction accuracy. While the proposed approach uses the susceptible–infected (SI) epidemic model to simulate cascading failures in infrastructures, future work can explore more advanced types of failure simulation models (and their stochastic variants) such as susceptible–infected–recovered (SIR). As another future direction, one can think of proposing a generalized probabilistic approach to evaluate infrastructure interdependencies in the case of uncertain and time-varying ICI networks.

acknowledgment↩︎

This material is based upon work supported by the National Science Foundation under Grant 1944559.

Data Availability↩︎

The data and software that support the findings of this article are openly available in Github at https://github.com/MirSaleh/PRE-Paper-Source-Code.

7 Space of Candidate Networks↩︎

Since the effect of reduction on the sampling space of network topology by considering additional constraints is universal among different networks, we only select one network \(I\in\mathcal{I}\) and calculate the total number of potential topologies with and without topological constraints [constraint1]–[constraint10] to demonstrate the significant decrease in the sampling space.

Before imposing additional topological constraints, the number of potential edges in the candidate network \(I\) is upper bounded by \(\binom{|{{\mathcal{V}}^I}|}{2}\). Since every edge in the above \(\binom{|{{\mathcal{V}}^I}|}{2}\) edges has two states (existent or non-existent), the total number of the candidate topology is \(2^{\binom{|{{\mathcal{V}}^I}|}{2}}\) whose complexity is an exponential of a quadratic term in terms of \(|{{\mathcal{V}}^I}|\).

After imposing topological constraints [constraint1]–[constraint10], the total number of potential topologies for the network \(I\) is equal to the total number of bipartite graphs consisting of \(|{{\mathcal{V}}^I}| = |\mathcal{V}_\text{s}^I| + |\mathcal{V}_\text{d}^I|\) non-isolated vertices (with \(|\mathcal{V}_\text{s}^I|\) vertices and \(|\mathcal{V}_\text{d}^I|\) vertices in the supply part and the demand part, respectively) that can be computed as follows: \[\begin{align} & \sum_{k=0}^{|\mathcal{V}_\text{d}^I|}~(-1)^k \binom{|\mathcal{V}_\text{d}^I|}{k}(2^{|\mathcal{V}_\text{d}^I|-k}-1)^{|\mathcal{V}_\text{s}^I|}. \label{eq:noisolated} \end{align}\tag{7}\]

It is worth noting that following a similar path (exchanging the supply nodes and the demand nodes), an alternative formula for 7 can be obtained as follows: \[\begin{align} & \sum_{k=0}^{|\mathcal{V}_\text{s}^I|}~(-1)^k \binom{|\mathcal{V}_\text{s}^I|}{k}(2^{|\mathcal{V}_\text{s}^I|-k}-1)^{|\mathcal{V}_\text{d}^I|}. \end{align}\]

Considering a network with \(2\) supply nodes on the top level and \(3\) demand nodes on the bottom level, the total number of its potential topologies is \(25\) according to 7 , which is a significant decrease from the original \(2^{\binom{5}{2}} = 1024\) different topology candidates.

8 Convergence Proof of Alg. \(1\)↩︎

The proof of Alg. \(1\) can be detailed as follows:

Proof. The limit theorem of Markov Chains states that if the Markov chain is Irreducible and Aperiodic, the chain will converge to a unique stationary distribution. Therefore, to prove the convergence of the Markov Chain in Alg. \(1\), we need to prove the above two properties for this Markov Chain. Let \(\boldsymbol{G}_0\), \(\boldsymbol{G}_1\), …, be the Markov chain \(\mathcal{G}\) of Alg. \(1\) where \(\boldsymbol{G}_i \in \mathcal{G}\) is the \(i\)-th state of the ICI networks and the movement from state \(\boldsymbol{G}_i\) to state \(\boldsymbol{G}_{i+1}\) is strictly governed by the infrastructure-dependent proposal as shown in Fig. 10.

Irreducible. A Markov chain is irreducible if every state can be reached from every other state, which can mathematically be formalized as \(\forall \boldsymbol{G}_i, \boldsymbol{G}_j \in \mathcal{G}, P(\boldsymbol{G}_j|\boldsymbol{G}_i) > 0\). To prove the irreducibility of \(\mathcal{G}\), we need to show that we can go from any state \(\boldsymbol{G}_i \in \mathcal{G}\) following the infrastructure-dependent proposal to any other state \(\boldsymbol{G}_j \in \mathcal{G}\) in finite steps. Denoting the edge sets of these two states as \(\mathcal{E}^{\boldsymbol{G}_i}, \mathcal{E}^{\boldsymbol{G}_j}\), their shared edge set as \(\mathcal{E}^{\boldsymbol{G}_i}\cap\mathcal{E}^{\boldsymbol{G}_j}\), and their difference as \(\mathcal{E}^{\boldsymbol{G}_i-\boldsymbol{G}_j} = \mathcal{E}^{\boldsymbol{G}_i} - \mathcal{E}^{\boldsymbol{G}_j}, \mathcal{E}^{\boldsymbol{G}_j-\boldsymbol{G}_i} = \mathcal{E}^{\mathbf{G}_j} - \mathcal{E}^{\boldsymbol{G}_i}\). We first construct a medium graph \(\boldsymbol{G}_{k}\) by adding all edges in \(\mathcal{E}^{\boldsymbol{G}_j - \boldsymbol{G}_i}\) consecutively to the graph \(\boldsymbol{G}_{i}\). Since \(\mathcal{E}^{\boldsymbol{G}_j - \boldsymbol{G}_i}\subset \mathcal{E}^{\boldsymbol{G}_j}\) and \(\mathcal{E}^{\boldsymbol{G}_j}\) satisfies constraints [constraint7]–[constraint8] , all edges in \(\mathcal{E}^{\boldsymbol{G}_j - \boldsymbol{G}_i}\) point from high-level nodes to low-level nodes in the feasible set \(\mathcal{S}\) and thus adding these edges during graph construction causes no cycles in \(\boldsymbol{G}_k\), i.e., no violation against constraints [constraint7]–[constraint10] . Besides, adding edges trivially satisfies constraints [constraint1]–[constraint6] . Therefore, the graph construction process from \(\boldsymbol{G}_i\) to \(\boldsymbol{G}_k\) is strictly governed by the infrastructure-dependent proposal. Then, we remove all edges in \(\mathcal{E}^{\boldsymbol{G}_i-\boldsymbol{G}_j}\) consecutively from the graph \(\boldsymbol{G}_{k}\) to obtain graph \(\boldsymbol{G}_j\). Since \(\mathcal{E}^{\boldsymbol{G}_j}\subset\mathcal{E}^{\boldsymbol{G}_k}\) and \(\mathcal{V}^{\boldsymbol{G}_j} = \mathcal{V}^{\boldsymbol{G}_k}\), the medium graph \(\boldsymbol{G}_k\) is formed by adding extra edges \(\mathcal{E}^{\boldsymbol{G}_i-\boldsymbol{G}_j}\) to the target graph \(\boldsymbol{G}_j\), which is also a spanning tree of \(\boldsymbol{G}_k\). Because \(\mathcal{E}^{\boldsymbol{G}_i-\boldsymbol{G}_j}\cap\mathcal{E}^{\boldsymbol{G}_j} = \varnothing\), removing edges in \(\mathcal{E}^{\boldsymbol{G}_i-\boldsymbol{G}_j}\) will not delete any edge in the spanning three \(\boldsymbol{G}_j\) of \(\boldsymbol{G}_k\) and this well-preserved spanning three \(\boldsymbol{G}_j\) guarantees that the medium graph \(\boldsymbol{G}_{k}\) always satisfies topological constraints [constraint1]–[constraint10] during edge removal. Therefore, the edge removal process from \(\boldsymbol{G}_k\) to \(\boldsymbol{G}_j\) is also strictly governed by the infrastructure-dependent proposal. Concatenating the process of edge addition and edge removal gives a way to go from \(\boldsymbol{G}_i\) to \(\boldsymbol{G}_j\) under the supervision of the infrastructure-dependent proposal as Fig. 11, thus the MC \(\mathcal{G}\) is irreducible.

Aperiodic. A Markov chain is aperiodic if every state is aperiodic. One easy way to determine whether a state \(i\) is aperiodic is to check whether state \(i\) has a self-loop, i.e., whether MC remains at the same state after moving one step forward from that state. Since every proposal in the M-H algorithm is rejected probabilistically, every state \(\boldsymbol{G}_i \in \mathcal{G}\) has a self-loop and therefore \(\mathcal{G}\) is aperiodic.

Combining the irreducible and aperiodic properties, we conclude that \(\mathcal{G}\) will converge to a unique stationary distribution. ◻

9 Proof of Theorem 1↩︎

First, we propose Lemma 1 that will be used to prove Theorem 1 in Section 4.

References↩︎