Steady state

At each time step, we measure the temporal load correlation, defined as:

where is the load averaged over all nodes at time t. The brackets represent an average over the last ten consecutive time intervals of length 100T p , i.e., t′ = t − n100T p with n = 0, 1, …, 10.

Starting with all loads equal to zero, we observe that decays in time towards a steady state in which the dissipation balances the total incoming load. When drops below 1% of the relative standard deviation of the loads within the network, we assume that the steady state is reached. In the steady state, the load at each node has a well-defined average value with small fluctuations. The spatial distribution for a given realization of voltages in the European power grid is shown in Fig. 1. The size of the dots represents the average load over a time window of 1000T p measured in the steady state. Most nodes accumulate negligible load (ρ ≈ 0), while surprisingly a small fraction of the nodes are overloaded (ρ > 10ρ 0 , where ρ 0 is the magnitude of each perturbation).

Figure 1 Spatial distribution of load on the European power grid in the steady state. Distribution is shown for a particular realization of the voltage distribution. The size of the nodes corresponds to the average load allocated on them. The smallest size refers to zero load, where the largest one to the maximum load. Full size image

The observed load exemplarily shown in Fig. 1 corresponds to one realization of the voltage distribution and thus for one configuration of the direction of the edges. Assuming small temporal changes in the network (e.g., fluctuations of voltage in the power grid or number of cars entering a road junction), the direction of the edges changes in time. Thus, we also consider different realizations of the voltage distribution and, for each realization, we determine the steady state load distribution. Figure 2 shows the relative load distribution averaged over 5000 realizations. In order to compare the load distribution of different networks (power grid, Watts–Strogatz and scale-free networks), loads in each curve are rescaled by the magnitude of each perturbation (ρ 0 ). The strongly inhomogeneous behavior of the steady-state load seen in Fig. 1 is also visible in the load distribution. The distributions are bimodal defining two different types of nodes: those with a negligible load compared to the perturbation (ρ < ρ 0 ) and those with a larger load (ρ > ρ 0 ). The latter ones are typically overloaded in the steady state, suggesting that the incoming perturbations interfere constructively at them.

Figure 2 Relative load distribution in the steady state for different network topologies. The considered topologies are: power grid (green dots), scale-free (red squares) and Watts–Strogatz networks (blue triangles). The power grid and the scale-free network have N = 1254 nodes and M = 1811 edges, while the Watts–Strogatz network has N = 1254 nodes and an average degree of 〈k out 〉 = 2. All loads are in units of the amplitude of the perturbation. Each curve is an average over 5000 voltage realizations and, in the case of the model networks, also an average over 100 different networks. The magnitude of the standard deviation of the curves is comparable to the size of the symbols. Full size image

The plots for the two network models (Watts–Strogatz and scale-free) in Fig. 2 are obtained for networks with the same number of nodes as the power grid. The average degrees are also kept close to the power grid, with the same number of edges in the scale-free network and 〈k out 〉 = 2 in the Watts–Strogatz graph. In both cases, a bimodal distribution is also observed. The power-grid network is constructed from real data and its size corresponds to the real network size. Thus, a finite-size study is not possible. Yet, in the case of the model networks one can systematically study the effect of the network size on the load distribution. Figure 3A shows the load distribution for Watts–Strogatz networks of different network sizes. The majority of the nodes (more than 90%) has always a negligible load, while the load of the remaining nodes follows a broad distribution, characterized by a decay in the relative frequency with increasing load and a cut-off for values of load close to ρ 0 . The bimodal distribution smoothens out for larger network sizes. For scale-free networks the qualitative picture is slightly different. As shown in Fig. 3B, for all network sizes, one observes two power-law regimes, with a crossover at ρ 0 . Nevertheless, note that for both network models, there is always a significant fraction of nodes (around 10%) with a non-negligible load. The load value of the cutoff suggests that, at large system sizes, consecutive shock waves that enter the network are separated so that they attenuate their amplitude before being able to interfere.

Figure 3 Size-dependence of the relative load distribution in the steady state for two network models. The considered topologies are: (A) Watts–Strogatz and (B) scale-free networks with 〈k〉 = 2. Insets show the respective data collapse, where ρ and P denote the load and the relative frequency, respectively, N is the size of the network. Loads are divided by the magnitude of the applied perturbation to ensure comparability. Each data is averaged over at least 100 different graphs and 100 voltage realizations. The breaks in the distributions around ρ = 10−3 are due to the logarithmic binning. Full size image

