Competitive percolation strategies for network recovery

We are interested in the problem of restoring the operation of a critical infrastructure system after sustaining large-scale damage. The infrastructure network is represented by a graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$, where $\mathcal{N}$ is the set of $N$ nodes corresponding to substations and $\mathcal{E}$ is a set of $E$ edges corresponding to transmission lines, e.g., power lines, water or gas pipes. We introduce the parameter $d_{i}$ representing the commodity demand of node $i$: if $d_{i}<0$, node $i$ is a net consumer; if $d_{i}>0$, node $i$ is a net supplier. We normalize $d_{i}$ such that the total consumption (or production) sums to unity, i.e., $\frac{1}{2}\sum_{i=1}^{N}\left|d_{i}\right|=1$. Following a catastrophic event, a subset of the network $\mathcal{G}^{\prime}=(\mathcal{N}^{\prime},\mathcal{E}^{\prime})$ becomes damaged. We study a discrete time reconstruction process: in each timestep we fix one damaged component, and the process ends once the entire network is functional. Our goal is to identify a sequence in which to repair the elements such that the total cost of recovery is minimized. In this manuscript, we focus on the fundamental case where all links are damaged but nodes remain functional, i.e., $\mathcal{E}^{\prime}=\mathcal{E}$ and $\mathcal{N}^{\prime}=\emptyset$.

Optimization frameworks are often used in order to explore the space of possible repair sequences and identify best solutions. We implement the time-dependent network design problem (td-NDP) algorithm [22, 23, 30], which is a well known example of an optimization algorithm for network recovery developed by the civil engineering community. Out of the recovery processes we examine in this paper, td-NDP is the most realistic, and therefore the most computationally complex. It is formulated as a mixed integer problem which optimizes a cost that includes reconstruction costs of network components, operational costs, and penalties incurred for unsatisfied demand, while taking constraints on flows of commodities into account. In general mixed integer program are known to be NP-hard except in special cases. For the td-NDP this means problems become exponentially harder as the size of the network to be reconstructed increases; therefore, it is common practice to break up the recovery into time windows of length $T$, and find the locally optimal solution in each window. A formal definition of td-NDP is provided in the Methods section.

To uncover the key driving factors and properties of the recovery process, we apply the td-NDP to a representative example, the power grid of Shelby County, Tennessee, which consists of $9$ suppliers and $37$ consumers (with $14$ junction nodes where $d_{i}=0$), connected by $M=76$ transmission lines. Network topology and necessary parameters where obtained from Refs. [22, 23]. Figure 1a shows the total repair cost as a function of time as we perform the td-NDP with $T=5$ on a network that was initially completely destroyed, with cost broken down by the type of expense. We see that deficit cost (i.e., the penalty accrued for unsatisfied demand) is overwhelmingly dominant and exponentially decreasing throughout the initial stages of recovery. Investigating how this impacts the growth of the components, Figure 2a shows the commodity deficit or supply of each connected component throughout the recovery process, with circle sizes representing the component size, and colors representing if a component is over- (blue) or undersupplied (red). We see that the td-NDP process results in many small components initially, with relatively small surplus/deficit, and only towards the end of the process all components are joined. This is consistent with islanding techniques discussed in engineering practice. Our goal is to develop a simple and computationally efficient model of the recovery process that captures these key features.

Previous topology-based recovery processes prioritize the rapid growth of the Largest Connected Component (LCC) [11, 12, 10, 13]. Models vary in the details, such as the type of failure (random, localized, catastrophic) and additional secondary objectives (such as prioritizing nodes based on population), but the metric for the quality of the solution is directly related to how quickly the LCC grows.

As a representative example of topology driven recovery strategies, we implement an Achiloptas process using a selection rule that maximizes the sum of the resulting component, which we refer to as LCC percolation. In this process, we randomly select $M>0$ candidate edges out of the set of damaged edges at each discrete timestep. We then examine the impact that each individual edge would have and select the edge that, when added, creates the largest connected component. More specifically, let $s_{i}$ denote the size of the component to which node $i$ belongs. If nodes $i$ and $j$ belong to separate components, repairing edge $(i,j)$ creates a component with size $S_{ij}=s_{i}+s_{j}$; if they belong to the same component, the size of the component does not change and we set $S_{ij}=0$. Out of the $M$ candidate edges, we repair the one that maximizes $S_{ij}$; if multiple candidate links have the same maximal $S_{ij}$, we select one of them uniformly at random. If $M=1$, the process is equivalent to traditional percolation. If $M=E$, the process is largely deterministic, we always repair an edge that is optimal with respect to the selection criteria.

