The reliable operation of supply networks is crucial for the proper functioning of many systems, ranging from biological organisms such as the human blood transport system or plant leaves to man-made systems such as power grids or gas pipelines. Whereas the failure of single transportation links has been analysed thoroughly, the understanding of multiple failures is becoming increasingly important to prevent large scale damages. In this publication, we examine the collective nature of the simultaneous failure of several transportation links. In particular, we focus on the difference between single link failures and the collective failure of several links. We demonstrate that collective effects can amplify or attenuate the impacts of multiple link failures—and even lead to a reversal of flows on certain links. A simple classifier is introduced to predict the overall strength of collective effects that we demonstrate to be generally stronger if the failing links are close to each other. Finally, we establish an analogy between link failures in supply networks and dipole fields in discrete electrostatics by showing that multiple failures may be treated as superpositions of multiple electrical dipoles for lattice-like networks. Our results show that the simultaneous failure of multiple links may lead to unexpected effects that cannot be easily described using the theoretical framework for single link failures.
1.Introduction
The failure of links can impede the operation of supply networks leading to potentially critical events [1–3]. Cascading failures in power grids can cause power outages affecting millions of households [4–7], and embolism in humans and plants may result in strokes [8, 9] or leaf death [10]. Such events are typically caused by the failure of one or few transportation links [5].
In this article, we analyse the impacts of multiple link failures in linear flow networks. We focus on electric power grid operation while the mathematical results apply to a variety of networks. In case of the power grid, cascades triggered by a single failure are in most cases prevented by the transmission grid operators. This is typically achieved by running the grid N − 1 secure which means that a single failing transmission or generation element does not prevent stable operation of the power grid [11]. However, grid operators are now encouraged to also take specific dangerous N− 2 contingencies into account [12] due to an increased vulnerability of the grid. For instance, an increased risk of extreme weather events caused by climate change raises the risk of several transmission elements failing, thus leading to power outages [13]. More specifically, an increase in correlation between transmission line outages, e.g. through more extreme weather events, was recently associated with an increased risk of cascading outages, raising the relative contribution of multiple link failures to such cascades [14]. In addition to that, future power systems with a high share of renewable energy sources will have to transport power over long distances using long transmission lines, thus also increasing the risk of dependent link failures.
Many computational approaches towards studying and classifying N − 2 outages have been developed in order to identify contingencies that result in additional overloads [15–18]. Nevertheless, such outages still lack a fundamental theoretical understanding. Basic mathematical tools have been developed extending the concept of Line Outage Distribution Factors (LODFs) originally used for single link contingencies [11] to include multiple link failures thereby allowing for a mathematical description of these contingencies [19–21]. These tools demonstrate that the nature of multiple outages may be fundamentally different from the outage of single lines, thus making a direct transfer of understanding and intuition developed for single link failures [22] difficult. In particular, multiple outages can enhance or attenuate each other in a counterintuitive manner. Which topologies drive such phenomena and which ones prevent them from happening is at present not fully understood.
In this article we analyse the collective nature of N − 2 failures in linear flow networks, which describe different systems including AC power grids in the linear approximation [11]. We demonstrate that two simultaneous failures can cause a disturbance that strongly differs from the sum of the disturbances induced by individual failures; they can amplify or attenuate the flow changes in a grid. In addition to that, we introduce a predictor which allows us to understand under which circ*mstances these collective effects play an important role and when they can be neglected. We then apply the predictor to different test grids and reveal its performance in forecasting collective effects for multiple link failures quantitatively, outperforming also distance measures proven to be good predictors in the case of single link failures. Finally, we extend on previous work [22] that successfully established an analogy between flow rerouting after single link failures and the fields of electromagnetic dipoles in regular grids by demonstrating that flows after multiple link failures may be treated as a superposition of multiple individual dipole fields in such grids in the continuum limit.
2.Link failures in linear flow networks
2.1.Fundamentals of linear flow networks
Linear flow networks describe the operation of various types of systems including AC power grids [11, 23, 24], DC electric circuits [25–27], hydraulic networks [28, 29], and vascular networks of plants [30]. In such networks, the flow 
over a link (m, n) is assumed to be linear in the potential or pressure drop along this link,
In this article, we focus on applications to AC power grids, where 
is the real power flow, 
is the voltage phase angle at node n and 
is proportional to the link's susceptance. We assume that the susceptance is independent of the direction of the link, 
, and that it vanishes if no link (m, n) exists. In this context, the linear description is commonly referred to as the DC approximation due to its formal equivalence with DC resistor networks [11, 23, 24]. This approximation is typically good for transmission grids with weak link loading, see [23] for details. In hydraulic or vascular networks, 
denotes the pressure at node n while the transmission capacity bmn depends on the geometry of a pipe or vein [28–30]. The flows are subject to the continuity equation which means that at each node of the grid the sum of the network flows must equal the inflow to the grid;
The inflow Pm is positive if a current, power, or fluid is injected to the node and negative if it is withdrawn from the node. In the following, we assume these in- and outflows to be balanced, 
.
Equations (1) and (2) fully describe the state and the flow of the network—up to a constant phase shift applied to all voltage phase angles—once the link parameters bmn and the injections Pm are given. We introduce a compact vectorial notation summarising the nodal potentials or voltage phase angles in the vector 
and the nodal injections in the vector 
. Here and in the following sections, the superscript '⊤' denotes the transpose of a vector or matrix. We further label all lines in the grid by 
and fix an orientation for each link. Summarising all link flows in the vector 
, equation (1) reads as
with the diagonal matrix of link strengths 
. Furthermore, we made use of the node-edge incidence matrix 
commonly used in graph theory. It establishes a correspondence between the nodes in the graph and the edges connecting them and has the components [31]
In the following, we use this matrix to assign an (arbitrary) orientation to each link in the network such that 
. Using the node-edge incidence matrix we can further rewrite the continuity equation (2) in the compact form
The matrix 
is commonly referred to as the nodal susceptance matrix in power engineering. Mathematically, 
is a weighted Laplacian matrix [31, 32] with components
For a connected network, this matrix has one zero eigenvalue 
with eigenvector 
such that 
. For this reason, it is not invertible. However, the matrix inverse appears naturally in many different contexts involving the spreading of failures in networks. To be able to nevertheless study these processes, one typically considers the Moore–Penrose pseudoinverse 
which has properties similar to the actual inverse, see e.g. [33, 34] for details. We are now ready to extend the notation to cover link failures as well.
2.2.Single and double link failures
Assume that a single link k in the network fails, thus losing its ability to carry any flow. Since the network after the failure is still subject to the continuity equation (2), the failure will cause the flows on other links to change to account for the remaining necessary transport. Assume that the new flows are given as 
, where 
is the vector of flow changes and 
is the preoutage flow. In general, we will use the superscript (0), i.e. round brackets, to indicate a flow before an outage and the superscript 
, i.e. square brackets, to indicate a flow after the failure of link k. Then the new vector fulfils the continuity equation (3),
In power engineering, the changes of flows are typically captured in a matrix of LODFs whose element Ll,k describes the flow changes monitored on a link l after another link k fails. Suppose that 
is the flow on link k before the outage. Then the LODF is defined by its elements [11]
For consistency, the effect of an outage of a line k onto itself is typically defined as 
. Importantly, the LODF may be expressed in purely algebraic form using the (pseudo-) inverse of the graph Laplacian 
[11],
Here, we abbreviate the line susceptance 
of a link 
by bl. Furthermore, we defined a vector 
that characterises a link k=(r, s) and has the entries +1 at position r, −1 at position s and zero otherwise. Using the standard basis vectors in 
, this vector may be written as 
. In power engineering, the expression in the numerator is also referred to as the Power Transfer Distribution Factor (PTDF) [11]. A PTDF between links l and k is calculated as
and describes the flow changes on link l upon a power transfer from one end of link k to the other one. PTDFs are typically defined for power injections and withdrawals at arbitrary nodes in the network [11]. In the context of link failures, however, it is useful to restrict them to power injection and withdrawal taking place at the two ends of a link. For this reason, power injection vectors 
have to correspond to the columns of the incidence matrix in our setup such that 
, where 
is again the standard basis vector.
On the other hand, the link failure may also be described on the nodal level. If we collect all changes in voltage phase angles after the failure of link k in the vector 
, denoting phase angles after the failure by 
, a Poisson-like equation describing the outage in terms of the phase differences may be derived [22],
Here, 
is the Laplacian of the network after removal of link k and 
is the pre-outage flow on link k. This equation was studied in the past in different settings [35, 22, 36, 33]. The failure of single links is thus comparably well understood [22, 37, 38], whereas the simultaneous failure of multiple links was not yet studied to the same extend on a theoretical level.
For this reason, we now turn to the case of multiple link failures and derive an expression for the flow changes on the remaining lines in the grid. We will focus on the case of two outages for now, but extend the results to more than two outages in section 4. Naively, we could just superpose the flow changes caused by the two individual outages as described by equation (5). Assuming that two arbitrary links o and k fail, this naive approach yields the following expression for the flow changes on link l,
However, this approach neglects the effect of the outage of link k on link o and vice versa. To arrive at the correct formula, we need to consider this interaction as follows; the outage of link o changes the flow on link k by
where we denote by 
the effective flow on link k as a result of the simultaneous failure of o and k. Note that this is not the actual flow on link k, but rather an effective quantity used for calculation here since link k fails and thus carries no flow after the outage. An analogous expression holds for the effect of the outage of link k on link o. Inserting these corrected flows into equation (8), we arrive at the following result for the flow changes on link l [17],
Here, the superscript [o,k] is used to indicate the flow or flow changes on a link after the failure of both links o and k. Finally, expanding this expression results in the following equation encoding the collective flow changes in a compact form,
The resulting expression for 
is different from the simple linear combination (8) due to the interaction of the two failing lines o and k, which is encoded in the matrix in the centre. More precisely, the collective effects are governed by the LODFs of the interacting lines Lo,k and Lk,o forming the off-diagonal elements of the matrix. The inverse matrix in this formula may be calculated as
Thus, the flow changes on a link l are given by
where 
is a symmetric prefactor.
The equation describing flow changes after the failure of two links thus differs from a naive overlay of the two individual outages. In the following sections, we will demonstrate in which cases these collective effects resulting from the interaction of both outages are important and in which cases they may be neglected.
2.3.Elementary examples
In this section, we elucidate different elementary examples that describe possible interactions between the individual outages and allow us to understand the role played by collective effects in more detail.
2.3.1.Amplifying single outages
To start with, we present a case where the naive superposition of two individual outages underestimates the collective effects such that
for some link l. An elementary example of a network where this is happening is shown in figure 1 where the topology is given by a network consisting of N=9 nodes and M=12 links connecting them in a square grid. This initial setup is shown in panel a. Panels b and c illustrate the flows on each link (numbers on links) after the failure of two different individual links (coloured red). The bold number indicates the link with maximal flow for the given setup. Each single outage leads to a maximal flow on the top right link of 
and 
after the failure shown in panels b and c, respectively. Naively, we would thus expect the failure of both links to lead to a flow of
by simply superposing the two individual outages. The actual outage of both links, however, results in a much larger flow on the link that reads as
Whereas the two individual outages separately lead to a moderate increase in flow on the most heavily loaded link, their interaction results in a larger flow potentially reaching the link limit. If the flow on all links was limited to, say, 
, the naive superposition would thus predict no overload caused by the two link failures, whereas in fact, the link loaded maximally in panel d breaks down in this case. This example demonstrates that the results predicted by the theory of single link failures may differ drastically from the correct calculation that takes into account the collective effects.
2.3.2.Effect of individual outages exceeding simultaneous outage
In addition to the effect presented in the last section, it may also happen that an additional outage is beneficial for the most heavily loaded link. A minimal example is shown in figure 2. In this setup, we have a unit inflow of power at the centre left node (coloured yellow), P=+1, and a unit outflow at the bottom right node (coloured purple), P=−1, whereas all other nodes neither consume nor create power, P=0. The initial setup is shown in panel a, where again (absolute) flows are indicated as numbers on the edges as well as colour coded. Panels b and c show the network after the single outage of two different links (coloured red). The edge with the highest absolute flow is indicated by a bold face number in both cases. In panel b, we have 
whereas for panel c the maximum flow reads as 
. The situation after the simultaneous outage of both links is shown in panel d. The edge with maximum flow now carries an absolute flow of 
, i.e. the collective effects attenuate the flow on the edge with the highest flow compared to each individual outage. This effect can be seen as a realisation of Braess' paradox [39–42] since the outage of an additional link is beneficial in terms of the maximum absolute flow in the network for each individual outage. Hence, cascades of failures may in some situations be prevented by the intentional removal of a second, carefully chosen link after a first link failure threatens stability [1, 43].
2.3.3.Sign inversion through double outages
In this section, we discuss a highly surprising phenomenon that appears in the case of multiple interacting outages: the collective effects may dominate in such a way that purely collective effects can cancel (as shown in the previous section) or even overcompensate the direct effects of individual link failures such that the flow changes resulting from the failure of both links have a different sign compared to the flow changes after each of the individual failures.
To study this in more detail, we will use the following notation in this section. Suppose that links o and k fail and flow changes are monitored on link l. Then we denote by 
or 
the flow changes on link l when link o or k fail, respectively. For the simultaneous outage of both links o and k, we denote by 
the actual flow change on link l. With this notation at hand, we will construct examples where flow changes caused by individual link failures 
and 
have the same sign, but the collective flow change 
has the opposite sign.
A small example of a network where such a situation occurs is shown in figure 3. In the initial setup, there is a small flow Fl=0.5 on link l (a, bold face number). For the failure of two individual links o and k shown in panels b and c, respectively (red links), the flow on this link is amplified showing positive flow changes in both situations, 
and 
, respectively. However, if both links fail simultaneously, the overall flow change has the opposite sign, 
, thus even inverting the direction of the flow 
with respect to both, the individual setup and the situation after the failure of each individual link.
In the next paragraph, we explain this surprising phenomenon on a theoretical level in more detail. For simplicity, let us assume the flow changes due to the individual outages to be positive and the flow change in the case of a simultaneous outage to be negative,
Plugging in equation (9), we can cast these three conditions into the following form based on LODFs,
We can assume that both the initial flows on the failing links and the LODFs between the failing links and the reference link are positive, 
, without loss of generality–this can always be accomplished by redefining the orientation of one or both of the initial flows. Then the right-hand side of condition (1) is negative. Hence, the condition can only be satisfied if the left-hand side is negative as well, which requires that the mutual LODFs are both negative, 
because they always have the same sign (see appendix A). Note that we do not consider cases where both of the LODFs Lk,o and Lo,k are equal to (minus) one, thus keeping 
finite. For notational convenience, let us now introduce positive constants α and β defined by the following quotients;
We thus incorporated the whole dependency on the link monitoring the flow changes l into the purely topological constant β, whereas any dependency on the flows, i.e. the specific power injections, is incorporated into α. Dividing the first condition (1) by its right-hand-side, we arrive at the following inequality
This inequality can only be fulfilled if there is a strong degree of heterogeneity between α and β, i.e. the ratio of initial flows on links k and o differs strongly from the ratio of their LODFs with respect to link l. Furthermore, the mutual LODFs Lo,k and Lk,o need to be both large in magnitude in order to reduce the size of the denominator compared to the numerator. We will see in the next sections that strong mutual LODFs also imply strong collective effects caused by the simultaneous failure of links k and o.
Condition(11) can be simplified drastically if 
which may be realised e.g. through a very small initial flow on link o compared to link k, such that 
. The above inequality then reduces to
We will now demonstrate how to design a network where this inequality is satisfied. The construction works as follows; we design a network topology where two links o and k influence each other heavily (measured in terms of LODFs) while a third link l is influenced very differently by each of the links. Formally, Lo,k and Lk,o both need to be comparatively large while Ll,o and Ll,k should be very different in size, thus leading to a small value of β. To observe flow sign reversal, we choose the power injections 
in such a way that the flow on link k is much larger than the one on link o, which results in 
. Furthermore, our choice of power injections also needs to make sure that flow changes on link l are positive 
. The resulting network fulfils the three conditions on the flow changes given in the inequalities(10).
Indeed, we can find networks where the inequalities are fulfilled as shown in figure 3. The parameters in this case are given by 
, 
, 
and 
. Inequality(11) then holds and reads as
As discussed previously, the mutual LODFs Lo,k and Lk,o are both relatively large in this case. This is indicative of strong collective interactions as illustrated in the following sections. The purely collective effects not only overshadow the individual outages' effects on link l, indeed they reverse the total flow over the link.
3.Collective effects in complex networks
As shown above, the impact of a double link failure is not given by the simple sum of the individual outages' effects, but strong collective effects may be present. Based on the intuition developed in the last section, we will introduce a quantifier in the following section that may be used to identify in which situations collective effects need to be taken into account and in which situations they may be neglected, thus being able to rely on results obtained for single link failures. We will test this predictor on different test grids, mainly on the ones shown in figure 4; the Scandinavian power grid extracted from the software package PyPSA-eur [44] (panel a) and the IEEE test case 118 [45] (panel b).
3.1.Quantifying the strength of collective effects
To understand the purely collective effects of a simultaneous outage of two given links o and k, we first calculate the difference between the real flow changes in case of an outage of both links 
and the naive prediction in terms of the sum of individual flow changes 
. The difference calculated according to equations (8) and (9) reads as
The overall prefactor 
is one if the product 
is zero and tends to infinity as the product approaches one. In order to write this expression more compactly, we introduce the matrix 
which has the row vectors
This matrix includes the topological properties of the rerouting problem and ignores the initial flows 
and 
, which are determined by the specific power injections that may be time-dependent. The approach thus allows to quantify the impact of collective effects purely based on the network topology. However, this comes at the price of potentially missing situations with very unusual flow patterns in which the approach presented here might not be valid any more to predict collective effects. To get an overall measure of the purely collective part of the failure of two specific links o and k, we define a single collectivity parameter ξ(o, k) by taking the ℓ2-norm 
of the matrix Ξ,
The collectivity parameter thus quantifies the overall collective effects in the whole network to be expected if links o and k fail and indicates whether links o and k interact strongly or not. Since interpreting the collectivity parameter ξ in this form is rather cumbersome and we are looking for an easily accessible criterion that tells us which pairs of links interact strongly, we further reduce this expression by making a few approximations. Since the LODFs are bounded by one, 
for all links 
, and are typically much smaller than one, in the order of 
, we expect terms of third order in LODFs to be negligible against terms of second order such that we can on average neglect the former ones. In doing so, we arrive at the following approximation for the collectivity parameter,
Furthermore, we can make the following observation allowing us to further simplify the expression: summing over all links l in a large network, Ll,o and Ll,k will vary a lot and may thus essentially be treated as random variables. We therefore expect the collectivity parameter ξ(o, k) to be predicted by the two non-varying quantities Lo,k and Lk,o characterising the interaction between the two failing links. Since LODFs are in general non-symmetric (see appendix A), both Lo,k and Lk,o need to be incorporated to successfully predict the collectivity parameter ξ(o, k). In addition to that, we expect the prefactor 
to be well approximated by one in general, 
since (absolute) LODFs are typically small.
Based on these considerations, we introduce a parameter that predicts the overall strength of collective effects ξ and is defined as follows,
This predictor takes into account the relative effect of the failing links o and k on one another: it is the geometric mean of the mutual LODFs between the two failing links. It is not only a good predictor for the collectivity parameter ξ(o, k), but can also be shown to bound it from below as summarised in the following theorem.
Theorem 1.Consider a connected network where two links 
and 
with non-vanishing mutual LODFs 
fail. Then the collectivity parameter 
) as defined in equation (12) is bounded from below by the predictor 
A proof is given in appendix C. Figure 5 illustrates the performance of the predictor in forecasting collective effects for the IEEE test case 118 (panel a) and the Scandinavian power grid (panel b) when averaging over all possible trigger links. The predictor Λ(o, k) (abscissa) has a Pearson correlation coefficient with the collectivity parameter ξ (ordinate) of ρ=0.998 for both grids, thus indicating a linear relationship between the two quantities.
The predictor performs equally well if we replace the Euclidean ℓ2-norm in the definition of the collectivity parameter ξ in equation (12) by other ℓp-norms Norms with p>2 are dominated more by large values compared to norms with p<2, therefore we also tested the p=10-norm and even up to the 
-norm, which simply takes the maximum value. The predictor performs very well in forecasting collective effects also for other test grids and norms as summarised in table 1. For all norms and all grids tested, we observe a very strong correlation between predictor and collectivity parameter, exceeding ρ=0.9 in most cases. We discuss the predictor and the different norms used to calculate it in more detail in appendix B.
Table 1.Pearson correlation ρ between predictor Λ(o, k) and collectivity parameter ξ(o, k) in the case of a double outage of links o and k for all possible pairs of inks o and k and different test grids. Values are given for a number of test grids, namely IEEE 'case30', 'case118' and 'pegase1354' [45, 46] as well as the Scandinavian grid, the German grid, [47, 48] a periodic square grid with 20×20 nodes and another one with a share of s=0.45 of its links removed, and different norms used to calculate the collectivity parameter ξ(o, k). While ξ(o, k) is predicted very well for all norms, the 2-norm consistently yields the best results, albeit by a small margin.
| Pearson correlation ρ Λ(o, k) versus ξ(o, k) | ||||
|---|---|---|---|---|
| test grid | 1-norm | 2-norm | 10-norm |  -norm |
| case30 | 0.959 | 0.98 | 0.951 | 0.946 |
| case118 | 0.947 | 0.984 | 0.972 | 0.97 |
| pegase1354 | 0.933 | 0.974 | 0.967 | 0.966 |
| Square grid | 0.994 | 0.999 | 0.995 | 0.992 |
| Sparse square grid | 0.946 | 0.968 | 0.981 | 0.976 |
| Scandinavia | 0.909 | 0.964 | 0.968 | 0.967 |
| Germany | 0.869 | 0.951 | 0.927 | 0.921 |
To summarise, we find that two links show the strongest collective interaction if their mutual LODF values are large, thus implying that a failure of one link has a strong effect on the flow going over the other one and vice versa.
3.2.Impact of network distance
Distance is known to play an important role for failure spreading in power grids and other types of linear flow networks [22, 49, 37, 50]. In this section, we will examine if it may also be used to successfully predict collective effects in multiple link failures. Typically, distances in networks are measured between two nodes with the most prominent distance measure being the geodesic distance. It is given by the sum of the lengths or weights of all edges along a shortest path between the respective nodes,
where the superscript 'u' or 'w' denotes the unweighted or weighted distance, the subscript '0' describes the distance in the initial graph before any kind of outage, v1 and v2 are the nodes whose distance is calculated, 
is a path from v1 to v2 and ℓe is the length or weight of edge e, which is set to unity when calculating unweighted distances. For our purposes, we use the inverse link strength as the length, le=1/be. Additionally, one can define the geodesic distance between edges as the smallest possible distance between the nodes incident to the corresponding edges plus half of each edge's length,
Here the subscript 'ge' denotes the geodesic distance while (r, s) and (m, n) are the respective edges given by the nodes they are incident in. As we demonstrated in a recent publication [22], this distance measure does not capture essential aspects of the flow rerouting after a link failure. Instead, we proposed the rerouting distance
given by the length of the shortest cycle crossing both edges (r, s) and (m, n). If no such cycle exists, the rerouting distance is defined to be 
. This distance measure is strongly correlated with the magnitude of the LODFs as shown in [22].
Figure 6 shows the scaling of ξ with distance between the failing links for a failure of two links. Here, we use the Kendall rank correlation τ to quantify the degree of general non-linear correlation between the two quantities. The rerouting distance performs slightly better in predicting the collectivity parameter ξ than the geodesic distance, where the former has a rank correlation of τ=−0.7 and τ=−0.62 with ξ and the latter a correlation of τ=−0.66 and τ=−0.51 in the test grid 'case118' and the Scandinavian grid, respectively. Thus both distance measures perform moderately in predicting the collectivity parameter although not nearly as well as the predictor introduced in the last section. Still, the distance seems to be an important factor in determining the simultaneous outages' effects–but contrary to the case of a single outage [22], other factors play an important role, too. We may thus deduce that links that are closer to each other in both, the rerouting distance and the simple edge distance tend to have a stronger collective response. This behaviour is expected given that the predictor performing best is given by the product of the mutual LODFs between the two failing links and the rerouting distance is known to be a good predictor for the LODF [22].
3.3.Impact of community structures
Besides the effect of distance on failure spreading analysed in the last section, community structure is also known to play an important role. We analyse this effect exemplarily for the Scandinavian grid whose community structure we present in figure 7(a). The communities are determined using spectral graph partitioning [51–53] based on the sign of the eigenvector 
of the Laplacian matrix corresponding to the second smallest eigenvalue 
, also referred to as Fiedler eigenvalue [31]. The affiliation of a node to one of the two communities is encoded in red and blue colouring representing the value of the eigenvector at the respective node, such that two nodes belong to the same community if they have the same colour. As shown in figure 7(b), we evaluate ξ for the case were both trigger links are located in the same community (right, green boxplots) and different communities (left, blue boxplots). Overall, we observe that community structure plays an important role for the overall strength of collective effects; if both links, are located in the same community, the collectivity parameter is on average three orders of magnitude larger evaluated in terms of the respective median.
Community structure affects ξ directly and indirectly through the distance, because the distance between two links tends to be higher if they are in different communities. To account for this effect, we also present the collectivity parameter ξ in dependence of the rerouting distance (c, cf equation (14)) and the geodesic distance (d). We consistently observe a stronger collectivity parameter for both distance measures when the trigger links are contained in the same communitiy (green, upper curve) compared to the case of trigger links contained in different communities (blue, lower curve). We thus conclude that community structure also determines if two links interact strongly collectively.
4.Extension to arbitrary link failures
Now that we analysed the simultaneous failure of two links in detail, we will extend the theoretical framework to more than two links failing. To this end, we will derive a formula that describes this type of contingencies on a nodal level and perform a continuum limit that is valid for infinitely large regular grids.
4.1.Derivation of generalized LODFs
Now consider the simultaneous outage of K links 
with 
. Then we define the projection matrix from the space of all links onto the subset of failing links 
via
where 
denotes the Kronecker delta. Consequently, let 
be the projection of the node-edge-incidence matrix 
onto the subset of failing links which reads as
With this definition the columns of this matrix are the vectors 
introduced in the definition of the LODF, as per equation (6). As a reminder, they are defined by their entries being +1 at the node corresponding to the start of the respective failing link, −1 at the node corresponding to the end of the failing link and 0 otherwise. Furthermore, we define the projected branch reactance matrix 
by
Using this matrix, we can also project the vector of all initial flows onto the failing links defined by
The failure of multiple links may be regarded as a perturbation to the graph Laplacian 
in the same way as for a single link, see [22],
where 
is the graph Laplacian after the failure of the K links. The corresponding perturbation matrix 
may then also be expressed using the projected node-edge-incidence matrix as
In addition to that, the failure causes the nodal potentials to change,
where 
is a vector of the changes in angles. Using the continuity equation (3) in the new grid,
subtracting from it the current balance for the old grid, and applying the Moore–Penrose-pseudoinverse to the resulting equation, the change in potential is calculated as
We can simplify this expression by making use of the Woodbury matrix identity [55] and arrive at the final result,
Here, we defined a projection of the PTDF matrix onto the subset of failing links 
given by
The change in phase angles may then be used to calculate the flow changes by making use of equation (1). The vector of flow changes reads as
In principle, we may now make use of equation (18) to calculate the flow changes after an arbitrary number of simultaneous contingencies. The immediate insight into the structure of the contingency problem from this equation is, however, limited. We will thus try to gain more insight into the interplay of multiple outages by rearranging the equation.
Starting with equation (17) expressing the change of voltage phase angles after the failure 
, we can derive the following Poisson-like equation similar to the case of a single link failure as presented in equation (7)
where we defined the vector of flows weighted by the dipole source terms
We may thus rewrite equation (19) for the change in nodal potentials as follows, making the correspondence to the Poisson equation more apparent;
with the dipole sources
and the 
being the rows of 
, see also section 2.2. In addition to this expression, we can derive an analogous equation for the graph 
from which all the failing links have been removed. Simply plugging equation (15) into equation (16), we arrive at the following equation
We are thus left with a discrete Poisson equation with potential 
, which is analogous to the result obtained in our previous work [22], sections3 and 4 for a single failing link. Instead of a single dipole source this equation is governed by K dipole sources. However, this equation differs from the naive approach obtained by simply superposing single dipole sources. To see this, consider the case of K=2 link failures. As we have seen in the previous section 3, collective effects play an important role in the interaction of the two links. In this case, a simple superposition of two dipoles results in the equation
where we defined the dipole sources resulting from the naive approach
On the other hand, using the exact approach in equation (20), the actual dipole sources read as
The naive approach thus underestimates the interaction between the two dipole sources encoded in the matrix inverse 
, as discussed in detail in section 3. In the following paragraph, we will demonstrate, however, that this collective effect can be neglected in the continuum limit, thus making the naive superposition approach exact in that case.
4.2.Continuum limit for regular square lattice
We will now demonstrate how one may derive an exact formula for the potential changes after an arbitrary number of link failures for the setup of an infinite square lattice extending on our previous work [22]. Consider the elementary example of a regular square lattice embedded in the plane 
. Label all nodes by their positions 
and let the lattice spacing be denoted by h. Now introduce continuous functions ψ and b such that ψ(x, y) is the potential of the node at (x, y) and 
is the weight of the link connecting the two nodes at (x, y) and 
and analogously for two nodes connected in y-direction.
For a small lattice spacing 
and an infinitely large grid, the left-hand side of the Poisson equation (20) evaluated at position (x, y) can be written in a continuum version as [22]
Then, the flow changes according to equation (18) are given by
where 
refers to the change in flow due to the link failures here and should not be confused with the continuous Laplace operator.
The right-hand side of the Poisson equation (20) may be calculated similarly noting that at most 2K nodes contribute when K links fail. Let any failing link 
connect the nodes sk and tk with positions 
and 
respectively. The discrete version of a single addend on the right-hand side reads:
We will now show how this equation may be interpreted in the continuum version. First, the flow on a failing link before the outage 
may be calculated as
where 
is the continuum version of the flow before the outage. Second, the vector 
can be formally interpreted in terms of the two-dimensional delta function 
and reads for a link lk oriented in x-direction
For links oriented in y-direction, we simply replace 
by 
. Finally, in order to calculate the continuum version of the inverse matrix elements 
for two arbitrarily chosen links lk and li, assume without loss of generality that both links are oriented in x-direction and that a continuum version 
of the Green's function 
exists. Then, the elements of the projected PTDF matrix may be calculated as
All off-diagonal entries are zero due to the delta functions' different arguments. Importantly, this observation is independent of the orientation of the two links under consideration. The inverted matrix is thus diagonal and can be calculated as
In total, we obtain after expanding the entire expression to lowest order in the continuum limit
We can now formally divide the left-hand side(21) and the right-hand side(22) by h2 and take the limit 
to obtain the final continuum limit of the Poisson equation,
where the source terms are 
, the unperturbed current field.
If the link weights are hom*ogeneous, b(x, y)=b, the solution is given by the superposition of K two-dimensional dipole fields
We thus obtain a fully analytic solution in the continuum limit. This solution reveals that in hom*ogeneous lattices the effects of multiple outages are given by the superposition of single outages.
5.Conclusion and outlook
In this article, we have shown that multiple link failures can lead to fundamentally different impacts than expected from a naive superposition of single link failures. We have also established a parameter, namely the predictor Λ(o, k), quantifying in which cases these effects have to be taken into account. In addition to that, we have extended on previous work demonstrating that multiple link failures correspond to the overlay of correspondingly many single dipoles in infinitely large regular grids, thus allowing for a description similar to single link failures in this case. However, the strength of the effective dipoles is strongly determined by the collective effects, i.e. the interplay of the failing links. Our results demonstrate that further understanding of multiple link failures is an important task for the development and security of future power systems, thus helping to understand in which cases additional link shutdowns can help or counteract overall system security.
We have presented several elementary examples which demonstrate the counterintuitive behaviour of collective effects in some particular cases. Furthermore, we have shown that additional outages can be beneficial for the overall grid loading, thus presenting another occurence of Braess' paradox in power grids. In addition to that, we have shown that collective effects may lead to a sign inversion of flow direction compared to the individual failure of each single link. Both phenomena are potentially of high relevance when operating power grids as they might help to resolve situations where a single link fails or a redispatch occurs. However, further work should be dedicated to understanding and predicting these particular collective effects on a more fundamental theoretical basis.
The predictor for collective flow changes introduced in this manuscript allows for an easier understanding of when collective effects become considerably important. Mostly, collective effects are small, they only become relevant in cases where both failing links have a strong effect on one another. This is for example the case if the links are in close proximity or if they are both bottlenecks. Conversely, this implies that the intuition developed for single link failures may in many cases also be applied to study multiple link failures if the possibility of collective effects is kept in mind.
Distance between the two failing links seems to play an important role for the overall collective effects. Previous work has addressed the role of distance in flow changes for single link failures where in particular the rerouting distance was shown to be a decisive measure in predicting the flow changes [22]. In predicting collective effects, the rerouting distance between the failing links still seems to be an important quantity but not to the same degree as it is important for single link failures. In the future, it would be interesting to extend the rerouting distance to more than two links which could potentially also predict collective effects better.
Further work should address the role particular flow patterns play in more detail. In the approach used here, we abstracted from individual flow patterns and focused on topological aspects. This should be a good approximation for many cases, in particular when dealing with large power grids. However, there might be cases in which specific lines are nearly always more heavily loaded than other ones which should imply a more important contribution of such lines to collective effects.
Acknowledgments
We gratefully acknowledge support from the German Federal Ministry of Education and Research (BMBF grant no. 03EK3055B) and the Helmholtz Association (via the joint initiative 'Energy System 2050—a contribution of the research field energy' and the grant no. VH-NG-1025 to D W).
Appendix A.: Symmetry of LODFs
The LODFs according to equation (6) are given by
In this section, we will study the symmetry of this matrix in terms of interchanging the role of the failing link k and the link l on which flow changes are monitored. This symmetry describes the extend to which the flow change on one link l due to another, failing link k corresponds to the opposite flow change on link k if link l fails instead and thus provides a measure of symmetry for the whole network. In particular, we analyse how the matrix becomes asymmetric with an increasing degree of asymmetry in the links surrounding the monitored, and failing link. This explains why both LODFs are important for predicting the strength of collective effects in the predictor(13).
If we assume hom*ogeneous edge weights for the time being such that 
, we notice that the numerator in equation (6) is symmetric with respect to interchanging l and k. This numerator is also referred to as PTDF in power engineering [11]. The symmetry can be seen by taking the matrix transpose of the expression
On the other hand, the denominator is non-symmetric even in the case that line susceptances are uniform, as it reads as 
for Ll,k and 
for Lk,l thus encoding the importance of the link that fails. The LODFs are only completely symmetric if both links have not only the same weights, but also the same topological structure around them. This is for example the case for the periodic square grid, see figure A1, dark blue dots corresponding to s = 0 in the legend. We analysed this expression in detail in our previous publication [22] and showed that it can be predicted using the minimum cut that disconnects the two vertices k1 and k2. In the case of uniformly distributed line susceptances, this expression simply reduces to the well-studied resistance distance [56, 57]. Thus, the asymmetry in LODFs is encoded in both, the asymmetry in connectivity in the network, i.e. the variance in the node degree, and the asymmetry in the line susceptances.
In addition to the observation made before, we can notice that the mutual LODFs Ll,k and Lk,l will always have the same sign. This is due to the fact that as discussed above, the numerator is the same for both expressions. On the other hand, the denominator is always positive or equal to zero
Whereas this is not obvious from the above expression, it follows from the definition of the PTDFs given by 
[11]. Therefore we may conclude that mutual LODFs always have the same sign.
In figure A1, we demonstrate how asymmetry in LODFs arises with an increasing degree in inhom*ogeneity in the nodal degrees for a periodic square lattice from which we randomly remove a fraction s of its total number of links according to the procedure described in [22]. In this figure, we plot Ll,k against Lk,l for all possible combinations of links l and k. Starting at s=0 for a perfect periodic square lattice with 50 × 50 nodes, the LODF is perfectly symmetric (dark blue dots). With increasing degree of sparsity s∈{0.1, 0.2, 0.3, 0.4, 0.45} (dots from dark blue to light blue), we observe an increasing spread of the LODFs indicating an increasing degree of asymmetry in the LODFs.
Appendix B.: A predictor for collective effects
To support the choice of the predictor for collective effects, consider figures B1(a) and (b). We show the predictor Λ(o, k) and the remaining term, referred to as 'other term' in the figure, for all possible combinations of trigger links for the test case 'pegase1354' [46]. The remaining term is constructed by factoring out the predictor Λ(o, k) defined in equation (13) from the expression for the collectivity parameter ξ(o, k) given in equation (12) and assuming both LODFs Lk,o and Lo,k to be non-zero,
Applying the approximations discussed in section 3 to this equation, this expression reduces to the following equation,
Based on this expression, we discuss certain limiting cases which explain the performance of the predictor. Assume that Lo,k is very small keeping Lk,o constant and much larger than Lo,k. In this case, the expression is dominated by 
which is predicted well by Lk,o. Performing the same approximation for small values of Lk,o keeping Lo,k constant and large, the expression is well predicted by Lo,k. For this reason, we need to keep both values in order to predict the overall collective effects. On the other hand, one can easily check that the approximation is equally valid if both LODFs are of the same order.
Importantly, considering an arbitrary ℓp-norm instead, the conclusions differ only slightly. An ℓp norm 
of an arbitrary vector 
is defined as
For the predictor, we then have
Therefore, the general expression to be considered does not change fundamentally when calculating the ℓp norm instead.
Appendix C.: Proof of theorem 1
Proof.The collectivity parameter defined by equation (12) reads as
We will demonstrate that ξ(o, k) is bounded from below by Λ(o, k). First, since all addends in the sum are greater than zero, neglecting any or all of them will not increase the expression's value. We can thus choose the addends with l=o and l=k. The expression then reads as
Now we can make use of the fact that Lo,o=Lk,k=−1. In order to show that this expression is bounded from below by Λ(o, k), we can square both sides of the inequality since all expressions considered here are positive. This yields
The last inequality follows from the following considerations
which completes the proof. For 
, the inequality is strict.□
Note that the proof makes use of the fact that Lo,o=−1 for all links o, but we expect the statement to hold even without this assumption.