The specific nodes that exhibit these high load values typically change from realization to realization. However, after averaging over different voltage distributions, we still find some nodes which are consistently overloaded. For each distribution of voltages, we classify as “overloaded nodes” the ones with a load at least ten times larger than the average. We define vulnerability of a node as the probability that it is an overloaded node. Figure 4A shows the spatial distribution of vulnerability in the European power grid where the color and size of the nodes denotes their vulnerability. The vulnerability of white nodes is lower than 0.1%, while the one of the black nodes is larger than 5%. All the other nodes (about 30% of the nodes, in gray color) have a vulnerability between 0.1% and 5%. In comparison, the highly vulnerable nodes are at least 50 times more frequently exposed to large incoming fluxes. In the case of random perturbations, vulnerable nodes are more likely to fail or be congested. It is therefore crucial to identify these nodes to improve their capacity and to mitigate the risk of failure. Figure 5 shows the vulnerability distribution corresponding to the map in Fig. 4A and to other network topologies. In the case of the Watts–Strogatz and scale-free networks, a vulnerability distribution for networks of size N = 105 are presented.

Figure 4 Spatial distribution of two properties on the power grid network. (A) Vulnerability of the nodes obtained solving the Burgers equation (probability of having in a voltage realization). (B) Node-basin sizes averaged over different voltage realizations. Color and size indicate the strength of the corresponding property: black corresponds to large vulnerabilities (> 5%) and basin sizes (> 10), while white nodes have negligible vulnerabilities (below 0.5%) and small basin sizes (close to 1). The data are the average over 5000 voltage realizations. Full size image

Figure 5 Distribution of vulnerabilities for three different network topologies. The distribution of the vulnerability (i.e., the probability of having a load ten times larger than the average) is shown for power grid (green circles), scale-free network (red squares) and Watts-Strogatz graph (blue triangles). The data for the model networks are the average over 100 network realizations. Full size image

Identifying vulnerable nodes

Next we will introduce a simple topological property of the nodes to identify the vulnerable spots without solving the dynamics of the Burgers equation. According to the Burgers equation, each node sheds the incoming shock wave among its out-edges and the total load agregated at the node at a time t is the sum of incoming loads. Hence, if we track the path of a given shock wave, it is fragmented at each node with multiple out-edges and will stop at any node that does not have any out-edges. Assuming that the time average of the load at a node is proportional to the number of incoming shock waves, we determine the basin corresponding to that node. For this purpose, let us consider one realization of the edge directions (see Fig. 6). For a given node (the one in red in Fig. 6), the corresponding basin is defined as the smallest subgraph of the network containing this node (as the sink), which is connected to the rest of the network by outgoing edges. Following the procedure as illustrated in Fig. 6, we go through the nodes following the opposite direction of the in-edges, starting from the red node (Fig. 6A) and add nodes to the basin that has only out-edges at the end of the process (marked by the blue region in Fig. 6B). The resulting subgraph is the basin of the red node and the contribution of the load of the nodes in the basin is simply the inverse of their out-degree (Fig. 6C), except for the initial (red) node, that contributes to the load with unity. This choice of the contribution of the nodes in the basin is based on the fact that Eq. (7) conserves the flux at each node. For simplicity, we assume here that the amplitude of the shock waves leaving a node is on average the same for each out-edge.

Figure 6 Calculation of the node-basin size. (A) The basin corresponding to the red node is considered. (B) We determine the smallest subgraph containing the red node and having only out-edges to the rest of the network. This is equivalent to a breadth-first search traversing the fraction of the graph reached through only in-edges. (C) When the basin is determined, each node in the basin contributes to the red node's basin size by the inverse of its out-degree. The contribution of the red node (i.e., the sink of its basin) is unity. Full size image

We calculated the size of the basins for each node, defined as the sum of contributions from all nodes inside the basin of the corresponding node. This basin size (which is determined for one voltage distribution) is then averaged over different voltage configurations. The resulting basin size distribution is depicted in Fig. 4B for the power-grid network. One sees that for this network, the distribution of the node-basin size is very similar to the distribution of the vulnerability. A quantitative comparison of the two properties can be given by their correlation. Thus, we plot the rank-rank scatter plots in Fig. 7 of the indices of the nodes after being sorted in ascending order by vulnerability and node-basin size. The plots are the average over 5000 voltage distributions and over 100 different topological realizations of model networks.

Figure 7 Scatter plots of the ranked vulnerability and node-basin sizes for different network topologies. The corresponding network models are: (A) power grid, (B) scale-free and (C) Watts–Strogatz network. The axes on the plots indicate the corresponding ranked property and the color denotes the number of nodes for which the two ranks were identical. In other words, the color of a dot at (x, y) corresponds to the number of nodes with vulnerability rank of x and node-basin size rank of y (see the colorbars). All plots are the average over 5000 voltage realizations. In the case of model networks, 100 different realizations are considered. Full size image

The corresponding product-moment correlations ρSp of the ranks are given above the plots, showing strong correspondence between the ranked vulnerability and node-basin size for the power grid and scale-free networks. The product-moment correlation is,