We now apply LCC percolation as a model for the recovery of the Shelby County power grid and compare the results to the benchmark td-NDP process. Figure 2b shows that if we use growth of the LCC as our objective, the LCC, represented by the largest circle, grows rapidly throughout the process as expected. However, the deficit/surplus of this component fluctuates greatly. As indicated by the magnitude of total commodity deficit $D$, i.e., the total unsatisfied demand in the network, in Figure 1b, such a recovery algorithm is costly and leaves large portions of the grid without power until the final steps of the recovery algorithm. To conclude, we find that a recovery process based on quick growth of the LCC is neither cost efficient nor effective for satisfying consumer’s demand as quickly as possible. As a result, we do not consider this algorithm for further studies.

We have shown that the key driving factor in recovery processes is reduction of the total commodity deficit, i.e., the total unsatisfied demand, and that optimizing LCC fails to capture this. To capture the essence of real recovery strategies, we propose a competitive percolation process which we refer to as recovery percolation that, instead of optimizing for LCC growth, aims to directly reduce the unsatisfied demand. In addition to network topology, this recovery process also takes into account the net demand or production of the individual nodes.

We define $D_{i}$ as the commodity deficit of the connected component to which node $i$ belongs. We assume that capacity constraints of the transmission lines are sufficient and thus do not limit the flow of commodities during the recovery process, a common practice in infrastructure recovery literature [14, 22, 23]. Therefore the commodity deficit of a component is the sum of demand or supply of individual nodes belonging to the component, i.e., $D_{i}=\sum_{j\in\mathcal{C}_{i}}d_{j}$, where $\mathcal{C}_{i}$ is the set of nodes belonging to the component containing node $i$.

We use commodity deficit of the components as a selection criteria for the competitive percolation model to account for the goal of balancing supply and demand. Similar to LCC percolation discussed above, we randomly select $M>0$ damaged candidate edges, from which one is chosen to be repaired and added to the network at each timestep as follows. We first consider how much demand would be met by adding each of the $M$ edges individually to the network. More specifically, if nodes $i$ and $j$ belong to components such that $D_{i}D_{j}<0$, then repairing edge $(i,j)$ reduces the total commodity deficit by $\Delta D=\min(\lvert D_{i}\rvert,\lvert D_{j}\rvert)$; if $D_{i}D_{j}\geq 0$, then there is no commodity deficit reduction, i.e., $\Delta D=0$. Out of the $M$ candidate edges, we repair the one that maximizes $\Delta D$; if multiple candidates have the same maximal $\Delta D$, we repair one of them chosen uniformly at random. As $M$ shrinks and approaches $1$, the process becomes more stochastic; while if $M=E$ the process is largely deterministic, as we always select an edge that is optimal with respect to the selection criteria.

Figure 1b shows that for the Shelby County power grid the total commodity deficit during the recovery percolation for $M=E=76$ well approximates the td-NDP, especially at the beginning of the recovery process when costs are much higher. We also see that even when $M=10$, corresponding to only $\sim 13\%$ of the total edges considered at each timestep, the approximation remains very effective. Figure 2c shows similar dynamical behavior in recovery percolation as in the td-NDP solution (cf. Figure 2a): larger components delay formation, and tend to have smaller commodity deficit.

We have shown that recovery percolation follows our benchmark td-NDP solution closely on a real-world topology. We also observed that the growth of connected components is suppressed via recovery percolation as compared to LCC percolation. To understand this behavior, we study large systems with $N=10^{4}$ nodes and we allow potential edges to exist between any node pair, removing underlying topology constraints. Note that the td-NDP process is intractable for networks of this size.

Figure 3 (main figure) shows the growth of the LCC for a range of $M$ values. For $M=1$, the model reduces to random percolation which has a second-order phase transition at $t/N=0.5$, and above this critical point the LCC becomes proportional to $N$. As we increase $M$, the apparent transition point shifts to higher values and approaches $1$, indicating that the appearance of large-scale connectivity is suppressed; however, once the transition point is reached, the growth of the LCC becomes increasingly abrupt. This observation is analogous to explosive percolation, where links are chosen to be constructed explicitly to delay component growth [29]. In contrast, in recovery percolation it is an indirect consequence of the restoration strategy.

To understand the underlying mechanism of component growth, we plot the average component sizes and their corresponding average undersupply at various points during the reconstruction process in the bottom row of plots in Figure 3. Note that average oversupply behaves in a similar manner, but is omitted for clarity. The left column of plots in Figure 3 shows the same quantities such that the size of the LCC is fixed. The main trend we observe at any given point in time is that for large enough $M$ there is negative correlation between component size and undersupply, and this correlation becomes stronger as $M$ increases. This means that as components grow their commodity deficit is reduced and therefore the likelihood of further growth is also reduced, ultimately suppressing the appearance of large scale connectivity.

The observed two features also describe islanding, an intentional behavior in resilience planning and recovery in real-world power grids. This islanding behavior is already observed in early stages of the restoration process, becoming more apparent as $t$ approaches the transition point.

So far we investigated the recovery process on an underlying graph without topological constraints. We also wish to analyze more realistic networks and turn our attention to synthetic power grids [31, 32]. This allows us to systematically investigate how typical structural features of power grids affect the efficiency of the recovery process.

Power grids are spatially embedded networks, and physical constraints limit the maximum number of connections a node can have; their degree distributions, therefore, have an exponential tail, in contrast to many complex networks that display high levels of degree heterogeneity. Power grids typically have average degree $\left\langle k\right\rangle$ between $2.5$ and $3$ [32, 33]. An important requirement of power grids is structural redundancy, meaning that the failure of a single link cannot cause the network to fall into disconnected components. A network without redundancy has tree structure, has average degree 2 and all node pairs are connected through a unique path. Any additional link creates loops and improves redundancy. Structural redundancy can be characterized locally by counting short range loops. For example, power grids have a high clustering coefficient $c$, typically ranging between $0.05$ and $0.1$ [32]. The algebraic connectivity, denoted by $\lambda_{2}$, is the second smallest eigenvalue of the network’s Laplacian matrix and captures a measure of global redundancy: it is related to the number of links that have to be removed in order to break the network into two similarly sized components, with high value corresponding to high redundancy. The exact value of $\lambda_{2}$ depends on system size, where for a given number of nodes $\lambda_{2}$ is minimal for tree structure, and monotonically increases as further links are added [34].

To generate networks that exhibit these features, we use a simplified version of a practical model developed by Schultz et al. [32]. The model generates spatially embedded networks mimicking the growth of real-world power grids. The process is initiated by randomly placing $N_{0}$ nodes on the unit square and connecting them with their minimum spanning tree. To increase redundancy, $qN_{0}$ ($0\leq q\leq 1$) number of links are added one-by-one, such that each link is selected to minimize the redundancy-cost trade-off function

where $i$ and $j$ are two nodes not connected directly, $d_{\textrm{net}}(i,j)$ is their shortest path distance in the network, and $d_{\textrm{euc}}(i,j)$ is their Euclidean distance. The $r\geq 0$ parameter controls the trade-off between creating long loops to improve redundancy and the cost of building power lines. After the initialization, we add $N-N_{0}$ nodes through a growth process. In each time step a new node is added: with probability $1-s$ the node is placed in a random position and connects to the nearest node; with probability $s$ a randomly selected link $(i,j)$ is split and a new node is placed halfway between nodes $i$ and $j$ and is connected to both of them. To increase redundancy, in each time step an additional link is added with probability $q$ connecting a randomly selected node $i$ to node $j$, such that $f(i,j)$ is minimized. Finally, a fraction of nodes $p_{\textrm{s}}$ are randomly selected to be suppliers, the rest are assigned to be consumers.

Changing parameters $q$, $r$, and $s$ allows us to systematically explore how these parameters impact the structure of these model power grids (Fig. 4): $q$ controls the average degree $\left\langle k\right\rangle=2(1+q)$ and adds redundancy to the network; $r$ controls how loops are formed, where small $r$ favors short distance connections leading to high $c$ and low $\lambda_{2}$, while large $r$ favors long loops leading to low $c$ and high $\lambda_{2}$; and $s$ increases typical distances in the network, lowering both $c$ and $\lambda_{2}$.

Here, we define our benchmark model for network recovery: the time-dependent network design problem (td-NDP). Our version follows the more general formulation developed by Gonzalez et al. [22, 23, 30]. The td-NDP takes a graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$, where $\mathcal{N}$ is a set of nodes, and $\mathcal{E}$ is the set of edges connecting nodes. At the beginning of the recovery process the td-NDP uses the destroyed graph, $\mathcal{G}^{\prime}=(\mathcal{N}^{\prime},\mathcal{E}^{\prime})$, where $\mathcal{N}^{\prime}$ and $\mathcal{E}^{\prime}$ represents the nodes and edges that are not functioning, respectively. The objective function (cf. Equation (2a)) minimizes the total reconstruction cost over a given time domain $\mathcal{T}$ with $t\in\mathcal{T}$, which includes the cost to repair nodes, $q_{it}$, cost to repair edges, $f_{ijt}$, cost of flow on each edge, $c_{ijt}$, and oversupply and undersupply penalties for each node, $M^{+}_{it}$ and $M^{-}_{it}$. These costs usually depend on multiple factors, such as the level of damage, the type and size of the components to be restored, their geographical accessibility, the amount of labor and resources required, and the social vulnerability of the affected areas, among others [22, 38, 39]. To keep track of demand satisfaction, each node also has a supply capacity (demand if negative), $b_{it}$. In the most general formalization of the problem, node supply $b_{it}$ can depend on time $t$, but in this paper we only consider constant values. The variables $\delta^{+}_{it}$ ($\delta^{-}_{it}$) account for oversupply (or undersupply) of node $i$. We refer to the sum of the absolute values of oversupply and undersupply ($\left|\delta^{+}_{it}\right|+\left|\delta^{-}_{it}\right|$) as the commodity deficit of node $i$. The td-NDP includes as decision variables the amount of flow on each edge, $x_{ijt}$, whether or not a node $i$ [edge $(i,j)$] is chosen to be recovered at timestep $t$, $\tilde{w}_{it}$ ($\tilde{y}_{ijt}$), and whether or not a node $i$ [edge $(i,j)$] is functional at timestep $t$, $w_{it}$ ($y_{ijt}$). Constraints 2b- 2o are imposed to ensure that conservation of flow properties are held and that only recovered and functional nodes can produce or consume flow.

The td-NDP formulation is a mixed integer program, which has been shown to be, in general, NP-hard (and becomes exponentially harder as $\mathcal{T}$ and $\left|\mathcal{G}^{\prime}\right|$ grows). The number of variables and constraints also become larger as the input graph becomes larger. For many reasonable size problems, computing a global optimal (i.e., where $\mathcal{T}$ contains the entire time horizon for recovery) is intractable. Therefore, heuristics are used to restrict the size of $\mathcal{T}$ by dividing the total recovery time into smaller windows, and finding the locally optimal solutions within these windows [22]. It has been shown that such heuristic finds solutions very close to the optimal; however, the computational complexity is still relatively high as a result of the underlying mixed-integer program.

For our computational experiments, we generate our demand distribution by following the load distribution of the European power grid [39]. This dataset was chosen due to its large system size ($N=1463,E=2199$) and its high resolution. Our goal is not to identify the true analytic form of the load distribution, but to generate statistically similar samples through bootstrapping. We found that an exponentiated Weibull distribution of the form $f(x,a,c)=ac(1-\exp{(-x^{c})}^{(a-1)}\exp{(-x^{c})}x^{(c-1)})$, where $a=3.59$ and $c=0.8$ well approximates the features of the demand distribution. Suppliers’ capacities are uniformly distributed to balance the total demand. We also get our ratio of suppliers ($0.3$) to consumers ($0.7$) from this dataset.