Non-local impact of link failures in linear flow networks

The robust operation of supply networks is essential for the function of complex systems across scales and disciplines. Almost all of our technical and economical infrastructure depends on the reliable operation of the electric power grid [1, 2]. Biological organisms distribute water and nutrients via their vascular networks, for instance in plant leaves [3], the human and animal circulatory system [4], or in protoplasmic veins of certain slime molds [5]. Huge amounts of money and assets are exchanged through a complex financial network [6]. Structural damages to such networks can have catastrophic consequences such as a stroke, a power outage or a major economic crisis.

In power grids, large scale outages are typically triggered by the failure of a single transmission or generation element [7–11]. The outages in the United States in 2003, Italy in 2003 and Western Europe in 2006 are very well documented and provide a detailed insight into the dynamics of a large scale network failure [12–14]. Each outage was triggered by the loss of a transmission line during a period of high grid load. Subsequently, the power flows were rerouted, causing secondary overloads and eventually a cascade of failures. In these three examples, the cascades propagated mostly locally—overloads took place in the proximity of previous failures. However, this is not necessarily the case during power outages (see, e.g. [15]), raising the question of how networks flows are rerouted after failures [16–23].

In biological distribution networks, robustness against link failure is a critical prerequisite that guards against potentially life-threatening events such as stroke [24] or embolism [25, 26], but also to function efficiently in the presence of fluctuations [3, 27, 28]. Thus, biological networks are often (but not in all cases, such as in the penetrating arterioles of the cortical vasculature [29]) endowed with highly resilient, redundant topologies that optimize rerouting of flow in case of link failure to the network [28] and are generated through adaptive developmental mechanisms [30]. For the understanding of such life-threatening conditions it is therefore crucial to investigate the behavior of the vascular network in the case of failure.

To understand the vulnerability of networks, we here provide a detailed analysis of the impact of link failures in linear flow networks. We focus on how the network topology determines the overall network response as well as the spatial flow rerouting. We consider linear supply network models, where the flow between two adjacent nodes is proportional to the difference of the nodal potential, pressure or voltage phase angle. Linear models are applied to hydraulic networks [31], vascular networks of plants and animals [28, 32–35], economic input–output networks [36] as well as electric power grids [37–42]. The linearity allows to obtain several rigorous bounds for flow rerouting in general network topologies and to solve special cases fully analytically.

The paper is structured as follows; first, we formally introduce linear flow networks in section 2 and present a framework for studying line outages in such networks in section 3. Afterwards, we establish a mathematical analogy of flow rerouting after line outages and electric dipole fields on square grids in the continuum limit in section 4. We then derive rigorous bounds on the strength of this effective dipole to describe how the flow is rerouted on arbitrary network topologies in section 5. Finally, we establish a new distance measure on networks, the rerouting distance, which is able to predict the flow redistribution much better than the ordinary geodesic distance in section 6. Furthermore, we study the effect of network sparsity on the dipole pattern of flow redistribution and quantify this scaling with distance from the failing link in the same section.

Consider a network consisting of N nodes (vertices) that are connected to each other via lines (edges) denoted by $(m,n)$ for a line going from node m to node n. We assume the network to be globally connected, otherwise consider each connected component of the network separately. Assign a potential or phase angle ${\theta }_{m}\in {\mathbb{R}}$ to each node m in the network. Then we assume the flow ${F}_{m\to n}$ between nodes m and n connected via line $(m,n)$ to be linear in the potential drop along the line

$$\begin{eqnarray}&&{F}_{m\to n}={b}_{{mn}}({\theta }_{m}-{\theta }_{n}).\end{eqnarray} \tag{ 1 }$$

Here, ${b}_{{mn}}={b}_{{nm}}$ is the transmission capacity assigned to the line $(m,n)$ that describes its ability to carry flow. This equation may for example be used to describe hydraulic networks [31, 43] or vascular networks of plants [28], where the ${\theta }_{n}$ denotes the pressure at some node n and the capacity b_mn scales with the diameter of a pipe or vein. Our main focus will be its application to electric power engineering, where this linear approximation of the power flow equations is referred to as the DC approximation [38–40]. In this case, ${F}_{m\to n}$ refers to the flow of real power along a transmission line $(m,n)$, ${\theta }_{n}$ is the voltage phase angle at node n and b_mn is proportional to the line's susceptance. For the sake of consistency, we refer to the ${\theta }_{n}$ as 'potentials' throughout this paper. Since only phase differences are involved in the flow calculation, these potentials are only defined up to a constant phase shift. Typically, an arbitray node is selected where the potential is set to zero, ${\theta }_{n}=0$.

In addition to that, we assume that Kirchhoff's current law holds at the nodes of the network which states that the inflows and outflows at any node m balance

$$\begin{eqnarray}&&\displaystyle \sum _{n=1}^{N}{F}_{m\to n}={P}_{m},\end{eqnarray} \tag{ 2 }$$

where the right-hand side denotes the inflow (${P}_{m}\gt 0$) or outflow (${P}_{m}\lt 0$) at node m, commonly called the 'power injection' in power engineering. Equations (1) and (2) describe the state and the flow of the network up to a constant phase shift as described above once the line parameters b_mn and the injections P_m are given.

These equations may be conveniently written using a vectorial notation. Define the vector ${\boldsymbol{\theta }}={\left({\theta }_{1},\ldots ,{\theta }_{N}\right)}^{\top }\in {{\mathbb{R}}}^{N}$ of the nodal potentials or voltage phase angles and the vector ${\boldsymbol{P}}={\left({P}_{1},\ldots ,{P}_{N}\right)}^{\top }\in {{\mathbb{R}}}^{N}$ of nodal injections. Here and in the following sections, the superscript '$\top$' denotes the transpose of a vector or matrix. We further label all lines in the grid by ${\ell }=1,\ldots ,\ L$ and summarize all line flows in a vector ${\boldsymbol{F}}={\left({F}_{1},\ldots ,{F}_{L}\right)}^{\top }\in {{\mathbb{R}}}^{L}$. Equation (1) may then be rewritten as

$$\begin{eqnarray}&&{\boldsymbol{F}}={{\boldsymbol{B}}}_{d}{{\boldsymbol{K}}}^{\top }{\boldsymbol{\theta }},\end{eqnarray} \tag{ 3 }$$

where ${{\boldsymbol{B}}}_{d}\in {{\mathbb{R}}}^{L\times L}$ is a diagonal matrix containing the capacities ${b}_{{\ell }}$ of all edges. Furthermore, we defined the node-edge incidence matrix ${\boldsymbol{K}}\in {{\mathbb{R}}}^{N\times L}$. To define this matrix in an undirected graph, one typically fixes an arbitrary orientation of the graph's edges such that its components read

$$\begin{eqnarray}\begin{array}{r}{K}_{n,\ell }=\left\{\begin{array}{cl}1 & {\rm{i}}{\rm{f}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{e}}\,\ell \,{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{s}}\,{\rm{a}}{\rm{t}}\,{\rm{n}}{\rm{o}}{\rm{d}}{\rm{e}}\,n,\\ -1 & {\rm{i}}{\rm{f}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{e}}\,\ell \,{\rm{e}}{\rm{n}}{\rm{d}}{\rm{s}}\,{\rm{a}}{\rm{t}}\,{\rm{n}}{\rm{o}}{\rm{d}}{\rm{e}}\,n,\\ 0 & {\rm{o}}{\rm{t}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{w}}{\rm{i}}{\rm{s}}{\rm{e}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

The node-edge incidence matrix also relates the injections to the flows incident at a node. More specifically, Kirchhoff's current law (2) may be rewritten as follows

$$\begin{eqnarray}&&{\boldsymbol{P}}={\boldsymbol{K}}{\boldsymbol{F}}={\boldsymbol{K}}{{\boldsymbol{B}}}_{d}{{\boldsymbol{K}}}^{\top }{\boldsymbol{\theta }}={\boldsymbol{B}}{\boldsymbol{\theta }}.\end{eqnarray} \tag{ 5 }$$

Here, we defined the matrix ${\boldsymbol{B}}={\boldsymbol{K}}{{\boldsymbol{B}}}_{d}{{\boldsymbol{K}}}^{\top }\in {{\mathbb{R}}}^{N\times N}$ commonly referred to as the nodal susceptance matrix in power engineering. Mathematically, ${\boldsymbol{B}}$ is a weighted Laplacian matrix [44, 45] with components

$$\begin{eqnarray*}\begin{array}{r}{B}_{mn}=\left\{\begin{array}{ll}{\displaystyle \sum }_{\ell \in {{\rm{\Lambda }}}_{m}}{b}_{\ell } & {\rm{i}}{\rm{f}}\,m=n;\\ -{b}_{\ell } & \,{\rm{i}}{\rm{f}}\,m\,\,{\rm{i}}{\rm{s}}\,{\rm{c}}{\rm{o}}{\rm{n}}{\rm{n}}{\rm{e}}{\rm{c}}{\rm{t}}{\rm{e}}{\rm{d}}\,{\rm{t}}{\rm{o}}\,n\,{\rm{b}}{\rm{y}}\,\ell .\end{array}\right.\end{array}\end{eqnarray*}$$

Here, ${{\rm{\Lambda }}}_{m}$ is the set of lines which are incident to m.

An important question in network security analysis is how the flows in the network change if a line fails. Denoting by F_ℓ the initial flow of the failing line ${\ell }\hat{=}(r,s)$, the flow change ${\rm{\Delta }}{F}_{e}$ at a transmission line $e\hat{=}(m,n)$ is written as

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{e}={{\rm{LODF}}}_{e,{\ell }}\,{F}_{{\ell }}.\end{eqnarray} \tag{ 6 }$$

Adopting the language of power system security analysis [37, 38], we call the factor of proportionality the line outage distribution factor (LODF). In the following, we present two alternative derivations as well as a physical interpretation of the linear flow rerouting problem.

3.1. Self-consistent derivation of LODFs

To derive an explicit expression for the LODFs one generally starts with a related problem. Consider an increase of the real power injection at node r and a corresponding decrease at node s by the amount ${\rm{\Delta }}P$. The new vector of real power injections is then given by

$$\begin{eqnarray}&&\hat{{\boldsymbol{P}}}={\boldsymbol{P}}+{\rm{\Delta }}P\,{{\boldsymbol{\nu }}}_{{rs}},\end{eqnarray} \tag{ 7 }$$

where the components of ${{\boldsymbol{\nu }}}_{{rs}}\in {{\mathbb{R}}}^{N}$ are +1 at position r, −1 at position s and zero otherwise. Here and in the following, we use a hat to indicate the state of the network after a line outage or a similar change of the topology. The change of the real power injections causes the following change in the real power flow

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{mn}=\mathop{\underbrace{{b}_{mn}{{\boldsymbol{\nu }}}_{mn}^{{\rm{\top }}}{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{rs}}}\limits_{=:{{\rm{P}}{\rm{T}}{\rm{D}}{\rm{F}}}_{(m,n),r,s}}{\rm{\Delta }}P.\end{eqnarray} \tag{ 8 }$$

Here, ${{\boldsymbol{B}}}^{\dagger }$ denotes the Moore–Penrose pseudo-inverse of the Laplacian matrix ${\boldsymbol{B}}$ and the factor of proportionality is referred to as the power transfer distribution factor (PTDF).

The LODFs can be expressed by PTDFs in the following way [38]. To consistently model the outage of line $(r,s)$, one assumes that the line is disconnected from the grid by circuit breakers and that some fictitious real power ${\rm{\Delta }}P$ is injected at node s and taken out at node r. The entire flow over the line $(r,s)$ after the opening thus equals the fictitious injections ${\hat{F}}_{{rs}}={\rm{\Delta }}P$. Using PTDFs, we also know that

$$\begin{eqnarray*}&&{\hat{F}}_{{rs}}={F}_{{rs}}+{{\rm{PTDF}}}_{(r,s),r,s}\,{\rm{\Delta }}P.\end{eqnarray*}$$

Substituting ${\hat{F}}_{{rs}}={\rm{\Delta }}P$, solving for ${\rm{\Delta }}P$ and inserting equation (8) yields

$$\begin{eqnarray}&&{{\rm{LODF}}}_{({mn}),({rs})}=\displaystyle \frac{{{\rm{PTDF}}}_{(m,n),r,s}}{1-{{\rm{PTDF}}}_{(r,s),r,s}}.\end{eqnarray} \tag{ 9 }$$

For consistency, one usually defines the ${\rm{LODF}}$ for the failing line as follows: ${{\rm{LODF}}}_{({rs}),({rs})}=-1$. In addition to that, we exclude cases where the failing line is a bridge, i.e. a line whose removal disconnects the graph, from our analysis in the following sections.

3.2. Algebraic derivation of LODFs

The LODFs can also be obtained in a purely algebraic way without any self-consistency argument [46]. As the line $\ell \hat{=}(r,s)$ fails, the nodal susceptance matrix of the network changes as

$$\begin{eqnarray}&&{\boldsymbol{B}}\to \hat{{\boldsymbol{B}}}={\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}},{\rm{where}}\,{\rm{\Delta }}{\boldsymbol{B}}={B}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top },\end{eqnarray} \tag{ 10 }$$

which causes a change of the nodal potentials or voltage phase angles respectively,

$$\begin{eqnarray}&&{\boldsymbol{\theta }}\to \hat{{\boldsymbol{\theta }}}={\boldsymbol{\theta }}+{\boldsymbol{\psi }}.\end{eqnarray} \tag{ 11 }$$

Equation (5) for the perturbed grid now reads

$$\begin{eqnarray}&&({\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}})({\boldsymbol{\theta }}+{\boldsymbol{\psi }})={\boldsymbol{P}}.\end{eqnarray} \tag{ 12 }$$

Subtracting equation (5) for the unperturbed grid, we see that the change of the voltage angles is given by

$$\begin{eqnarray}&&{\boldsymbol{\psi }}=-{\left({\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}}\right)}^{\dagger }{\rm{\Delta }}{\boldsymbol{B}}\,{\boldsymbol{\theta }}={\left({\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}}\right)}^{\dagger }\,{{\boldsymbol{\nu }}}_{{rs}}{F}_{{rs}}.\end{eqnarray} \tag{ 13 }$$

The change of flows after the outage of line $(r,s)$ and thus the LODFs are calculated from the change of the voltage angles which yields

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{F}_{{mn}} & = & {b}_{{mn}}({\psi }_{m}-{\psi }_{n})={b}_{{mn}}{{\boldsymbol{\nu }}}_{{mn}}^{\top }{\boldsymbol{\psi }}\\ & = & {b}_{{mn}}{{\boldsymbol{\nu }}}_{{mn}}^{\top }{\left({\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}}\right)}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}{F}_{{rs}}.\end{array}\end{eqnarray} \tag{ 14 }$$

In principle, we could now use these equations to calculate the flow changes and the LODFs. However, this would require to invert the matrix $\hat{{\boldsymbol{B}}}={\boldsymbol{B}}+{\rm{\Delta }}{\boldsymbol{B}}$ separately for every possible line $(r,s)$ in the grid, which is impractical. Nevertheless, we can simplify the expression using the Woodbury matrix identity,

$$\begin{eqnarray*}&&{\left({\boldsymbol{B}}+{B}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }\right)}^{\dagger }={{\boldsymbol{B}}}^{\dagger }-{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}{\left({B}_{{rs}}^{\dagger }+{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{\nu }}}_{{rs}}\right)}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }.\end{eqnarray*}$$

Thus we obtain

$$\begin{eqnarray}&&{\left({\boldsymbol{B}}+{B}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }\right)}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}={\left(1+{B}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}\right)}^{-1}\,{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}},\end{eqnarray} \tag{ 15 }$$

such that the flow change (14) reads

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{{mn}}=\displaystyle \frac{{b}_{{mn}}{{\boldsymbol{\nu }}}_{{mn}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}}{1-{b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}}\,\times {F}_{{rs}}\,,\end{eqnarray} \tag{ 16 }$$

which is identical to equation (9) obtained using the standard approach.

3.3. Electrostatic interpretation

A deeper physical insight into the network flow rerouting problem is obtained by the analogy to discrete electrostatics. Equation (13) can be rearranged into a linear set of equations for the change of the nodal potentials

$$\begin{eqnarray}&&\hat{{\boldsymbol{B}}}{\boldsymbol{\psi }}={F}_{rs}{{\boldsymbol{\nu }}}_{rs}.\end{eqnarray} \tag{ 17 }$$

Here, $\hat{{\boldsymbol{B}}}$ is the Laplacian of the grid after the failure, i.e. the grid where line $(r,s)$ has been removed. Alternatively, we can formulate the equation in terms of the original network topology, substituting equation (15) into equation (13). This yields the linear set of equations

$$\begin{eqnarray}&&{\boldsymbol{B}}{\boldsymbol{\psi }}={\boldsymbol{q}}\end{eqnarray} \tag{ 18 }$$

with the dipole source

$$\begin{eqnarray}&&{\boldsymbol{q}}={\left(1-{b}_{rs}{{\boldsymbol{\nu }}}_{rs}^{{\rm{\top }}}{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{rs}\right)}^{-1}{F}_{rs}{{\boldsymbol{\nu }}}_{sr}.\end{eqnarray} \tag{ 19 }$$

As noted before, ${\boldsymbol{B}}$ and $\hat{{\boldsymbol{B}}}$ are Laplacian matrices and the right-hand side of both equations (18) and (17) are non-zero only at positions r and s with opposite sign. Hence, these equations are discrete Poisson equations with a dipole source and ${\boldsymbol{\psi }}$ is a dipole potential, see [47, ch 15] for a detailed analysis of this equation. The main complexity of the line outage problem thus arises from the network topology encoded in the Laplacian ${\boldsymbol{B}}$, which can be highly irregular.

The two equations (18) and (17) yield the same potential ${\boldsymbol{\psi }}$, but are formulated on different topologies—either on the original network topology or the topology after the outage. To compare the impact of different failures it is beneficial to use the original topology, such that only the dipole inhomogeneity differs—not the electrostatic problem itself. Then, the strength of the dipole depends on the network topology via the prefactor ${\left(1-{b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}\right)}^{-1}$.

Using the analogy to electrostatics we can solve the flow rerouting problem for regular network topologies (section 4) and provide some general rigorous results (section 6.1). To understand flow rerouting in networks with complex topologies, we thus have to account for the spatial spreading pattern described by ${{\boldsymbol{B}}}^{\dagger }$ (see section 6.3) as well as the dipole strength, which quantifies the gross response of the grid (see section 5).

To obtain a first insight into the spatial aspects of flow rerouting, we consider an elementary example admitting a solution in the continuum limit. Consider a regular square lattice embedded in a plane as depicted in figure 1 and studied in a slightly different form in [48]. All nodes are labeled by their positions ${\boldsymbol{r}}={\left(x,y\right)}^{\top }$ in this two-dimensional embedding and the lattice spacing is denoted as h. We introduce continuous functions ψ and b such that $\psi (x,y)$ is the potential of the node at $(x,y)$ and $b(x+h/2,y)$ is the weight of the link connecting the two nodes at $(x,y)$ and $(x+h,y)$. The left-hand side of the Poisson equation (18) evaluated at position $(x,y)$ reads

$$\begin{eqnarray}\begin{array}{rcl}({\boldsymbol{B}}\psi )(x,y) & = & b(x+h/2,y)\left[\psi (x,y)-\psi (x+h,y)\right]+b(x-h/2,y)\left[\psi (x,y)-\psi (x-h,y)\right]\\ & & +b(x,y+h/2)\left[\psi (x,y)-\psi (x,y+h)\right]+b(x,y-h/2)\left[\psi (x,y)-\psi (x,y-h)\right]\\ & = & -{h}^{2}{\rm{\nabla }}\cdot \left(b(x,y){\rm{\nabla }}\psi \right)+{ \mathcal O }({h}^{3}).\end{array}\end{eqnarray} \tag{ 20 }$$

Here, we made use of the fact that the components of the gradient ${\rm{\nabla }}\psi ={\left({\partial }_{x}\psi ,{\partial }_{y}\psi \right)}^{\top }$ may be expressed as

$$\begin{eqnarray*}&&\displaystyle \frac{\partial \psi (x,y)}{\partial x}=\mathop{\mathrm{lim}}\limits_{h\to 0}\displaystyle \frac{\psi (x+h,y)-\psi (x,y)}{h},\end{eqnarray*}$$

but did not take the limit yet. The derivative with respect to y may be calculated analogously.

Zoom In Zoom Out Reset image size

Before we proceed to the right-hand side, we remark that the flow changes ${\rm{\Delta }}{\boldsymbol{F}}$ according to equation (14) are given by

$$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{\rm{\Delta }}{F}_{mn} & \equiv & h{\rm{\Delta }}{{\boldsymbol{F}}}_{x}(x+h/2,y)=b(x+h/2,y)\left(\psi (x+h,y)-\psi (x,y)\right),\\ & & h{\rm{\Delta }}{{\boldsymbol{F}}}_{y}(x,y+h/2)=b(x,y+h/2)\left(\psi (x,y+h)-\psi (x,y)\right),\end{array}\end{array}\end{eqnarray*}$$

where $(m,n)$ denotes the link where the flow changes are monitored which is either oriented parallel to the x-axis, thus considering ${\rm{\Delta }}{{\boldsymbol{F}}}_{x}$ or the y-axis, thus considering ${\rm{\Delta }}{{\boldsymbol{F}}}_{y}$. If we divide by h and take the continuum limit $h\to 0$ the overall continuous flow changes read thus

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{F}}(x,y)=b(x,y){\rm{\nabla }}\psi (x,y).\end{eqnarray} \tag{ 21 }$$

Note that the expression ${\rm{\Delta }}{\boldsymbol{F}}$ refers to the change in flow due to the link failure here and should not be confused with the continuous Laplace operator.

The right-hand side of the discrete Poisson equation (18) may be calculated similarly noting that only two nodes contribute with opposite signs. Let us assume that the failing link is parallel to the x-axis connecting nodes r and s located at ${{\boldsymbol{r}}}_{r}={\left({x}_{r},{y}_{r}\right)}^{\top }$ and ${{\boldsymbol{r}}}_{s}={\left({x}_{s},{y}_{s}\right)}^{\top }={\left({x}_{r}+h,{y}_{r}\right)}^{\top }$. The discrete version of the right-hand side reads

$$\begin{eqnarray*}&&{\boldsymbol{q}}=\displaystyle \frac{{F}_{{rs}}}{1-{b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}}{{\boldsymbol{\nu }}}_{{rs}}.\end{eqnarray*}$$

We will now derive the continuum version of this equation. First, the flow on the failing link before the outage F_rs may be calculated as

$$\begin{eqnarray*}\begin{array}{rcl}{F}_{{rs}} & \equiv & b({x}_{r}+h/2,{y}_{r})(\theta ({x}_{r}+h,{y}_{r})-\theta ({x}_{r},{y}_{r}))\\ & = & {hb}({x}_{r},{y}_{r})\displaystyle \frac{\partial \theta }{\partial x}+{ \mathcal O }({h}^{2})\\ & = & h{{\boldsymbol{F}}}_{x}({x}_{r},{y}_{r})+{ \mathcal O }({h}^{2}),\end{array}\end{eqnarray*}$$

where ${\boldsymbol{F}}(x,y)=b(x,y){\boldsymbol{\nabla }}\theta (x,y)$ is the continuum flow before the outage. Second, the vector ${{\boldsymbol{\nu }}}_{{rs}}$ can be formally interpreted in terms of the two-dimensional delta function $\delta (x,y)$ and reads for the given link failure

$$\begin{eqnarray*}\begin{array}{rcl}{{\boldsymbol{\nu }}}_{{rs}} & \equiv & \delta (x-{x}_{r}+h,y-{y}_{r})-\delta (x-{x}_{r},y-{y}_{r})\\ & = & h\displaystyle \frac{\partial \delta (x-{x}_{r},y-{y}_{r})}{\partial x}+{ \mathcal O }({h}^{2}).\end{array}\end{eqnarray*}$$

Finally, let us assume that a continuum version of the Green's function ${{\boldsymbol{B}}}^{\dagger }$ exists. Then the denominator may be calculated as

$$\begin{eqnarray*}\begin{array}{rcl}{b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}} & \equiv & {h}^{2}b({x}_{r}+h/2,{y}_{r})\left(\displaystyle \int {\partial }_{y}\delta (x-{x}_{r},y-{y}_{r}){b}^{\dagger }(x,y){\partial }_{x}\delta (x-{x}_{r},y-{y}_{r}){\rm{d}}x\,{\rm{d}}y\right)\\ & = & {h}^{2}b({x}_{r},{y}_{r})\displaystyle \frac{{\partial }^{2}{b}^{\dagger }({x}_{r},{y}_{r})}{\partial x\partial y}+{ \mathcal O }({h}^{3}),\end{array}\end{eqnarray*}$$

where ${b}^{\dagger }(x,y)$ is the aforementioned continuum version.

Thus, in total we obtain after expanding the entire right-hand side to lowest order in the continuum limit

$$\begin{eqnarray}&&q(x,y)={h}^{2}{\boldsymbol{F}}{\left({x}_{r},{y}_{r}\right)}^{\top }{\boldsymbol{\nabla }}\delta (x-{x}_{r},y-{y}_{r})+{ \mathcal O }({h}^{3}).\end{eqnarray} \tag{ 22 }$$

Here, ${\boldsymbol{F}}({x}_{r},{y}_{r})$ is assumed to be parallel to the dipole axis, i.e. the direction of the link failure, which is either the x- or the y-direction for the given setting.

We can now formally divide left-hand side (20) and right-hand side (22) by h² and take the limit $h\to 0$ to obtain the final continuum limit of the Poisson equation,

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left(b(x,y){\rm{\nabla }}\psi \right)=-{{\boldsymbol{q}}}^{\top }{\rm{\nabla }}\delta (x-{x}_{r},y-{y}_{r}),\end{eqnarray} \tag{ 23 }$$

where the source term is ${\boldsymbol{q}}({x}_{r},{y}_{r})={\boldsymbol{F}}({x}_{r},{y}_{r})$, the unperturbed current field. We note that we obtain the same continuum limit regardless of whether we use equation (17) or (18) to do the expansions. Thus, the non-locality that is encoded in equation (18) is lost in the continuum formulation.

If the link weights are homogeneous, $b(x,y)=b$, and the failing link is assumed to be located at the origin $({x}_{r},{y}_{r})=(0,0)$ the solution is given by the well-known two-dimensional dipole field

$$\begin{eqnarray}&&\psi ({\boldsymbol{r}})=\displaystyle \frac{{\boldsymbol{q}}\cdot {\boldsymbol{r}}}{\parallel {\boldsymbol{r}}{\parallel }^{2}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{F}}({\boldsymbol{r}})=b\cdot \left(\displaystyle \frac{{\boldsymbol{q}}}{\parallel {\boldsymbol{r}}{\parallel }^{2}}-2{\boldsymbol{r}}\displaystyle \frac{{\boldsymbol{q}}\cdot {\boldsymbol{r}}}{\parallel {\boldsymbol{r}}{\parallel }^{4}}\right).\end{eqnarray} \tag{ 25 }$$

We thus obtain a fully analytic solution in the continuum limit. This solution reveals that the impact of link failures decays algebraically in homogeneous lattices. We consider this decay along two different axes. Assume the dipole to be located at the origin in x-direction, such that ${\boldsymbol{q}}={\left(\varepsilon ,0\right)}^{\top }$ where $\varepsilon \ll 1$ is some small real number. First, consider the decay in x-direction where ${\boldsymbol{r}}={\left(x,0\right)}^{\top }$. In this case, we obtain for the decay of the potential and the flow changes

$$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}\psi \left({\left(x,0\right)}^{{\rm{\top }}}\right) & = & \displaystyle \frac{\varepsilon \cdot x}{{x}^{2}}\propto \displaystyle \frac{1}{x},\\ {\rm{\Delta }}{\boldsymbol{F}}\left({\left(x,0\right)}^{{\rm{\top }}}\right) & = & b{\left(\displaystyle \frac{\varepsilon }{{x}^{2}}-2x\displaystyle \frac{\varepsilon x}{{x}^{4}},0\right)}^{{\rm{\top }}}=-b{\left(\displaystyle \frac{\varepsilon }{{x}^{2}},0\right)}^{{\rm{\top }}}.\end{array}\end{array}\end{eqnarray*}$$

This decay in the flow changes may also be observed in the discrete version of equation (22) and is shown in figure 2(a), for a line failure in a discrete square grid. Along the same lines, we may quantify the decay in y-direction where ${\boldsymbol{r}}={\left(\varepsilon ,y\right)}^{\top }$ for the same dipole orientation. In this case, we obtain

$$\begin{eqnarray*}\begin{array}{rcl}\psi \left({\left(\varepsilon ,y\right)}^{\top }\right) & = & \displaystyle \frac{{\varepsilon }^{2}}{{y}^{2}}\propto \displaystyle \frac{1}{{y}^{2}}\\ {\rm{\Delta }}{\boldsymbol{F}}\left({\left(\varepsilon ,y\right)}^{\top }\right) & \approx & b{\left(\displaystyle \frac{\varepsilon }{{y}^{2}},-2\displaystyle \frac{{\varepsilon }^{2}}{| y{| }^{3}}\right)}^{\top }.\end{array}\end{eqnarray*}$$

Here, we assumed the position vector to be dominated by its y-component, $\varepsilon \ll y$ such that $| | {\boldsymbol{r}}| | \approx | y|$. In total, we observe a ${y}^{-3}$-scaling in the flow changes in y-direction perpendicular to the dipole source and a ${y}^{-2}$-scaling in y-direction parallel to the dipole source, see figures 2(a)–(c).

Zoom In Zoom Out Reset image size

We now turn to realistic networks with irregular topologies. The change in the nodal potentials or voltage phase angles ${\psi }_{n}$ and flows ${\rm{\Delta }}{F}_{m\to n}$ is determined by the discrete Poisson equation (18). We first consider the right-hand side of this equation, the dipole strength, which describes the gross response of the network flows to the outage. This response is proportional to the initial flow of the failing edge F_rs and the factor

$$\begin{eqnarray}&&{\left(1-{b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}\right)}^{-1}=: {\left(1-\eta (r,s)\right)}^{-1}.\end{eqnarray} \tag{ 26 }$$

The factor $1-\eta (r,s)$ describes the non-locality of the network response to a local perturbation at link $(r,s)$. To see this, consider a grid where the real power ${\rm{\Delta }}P$ is injected and withdrawn at the terminal nodes of the link $(r,s)$. The direct flow over the link is given by

$$\begin{eqnarray}&&{F}_{r\to s}={b}_{{rs}}{{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}\,{\rm{\Delta }}P=\eta (r,s){\rm{\Delta }}P,\end{eqnarray} \tag{ 27 }$$

whereas the total flow is just given by ${\rm{\Delta }}P$. The factor $\eta (r,s)$ thus measures the fraction of the flow which is transmitted directly and $1-\eta (r,s)$ is the fraction transmitted non-locally via other pathways. Hence, $1-\eta (r,s)$ can also be seen as a measure of redundancy. A high non-local flow indicates that there are strong alternative routes from r to s in addition to the direct link $(r,s)$. If no alternative path exists, the flow must be routed completely via the direct link such that $1-\eta (r,s)=0$.

We conclude that the properties of alternative and direct paths are decisive for the understanding of flow rerouting. Before we proceed, we thus review the formal definition of a path in graph theory.

Definition 1. A path from vertex r to vertex s is defined as an ordered set of vertices

$$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},{v}_{2},\ldots ,{v}_{k}=s),\end{eqnarray} \tag{ 28 }$$

where two subsequent vertices must be connected by an edge and no vertex is visited twice. Two paths are called independent if they share no common edge. The unweighted length of such a path is defined as the number of steps k, while the weighted path length is given by the sum of the edge weights along the path, ${\sum }_{j=1}^{k}{w}_{{v}_{j-1}{v}_{j}}$. In this work, the edge weights are given by the inverse susceptances ${w}_{{ij}}=1/{b}_{{ij}}$. The (weighted or unweighted) geodesic or shortest path distance of two vertices r and s is defined as the (weighted or unweighted) length of the shortest path from r to s.

The interpretation as a redundancy measure directly relates the factor $1-\eta (r,s)$ to the topology of the network. A first rough estimate can be obtained from the topological connectivity ${\lambda }_{T}(r,s)$, which is defined as the number of independent paths from node r to node s. A comparison for several test grids in figure 3 shows that $\eta (r,s)$ decreases with ${\lambda }_{T}(r,s)$ on average as expected, but that there is a large heterogeneity between the links.

Zoom In Zoom Out Reset image size

To obtain a better topological estimate for the locality factor we need to take into account the heterogeneity of the link weights. The topological connectivity ${\lambda }_{T}(r,s)$ counts the minimum number of edges which have to be removed to disconnect the nodes r and s. We can define a weighted analog ${\lambda }_{F}(r,s)$ as the minimum capacity which has to be removed to disconnect the nodes r and s. This is a classical problem in graph theory, where it is referred to as the minimum cut [50]. We will now elaborate this quantity in a definition. An $(r,s)$-cut can be defined as follows. Let $r\in S\subset V$ and $s\in V\setminus S$ be two vertices taken from the two disjoint sets. The $(r,s)$-cut is defined as the set of edges $\delta (S)=\{(u,v)\in E\ | \ u\in S,v\in V\setminus S$ or $v\in S,u\in V\setminus S\}$ connecting the two disjoint vertex-sets. The set of edges ${\delta }^{+}(S)=\{(u,v)\in E\ | \ u\in S,v\in V\setminus S\}$ is referred to as the forward edges of the cut. The capacity C of a cut $\delta (S)$ and the corresponding minimum capacity ${\lambda }_{F}(r,s)$ between r and s are then given by

$$\begin{eqnarray*}\begin{array}{rcl}C(\delta (S)) & = & \displaystyle \sum _{(i,j)\in {\delta }^{+}(S)}{b}_{{ij}},\\ {\lambda }_{F}(r,s) & = & \mathop{\min }\limits_{\{S\subset V\ | \ r\in S,s\in V\setminus S\}}C(\delta (S)).\end{array}\end{eqnarray*}$$

By virtue of the max-flow-min-cut theorem [51], ${\lambda }_{F}(r,s)$ is equivalent to the maximum flow which can be transmitted from r to s respecting link capacity limits:

$$\begin{eqnarray}&&\begin{array}{l}{\lambda }_{F}(r,s)=\mathop{\max }\limits_{{\boldsymbol{F}}}\displaystyle \sum _{n=1}^{N}{F}_{r\to n}\\ {\rm{}}\,{\rm{such}}\,{\rm{that}}\,| {F}_{{mn}}| \leqslant {b}_{{mn}}\,\forall \,{\rm{edges}}\,(m,n)\\ {\rm{and}}\ \ \displaystyle \sum _{n=1}^{N}{F}_{{mn}}=0\,\forall \,m\ne r,s.\end{array}\end{eqnarray} \tag{ 29 }$$

Numerous efficient algorithms exist to calculate this maximum flow without performing the optimization explicitly [51]. The ratio ${b}_{{rs}}/{\lambda }_{F}(r,s)$ then gives the ratio of direct flow to total flow from r to s and thus provides an adequate topology-based estimate for the locality factor $\eta (r,s)$. Indeed, we can prove that it provides a rigorous lower bound.

Proposition 1. The algebraic locality factor $\eta (r,s)$ is bounded by

$$\begin{eqnarray}&&\displaystyle \frac{{b}_{{rs}}}{{\lambda }_{F}(r,s)}\leqslant \eta (r,s)\leqslant 1.\end{eqnarray} \tag{ 30 }$$

A proof is given in appendix A. Numerical simulations for several test grids reported in figure 4 reveal that the topological estimate not only provides a lower bound, but a high-quality estimate for the algebraic locality factor. The Pearson correlation coefficient ρ between $\eta (r,s)$ and ${b}_{{rs}}/{\lambda }_{F}(r,s)$ exceeds 0.92 for the three grids under consideration.

Zoom In Zoom Out Reset image size

We arrive at the conclusion that the dipole strength given by ${F}_{{rs}}{\left(1-\eta (r,s)\right)}^{-1}$ generally decreases with the redundancy measures ${\lambda }_{T}(r,s)$ and ${\lambda }_{F}(r,s)$.

An upper limit for the locality factor $\eta (r,s)$ can be obtained from an elementary topological distance measure. We consider the weighted geodesic distance of the two nodes r and s after the failure of the direct link $(r,s)$, which we denote by ${\mathrm{dist}}_{1}^{{\rm{w}}}(r,s)$. The superscript w stands for weighted distance, the subscript 1 for the distance measured in the graph after removal of the link $(r,s)$. We then have the following upper bound.

Proposition 2. The algebraic locality factor $\eta (r,s)$ is bounded from above by

$$\begin{eqnarray}&&\eta (r,s)\leqslant {\left[1+\displaystyle \frac{1}{{b}_{{rs}}\times {\mathrm{dist}}_{1}^{{\rm{w}}}(r,s)}\right]}^{-1}.\end{eqnarray} \tag{ 31 }$$

A proof is given in appendix B. Numerical simulations for several test grids reported in figure 5 reveal that the estimate in terms of the shortest path length not only provides an upper bound, but a high quality estimate for the algebraic locality factor. The Pearson correlation coefficient ρ exceeds 0.94 for the three grids under consideration.

Zoom In Zoom Out Reset image size

We further note that the factor ${{\boldsymbol{\nu }}}_{{rs}}^{\top }{{\boldsymbol{B}}}^{\dagger }{{\boldsymbol{\nu }}}_{{rs}}$ can also be interpreted as a distance measure—the resistance distance [52, 53]. We come back to the quantification of distances in flow networks later in section 6.3.

We now turn to the spatial aspects of flow rerouting in general network topologies. We first discuss some rigorous results, showing how the network topology determines the rerouting flows. Then, we return to the regular tilings and study the effect of increasing sparsity in these topologies on the dipole scaling. Finally, we suggest a new measure of distance for flow rerouting and examine its performance on realistic network topologies taken from power grids.

6.1. Rigorous results

To start off, we first present a lemma due to Shapiro [54], relating the flow changes after a link failure in an unweighted graph solely to the topology of the underlying network.

Lemma 1. Consider an unweighted network with a unit dipole source along the edge $(r,s)$, i.e. a unit inflow at node r and unit outflow at node s. Then the flow along any other edge $(m,n)$ is given by

$$\begin{eqnarray*}&&{F}_{m\to n}=\displaystyle \frac{{ \mathcal N }(r,m\to n,s)-{ \mathcal N }(r,n\to m,s)}{{ \mathcal N }},\end{eqnarray*}$$

where ${ \mathcal N }(r,m\to n,s)$ is the number of spanning trees that contain a path from r to s of the form $r,\ldots ,\ m,n,\ldots ,\ s$ and ${ \mathcal N }$ is the total number of spanning trees of the graph.

This lemma exactly gives the LODFs in terms of purely topological properties—the number of spanning trees containing certain paths. A generalization of this theorem to weighted graphs was recently presented in [55].

However, counting spanning trees is typically a difficult task such that these results are of limited use for practical applications. Nevertheless, they reveal the importance of certain paths through networks which we will analyze numerically in more detail below. Before we turn to this issue, we derive some weaker, but more easily applicable rigorous results.

We expect that the flow changes ${\rm{\Delta }}{F}_{{mn}}$ decay with distance as for the case of the square lattice analyzed in section 4. Can we establish some rigorous results on the decay with distance for arbitrary networks? Consider the outage of a single edge and assume that the network remains connected afterwards. We label the failing link as $(r,s)$ such that ${F}_{r\to s}\gt 0$ w.l.o.g. We first consider the change of the nodal potential or voltage phase angle ${\psi }_{n}$ and its decay with distance to the failing link $(r,s)$. More specifically, we define the maximum and minimum values of ${\psi }_{n}$ attained at a given distance:

$$\begin{eqnarray*}\begin{array}{rcl}{u}_{d} & = & \mathop{\max }\limits_{n,{\mathrm{dist}}_{0}^{{\rm{u}}}(n,r)=d}{\psi }_{n}\\ {{\ell }}_{d} & = & \mathop{\min }\limits_{n,{\mathrm{dist}}_{0}^{{\rm{u}}}(n,s)=d}{\psi }_{n}.\end{array}\end{eqnarray*}$$

Here, ${\mathrm{dist}}_{0}^{{\rm{u}}}(n,r)$ denotes the geodesic distance between two nodes n and r in the initial unweighted graph (indicated by the superscript u for unweighted and subscript 0 for the initial pre-contingency network). We then find the following result.

Proposition 3. Consider the failure of a single link $(r,s)$ with ${F}_{r\to s}\gt 0$ in a flow network. Then the maximum (minimum) value of the potential change ${\psi }_{n}$ decreases (increases) monotonically with the distance to nodes r and s, respectively:

$$\begin{eqnarray*}&&\begin{array}{l}{u}_{d}\leqslant {u}_{d-1},\quad 1\leqslant d\leqslant {d}_{\max }.\\ {{\ell }}_{d}\geqslant {{\ell }}_{d-1},\quad 1\leqslant d\leqslant {d}_{\max }.\end{array}\end{eqnarray*}$$

A proof is given in appendix C. We thus find that potential changes generally decrease with the distance in magnitude and so do the flow changes.

Furthermore, we can exploit the analogy to electrostatics to gain an insight into the scaling of flow changes with distance. As the flows are determined by a discrete Poisson equation, a discrete version of Gauss' theorem follows immediately. We note that we formulate this result in terms of the original network topology, see equation (18).

Lemma 2. Consider the failure of a single link $(r,s)$ in a flow network and denote by V the set of vertices in the network. For every decomposition of the network $V={V}_{1}+{V}_{2}$ with $r\in {V}_{1}$ and $s\in {V}_{2}$ we have

$$\begin{eqnarray}&&\displaystyle \sum _{m\in {V}_{1},n\in {V}_{2}}{\rm{\Delta }}{F}_{m\to n}={F}_{{rs}}\,{\left(1-\eta (r,s)\right)}^{-1}\,.\end{eqnarray} \tag{ 32 }$$

That is, for each decomposition the total flow between the two parts V₁ and V₂ equals the dipole strength.

This lemma supports the intuitive expectation that on average flow changes decay with distance in meshed networks: choose V₁ to include all nodes which are closer to r than to s and have a distance to r smaller than a given value

$$\begin{eqnarray*}&&{V}_{1}=\{n\in V| {{\rm{dist}}}_{0}^{{\rm{u}}}(r,n)\leqslant d;{{\rm{dist}}}_{0}^{{\rm{u}}}(r,n)\leqslant {{\rm{dist}}}_{0}^{{\rm{u}}}(s,n)\}.\end{eqnarray*}$$

With increasing value of d the number of nodes in V₁ increases and typically the number of edges between V₁ and V₂ increases, too. The total flow over these links remains constant according to lemma 2, such that the average flow will generally decrease. The exact scaling of the number of edges between V₁ and V₂ of course depends on the topology of the network.

One can furthermore show that a sufficient connectivity is needed for perturbations to spread. Generally, flow can be rerouted via an edge $(m,n)$ only if it can enter and leave the link via two independent paths. One can thus prove the following statement [55, 56].

Proposition 4. The line outage distribution factor ${{\rm{LODF}}}_{e,{\ell }}$ between two edges $e=(m,n)$ and ${\ell }=(r,s)$ vanishes if there are less than two independent paths between the vertex sets $\{r,s\}$ and $\{m,n\}$.

6.2. Impact of network topology

Now that we derived rigorous results on the scaling of LODFs, we want to study the influence of network connectivity on the scaling in more detail.

To do so, we first compare the scaling obtained for the square grid to the one in the other two regular tilings of two-dimensional space, namely the hexagonal grid and the triangular grid. In perfect realizations of these grids, each node has a degree of ${\deg }_{\mathrm{hex}}=3$ and ${\deg }_{\mathrm{tri}}=6$, respectively, whereas the degree for the square grid reads ${\deg }_{\mathrm{sg}}=4$. In figure 2(d), the LODFs are evaluated for these three topologies with increasing geodesic distance from the failing edge located again in the center of the networks between the nodes at $({x}_{r},{y}_{r})=(0,0)$ and $({x}_{s},{y}_{s})=(1,0)$. The quadratic scaling with the geodesic distance in x-direction $\parallel x{\parallel }^{-2}$ (black, dotted line) is preserved for all three topologies, i.e. the triangular grid (green triangles, bottom), the square grid (red squares, center) and the hexagonal grid (blue hexagons, top). The grids used here were of size 1000 × 1000 and 1001 × 500 nodes for the square grid and the triangular grid, respectively, and 150 × 150 hexagons for the hexagonal grid.

Thus, the quadratic scaling is robust throughout different regular networks. However, real networks are in general not regular. For this reason, we proceed by studying the effect of increasing sparsity in these regular tilings. Define the sparsity $\xi \in [0,1]\subset {\mathbb{R}}$ as the fraction of edges removed from the original graph. We make use of two different methods to achieve increasing sparsity. Our first method is a completely random removal of edges in the graph followed by measuring the LODFs along a specified path. If an edge along the path does no longer exist, we simply skip the edge. The results obtained from this method are shown in figures 2(a)–(c). There is no change visible in the scaling of LODFs, except for the direction perpendicular to the dipole in panel (c). In particular, only small values of sparsity ξ can be studied using this method, since a random removal of edges may easily result in disconnected graphs. For this reason, we make use of another method.

For the second method, we first construct an arbitrary spanning tree of the network after removal of the failing edge. Then, we subsequently remove random edges from the graph that are not part of the tree until a fraction ξ of its original edges is removed from the graph. This way, we make sure that the whole graph stays connected at all times. We continue by constructing the shortest path from the failing edge $((0,0),(1,0))$ to the node located at $({x}_{\max },0)$ and quantify the LODFs along this path. Note that using this method to make a graph sparser, we need to take into account the graph-specific maximal sparsity ${\xi }_{\max ,G}$, i.e. the fraction of edges whose removal would disconnect the graph. Assuming the initial tree to be minimal, this fraction may be calculated as ${\xi }_{\max ,\mathrm{hex}}=1/3$, ${\xi }_{\max ,\mathrm{sg}}=1/2$ and ${\xi }_{\max ,\mathrm{tri}}=2/3$ for the hexagonal grid, square grid and triangular grid, respectively.

Using this procedure, we can quantify the scaling of LODFs in grids with increasing sparsity. The direct assessment of a scaling exponent is difficult for sparser graphs due to the large spread in LODF values, see figure 6(a). This is why we construct a different measure to quantify this scaling. We consider the effective exponent $k(\xi )$, where ξ is the graph's sparsity, and assume a scaling of the form

$$\begin{eqnarray*}&&| \mathrm{LODF}(r,\xi )| \propto {r}^{-k(\xi )}\end{eqnarray*}$$

in some region of the geodesic distance $r=\parallel {\boldsymbol{r}}\parallel$ from the link failure. This effective exponent is calculated as follows

$$\begin{eqnarray*}&&k(\xi )=-{\mathrm{log}}_{5}\left(\displaystyle \frac{{\displaystyle \sum }_{r\in [5\times {10}^{1}-w,5\times {10}^{1}+w]}| \mathrm{LODF}(r,\xi )| }{{\displaystyle \sum }_{r\in [{10}^{1}-w,{10}^{1}+w]}| \mathrm{LODF}(r,\xi )| }\right),\end{eqnarray*}$$

where $w\in {\mathbb{N}}$ is a window specifying the range to average over in order to smooth the LODF values considered. We chose a window size of $w=2$ when calculating the effective exponent in practice which we found to result in a good compromise between smoothing and completely removing the trend. However, we did not observe a strong effect of the window size on the results. In addition to that, we chose to compare the LODFs at values centered around r = 10¹ and r = $5\times {10}^{1}$ when calculating the effective exponent since using this range allows us to capture only the intermediate range of the curve. For larger distances from the failing link, finite size effects prevent the assessment of the exponent whereas for smaller distances, the LODF does not yet decay when considering high values of sparsity due to a lack of alternative paths, as may also be observed in figure 6(a). For a perfect inverse square law $| \mathrm{LODF}| \propto \parallel {\boldsymbol{r}}{\parallel }^{-2}$ and a vanishing window $w=0$, this parameter yields $k=-{\mathrm{log}}_{5}({5}^{-2})=2$ as required. In figure 6(b), it can be observed that this effective exponent stays approximately constant at $k\approx 2$ over different values of sparsity and the three different topologies considered, where results for each value of sparsity were obtained using 100 random realizations of edge removals and with the same grid sizes as stated previously.

Zoom In Zoom Out Reset image size

To further quantify the effect of increasing sparsity in regular networks, we make use of another measure which we refer to as the LODF ratio ${R}_{w}(\xi )$. It is simply calculated as the logarithmic ratio between the LODFs with and without sparsity, again averaged over a fixed window of distances

$$\begin{eqnarray*}&&{R}_{w}(\xi )={\mathrm{log}}_{10}\left(\displaystyle \frac{{\displaystyle \sum }_{r\in [{10}^{1}-w,{10}^{1}+w]}| \mathrm{LODF}(r,\xi )| }{{\displaystyle \sum }_{r\in [{10}^{1}-w,{10}^{1}+w]}| \mathrm{LODF}(r,\xi =0)| }\right).\end{eqnarray*}$$

Note that we evaluate this parameter at a distance of 10¹ but we found the parameter to yield similar values for all distances considered. A parameter of ${R}_{w}(\xi )=1$ then represents a tenfold increase in the LODFs as compared to the network without any edges removed. In figure 6(c), this parameter is shown for the different topologies and sparsities. Here, a window size of $w=5$ was used. An increase with increasing sparsity is clearly visible. In particular, the LODFs increase on average more than tenfold close to the highest possible values of sparsity.

In total, we observe that the scaling exponent derived from the dipole analogy in section 5 holds for the regular networks even when removing a large fraction of their edges. On the other hand, the LODF values at a certain distance from the failing link show an increase with increasing sparsity, such that the actual effect of a link failure can be up to tenfold stronger than for the corresponding regular grid with no links removed. Thus, the overall effect of a link failure is more long-ranged in a sparser network, although no change in the effective exponent can be observed.

6.3. Scaling with distance

The impact of a link failure generally decays with distance. While the definition of distance is straightforward in regular lattices, different measures are meaningful in networks with complex topologies. The geodesic distance of two links follows from definition 1 for two vertices

$$\begin{eqnarray*}&&{\mathrm{edist}}_{\mathrm{ge}}^{{\rm{w}}}[(r,s),(m,n)]=\mathop{\min }\limits_{{v}_{1}\in \{r,s\},{v}_{2}\in \{m,n\}}{\mathrm{dist}}_{0}^{{\rm{w}}}({v}_{1},{v}_{2})+\displaystyle \frac{{w}_{{rs}}+{w}_{{mn}}}{2}.\end{eqnarray*}$$

Here, w_rs$\,=\,1/{b}_{rs}$ is the edge weight assigned to the edge $(r,s)$. When considering the unweighted analog, the edge distance is defined analogously setting all edge weights to one. The additional term $\tfrac{{w}_{{rs}}+{w}_{{mn}}}{2}$ ensures that neighboring edges have non-zero distance, e.g. unity distance ${\mathrm{edist}}_{\mathrm{ge}}^{{\rm{u}}}=1$ in the unweighted case. However, this distance is a bad indicator for flow rerouting in real-world irregular topologies. An example shown in figure 8 demonstrates that this simple distance is only weakly correlated with the magnitude of the LODFs for a real-world power grid test case.

Instead, we need a distance measure based on flow rerouting. If a link $(r,s)$ fails, the flow must be rerouted through other pathways, as described by the electrical lemma 1. However, it is not feasible to take into account all spanning trees which govern the flow rerouting. In order to still be able to estimate the impact on another link $(m,n)$, we will thus consider a path from r to s that crosses this link. The main difference to the ordinary graph theoretical distance is that we have to take into account a path back and forth. We are thus led to the following definition.

Definition 2. A ${\rm{rerouting}}\ {\rm{path}}$ from vertex $r$ to vertex $s$ via the edge $(m,n)$ is a path

$$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,{v}_{i}=m,{v}_{i+1}=n,{v}_{i+2},\ldots ,{v}_{k}=s)\end{eqnarray} \tag{ 33 }$$

or

$$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,{v}_{i}=n,{v}_{i+1}=m,{v}_{i+2},\ldots ,{v}_{k}=s),\end{eqnarray} \tag{ 34 }$$

where no vertex is visited twice. The ${\rm{rerouting}}\ {\rm{distance}}$ between two edges $(r,s)$ and $(m,n)$ denoted by ${\mathrm{edist}}_{\mathrm{re}}^{{\rm{u}}/{\rm{w}}}[(r,s),(m,n)]$ is the length of the shortest rerouting path from $r$ to $s$ via $(m,n)$ plus the length of edge $(r,s)$. Equivalently, it is the length of the shortest cycle crossing both edges $(r,s)$ and $(m,n)$. If no such path exists, the rerouting distance is defined to be $\infty$.

The definition of a rerouting path is illustrated in figure 7. Again, we consider a weighted and an unweighted version of this distance indicated by the superscript w and u, respectively. We note that the length of the edge $(r,s)$ is included in order to make the distance measure symmetric. In appendix D, we show explicitly that this definition satisfies the axioms of a metric and discuss how to compute the shortest rerouting path.

Zoom In Zoom Out Reset image size

An example of rerouting distances in comparison to the LODFs is shown in figure 8 for a small test grid. We observe a much better correlation in comparison to the ordinary geodesic distance defined above. The limitation of geodesic distances becomes especially clear for situations described by proposition 4. If exactly one independent path exists between two links, the rerouting distance is $\infty$, while the geodesic distance is finite. Hence, the latter fails to explain why the LODF between the two links vanishes.

Zoom In Zoom Out Reset image size

To further investigate the importance of distance, we simulate all possible link failures in four test grids of different size. For every failing link $(r,s)$ we evaluate the geodesic distance as well the rerouting distance to all other links in the grid. To quantify to which extend the distance predicts the magnitude of the LODFs, we then calculate the Kendall rank correlation coefficient τ [58]. This coefficient is used on ordinal data and assumes values in the interval $[-1,1]$. A value of (minus) one indicates perfect (anti)correlation, whereas a zero value implies no correlation between the data. Table 1 shows the results, averaging over all trigger links $(r,s)$ in the respective grid discarding bridges. The rank correlation is negative as LODFs generally decay with distance. The magnitude of the rank correlation is significantly higher for the rerouting distance. In particular for the test grid 'case1354pegase' we see that the ordinary geodesic distance has a very limited predictive power for the LODFs ($| \tau | \lt 0.25$), while the rerouting distance is strongly correlated to the magnitude of the LODFs ($| \tau | \gt 0.83$). Figure 9 illustrates this discrepancy in the distribution of τ values for the different distance measures for the test grids 'case118' and 'case1354pegase'. We are thus led to the conclusion that geodesic distances are of limited interest when considering the impact of link failures and should be replaced by other measures such as rerouting distances. Notably, we observe no major difference when comparing weighted and unweighted distances.

Zoom In Zoom Out Reset image size

Table 1.  Average of the Kendall τ rank correlation values for magnitude of LODF versus different distance measures. The four different IEEE test cases consistently show a higher degree of correlation between rerouting distances and LODF than between geodesic distances and LODF in both weighted and unweighted cases, while the unweighted rerouting distance slightly outperforms the weighted one. For examplary distributions of the τ values see figure 9.

test grid	Rank correlation τ for $| \mathrm{LODF}|$ versus distance
	Geodesic distance		Rerouting distance
	unw.	weighted	unw.	weighted
case30	−0.4027	−0.4015	−0.8528	−0.8440
case118	−0.6069	−0.5233	−0.8211	−0.7920
case1354pegase	−0.2269	−0.1341	−0.8664	−0.8438
case2383wp	−0.3604	−0.2318	−0.7213	−0.6066

Link failures represent major threats to the operation of complex supply networks across disciplines. In this article, we examined the impact of such failures in terms of the induced flow changes, which are commonly described by LODFs. We provide mathematically rigorous results and extensive numerical simulations with a focus on the gross network response (i.e. the dipole strength), the scaling of flow changes with distance and the role of network topology. These quantities are crucial to understand the global robustness of supply networks as each failure can trigger a cascade of secondary failures with potentially catastrophic consequences.

First, we demonstrated rigorously that the flow changes created by a single failure in a square lattice correspond to the field of an electromagnetic dipole. Hence the effects of a failure decay with the distance following an inverse-square law. The dipole analogy developed here allows for an analytical expression describing the spreading of link failures. Although this treatment is rigorously valid only in the continuum limit, we showed that the observed scaling extends to the other regular tilings of two-dimensional space even after removing a fraction of links. Thus, we conclude that the scaling may be expected to hold also for realistic topologies.

Increasing the sparsity of a network promotes more long-ranged effects up to the point where two links are only related by one independent pathway. Then, a rerouting between the two links becomes impossible and a failure of one link does not affect the other. However, this also implies a lack of redundancy such that a link failure can have catastrophic consequences locally. Our results thus suggest that sparsity promotes non-local responses to line failures. This is of potential relevance to the understanding of cascading failures, where previous outages increase sparsity, and deserves further study.

In real-world irregular networks, the gross response of a failure depends on the loading of the link as well as the local network structure. Rigorous upper and lower bounds were given for the dipole strength relating it to the redundancy of the failing link. Furthermore, the common notion of a geodesic graph distance is of limited use to predict flow rerouting. We thus introduced a rerouting distance which we showed to be much more meaningful to predict the impact of failures.

Whereas the classical analysis of link failures relies heavily on simulation results, our results provide heuristic methods and rigorous bounds which allow for an analytical insight into the relationship between the structure of a network and its robustness towards link failures. In particular for large networks where simulations are difficult, our results allow for an a priori analysis of link failures and might also be used to identify critical links, for instance in terms of the locality factor which quantifies the response of a network to a single failure. This type of analysis is aided by the general results on decay of maximal flow changes with geodesic and rerouting distances. We expect that these results fit the better, the more heterogeneous or disordered a network is. Previous studies [59] have shown that a strong heterogeneity of link parameters leads to a concentration of flows along the shortest path. In this limit, flow rerouting should be fully dominated by the shortest rerouting path.

We expect our results to be applicable far beyond power grids since the linearized treatment extends to other phenomena such as hydraulic or biological networks. The rerouting distance along with the bounds on the locality factor may greatly simplify the study of link failures in all kinds of supply networks and makes them more accessible. We expect our results on the scaling of LODFs for networks with increasing sparsity along with this distance measure to help identifying critical parts and paths and improving the overall robustness of supply networks.

We gratefully acknowledge support from the German Federal Ministry of Education and Research (BMBF grant no. 03SF0472) 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.).

Proof. By definition, $\eta (r,s)$ is given by the flow ${F}_{{rs}}/{\rm{\Delta }}P$ when the power ${\rm{\Delta }}P$ is injected at node r and withdrawn at node s, while there is no injection at any other node,

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{\displaystyle \sum }_{n}{F}_{r\to n} & = & {\displaystyle \sum }_{n}{F}_{n\to s}={\rm{\Delta }}P,\\ {\displaystyle \sum }_{n}{F}_{m\to n} & = & 0\,{\rm{\forall }}\,m\ne r,s,\end{array}\end{array}\end{eqnarray} \tag{ A.1 }$$

such that

$$\begin{eqnarray}&&\eta (r,s)=\displaystyle \frac{{F}_{r\to s}}{{\sum }_{n}{F}_{r\to n}}.\end{eqnarray} \tag{ A.2 }$$

For the sake of simplicity, we choose ${\rm{\Delta }}P$ such that ${\theta }_{r}-{\theta }_{s}=1$ w.l.o.g. Then, the inverse of $\eta (r,s)$ may be calculated using the basic relation ${F}_{r\to s}={b}_{{rs}}({\theta }_{r}-{\theta }_{s})$ as

$$\begin{eqnarray}&&\displaystyle \frac{{b}_{rs}}{\eta (r,s)}={\sum }_{n}{F}_{rn}.\end{eqnarray} \tag{ A.3 }$$

We can now use that the potential drop over all other links in the network is smaller than for the link $(r,s)$

$$\begin{eqnarray}&&| {\theta }_{m}-{\theta }_{n}| \leqslant {\theta }_{r}-{\theta }_{s},\end{eqnarray} \tag{ A.4 }$$

see proposition 3. If ${\theta }_{r}-{\theta }_{s}=1$ we thus know that

$$\begin{eqnarray}&&| {F}_{{mn}}| =| {b}_{{mn}}({\theta }_{m}-{\theta }_{n})| \leqslant {b}_{{mn}}.\end{eqnarray} \tag{ A.5 }$$

We thus obtain

$$\begin{eqnarray}&&\displaystyle \frac{{b}_{{rs}}}{\eta (r,s)}={\sum }_{n}{F}_{{rn}}\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\rm{such}}\,{\rm{that}}\,{\theta }_{r}-{\theta }_{s}=1,\\ {F}_{{mn}}={b}_{{mn}}({\theta }_{m}-{\theta }_{n}),\\ | {F}_{{mn}}| \leqslant {b}_{{mn}}\,\forall \,{\rm{edges}}\,(m,n),\\ \displaystyle \sum _{n=1}^{N}{F}_{{mn}}=0\,\forall \,m\ne r,s.\end{array}\end{eqnarray} \tag{ A.7 }$$

Comparing to the expression (29) for ${\lambda }_{F}(r,s)$ we see that two additional constraints have to be satisfied. Additional constraints can only decrease the flow-value with respect to the maximum in equation (29) such that we have

$$\begin{eqnarray}&&\displaystyle \frac{{b}_{{rs}}}{\eta (r,s)}\leqslant {\lambda }_{F}(r,s)\Rightarrow \eta (r,s)\geqslant \displaystyle \frac{{b}_{{rs}}}{{\lambda }_{F}(r,s)}.\end{eqnarray} \tag{ A.8 }$$

□

Proof. Consider first a reduced network consisting only of the link $(r,s)$ and the shortest alternative path between the two nodes, which we denote as $({j}_{1}=r,{j}_{2},{j}_{3},\cdots ,\ {j}_{n}=s)$. Fixing the nodal potentials such that ${\theta }_{r}-{\theta }_{s}=1$ as in A, the direct flow over the link $(r,s)$ is given by

$$\begin{eqnarray}&&F{{\rm{{\prime} }}}_{r\to s}={b}_{rs},\end{eqnarray} \tag{ B.1 }$$

whereas the indirect flow over the shortest alternative path is given by

$$\begin{eqnarray}\begin{array}{rcl}F{{\prime} }_{r\to {j}_{2}} & = & F{{\prime} }_{{j}_{2}\to {j}_{3}}=\cdots =F{{\prime} }_{{j}_{n-1}\to {j}_{n}}\\ & = & {\left[{b}_{{j}_{1},{j}_{2}}^{-1}+{b}_{{j}_{2},{j}_{3}}^{-1}+\cdots +{b}_{{j}_{n-1},{j}_{n}}^{-1}\right]}^{-1}\\ & = & \displaystyle \frac{1}{{\mathrm{dist}}_{1}^{{\rm{w}}}(r,s)}.\end{array}\end{eqnarray} \tag{ B.2 }$$

Now consider the initial network, to which all edges have been reintroduced, where we keep the same difference in nodal potentials ${\theta }_{r}-{\theta }_{s}=1$ which might require a different power injection ${\rm{\Delta }}P$. The direct flow thus remains the same while the total flow can only increase because new alternative paths may be present such that

$$\begin{eqnarray}&&\displaystyle \sum _{n}{F}_{r\to n}\geqslant F{{\prime} }_{r\to s}+F{{\prime} }_{r\to {j}_{2}}={b}_{{rs}}+\displaystyle \frac{1}{{\mathrm{dist}}_{1}^{{\rm{w}}}(r,s)}.\end{eqnarray} \tag{ B.3 }$$

Thus we obtain (see equation (A.2))

$$\begin{eqnarray}&&\eta (r,s)=\displaystyle \frac{{F}_{r\to s}}{\displaystyle \sum _{n}{F}_{r\to n}}\leqslant {\left[1+\displaystyle \frac{1}{{b}_{{rs}}\times {\mathrm{dist}}_{1}^{{\rm{w}}}(r,s)}\right]}^{-1}.\end{eqnarray} \tag{ B.4 }$$

□

In this appendix we first give the proof for proposition 3 and then show when the decay becomes strictly monotonous.

Proof. The proof is carried out by induction starting from $d={d}_{\max }$. We only give the proof for the maximum, the proof for the minimum proceeds in an analogous way. We assume that the network is large enough such that ${d}_{\max }\geqslant 2$, otherwise the statement is trivial anyway.

• (1)  
  Base case $d={d}_{\max }$: Consider the node n of the network for which $\mathrm{dist}(n,r)={d}_{\max }$ and ${\psi }_{n}$ assumes its maximum ${\psi }_{n}={u}_{{d}_{\max }}$. By assumption we have $\mathrm{dist}(n,r)\geqslant 2$ such that the node n cannot be adjacent to the perturbed edge such that ${q}_{n}=0$. The n-th component of equation (18) yields

  $$\begin{eqnarray*}\begin{array}{rcl}{B}_{{nn}}{\psi }_{n} & = & -\displaystyle \sum _{m\ne n}{B}_{{nm}}{\psi }_{m}\\ & = & -\displaystyle \sum _{m\ne n,\mathrm{dist}(m,r)={d}_{\max }}{B}_{{nm}}{\psi }_{m}-\displaystyle \sum _{m\ne n,\mathrm{dist}(m,r)={d}_{\max }-1}{B}_{{nm}}{\psi }_{m}.\end{array}\end{eqnarray*}$$

  We define the abbreviations

  $$\begin{eqnarray}&&{{ \mathcal B }}_{d}=-\sum _{m\ne n,\mathrm{dist}(m,r)=d}{B}_{{nm}}\end{eqnarray} \tag{ C.1 }$$

  and use some important properties of the matrix ${\boldsymbol{B}}$:

  $$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{B}_{nm} & \leqslant & 0\,{\rm{f}}{\rm{o}}{\rm{r}}\,n\ne m\,\,\Rightarrow \,\,{{B}}_{d}\geqslant 0,\\ {B}_{nn} & \geqslant & {{B}}_{{d}_{max}}+{{B}}_{{d}_{max-1}}.\end{array}\end{array}\end{eqnarray} \tag{ C.2 }$$

  We can furthermore bound the values of ${\psi }_{m}$ in equation (C) by ${u}_{{d}_{\max }}$ or ${u}_{{d}_{\max }-1}$, respectively, such that we obtain

  $$\begin{eqnarray}&&{u}_{{d}_{\max }}={\psi }_{n}\leqslant \displaystyle \frac{{{ \mathcal B }}_{{d}_{\max }}{u}_{{d}_{\max }}+{{ \mathcal B }}_{{d}_{\max }-1}{u}_{{d}_{\max }-1}}{{{ \mathcal B }}_{{d}_{\max }}+{{ \mathcal B }}_{{d}_{\max }-1}}\Rightarrow {u}_{{d}_{\max }}\leqslant {u}_{{d}_{\max }-1}.\end{eqnarray} \tag{ C.3 }$$
• (2)  
  Inductive step $d\to d-1$: We consider the node n of the network with ${\rm{dist}}(n,r)=d$ and ${\psi }_{n}={u}_{d}$. Starting from equation (18) and using the same estimates as above, we obtain

  $$\begin{eqnarray}&&{u}_{d}={\psi }_{n}=\displaystyle \frac{{q}_{n}-\displaystyle \sum _{m\ne n}{B}_{{nm}}{\psi }_{m}}{{B}_{{nn}}}\leqslant \displaystyle \frac{{{ \mathcal B }}_{d-1}{u}_{d-1}+{{ \mathcal B }}_{d}{u}_{d}+{{ \mathcal B }}_{d+1}{u}_{d+1}}{{{ \mathcal B }}_{d-1}+{{ \mathcal B }}_{d}+{{ \mathcal B }}_{d+1}}.\end{eqnarray} \tag{ C.4 }$$

  Note that the inhomogeneity ${q}_{n}\leqslant 0$ for all nodes except for $n=r$. With the induction hypothesis ${u}_{d+1}\leqslant {u}_{d}$ this yields

  $$\begin{eqnarray}&&{u}_{d}\leqslant \displaystyle \frac{{{B}}_{d-1}{u}_{d-1}+({{B}}_{d}+{{B}}_{d+1}){u}_{d}}{{{B}}_{d-1}+{{B}}_{d}+{{B}}_{d+1}}\Rightarrow {u}_{d}\leqslant {u}_{d-1},\end{eqnarray} \tag{ C.5 }$$

  which completes the proof.

□

The rerouting distance introduced in definition 2 is a proper distance measure in the sense that it satisfies the axioms of a metric as shown in the following lemma. It can be calculated by mapping it to the two-edge disjoint shortest path problem, which can be solved by Suurballe's algorithm [60]. The mapping is provided by the lemma 4.

Lemma 3. Consider an undirected graph with non-negative (all-equal) edge weights. Then the rerouting distance ${\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(r,s),(m,n)]$ of two edges $(r,s)$ and $(m,n)$ satisfies the following properties.

• (i)  
  Positive definiteness:

  $$\begin{eqnarray*}&&{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(r,s),(m,n)]\geqslant 0.\end{eqnarray*}$$
• (ii)  
  Symmetry:

  $$\begin{eqnarray*}&&{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(r,s),(m,n)]={\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(m,n),(r,s)].\end{eqnarray*}$$
• (iii)  
  Triangular inequality:

  $$\begin{eqnarray*}&&{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(a,b),(r,s)]\leqslant {\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(a,b),(m,n)]+{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(m,n),(r,s)],\end{eqnarray*}$$

both in the weighted and unweighted case.

Proof. 

• (1)  
  Positive definiteness: as long as all edge weights are non-negative, all paths lengths and hence also the rerouting distances are non-negative.
• (2)  
  Symmetry: suppose

  $$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,{v}_{i}=m,{v}_{i+1}=n,\ldots {v}_{k}=s)\end{eqnarray} \tag{ D.1 }$$

  is the shortest rerouting path from r to s via $(m,n)$. Then

  $$\begin{eqnarray}&&({v}_{i+1}=n,{v}_{i+2},\ldots ,{v}_{k}=s,{v}_{0}=r,{v}_{1},\ldots ,\ {v}_{i}=m)\end{eqnarray} \tag{ D.2 }$$

  is also a rerouting path from n to m via $(r,s)$. One can then show that this must be the shortest such rerouting path via contradiction. So suppose that another path from n to m via $(r,s)$,

  $$\begin{eqnarray}&&({u}_{j+1}=n,{u}_{j+2},\ldots {u}_{{\ell }}=s,{u}_{0}=r,{u}_{1},\ldots ,\ {u}_{j}=m),\end{eqnarray} \tag{ D.3 }$$

  is shorter. Then the path

  $$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,{v}_{i}=m,{v}_{i+1}=n,\ldots {v}_{{\ell }}=s)\end{eqnarray} \tag{ D.4 }$$

  is a rerouting path from r to s via $(m,n)$ and it is shorter than than the one defined in equation (D.1). This contradicts our initial assumption such that the path defined in equation (D.1) is the shortest rerouting path from n to m via $(r,s)$ and we obtain

  $$\begin{eqnarray}&&{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(r,s),(m,n)]={\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(m,n),(r,s)].\end{eqnarray} \tag{ D.5 }$$
• (3)  
  Triangle inequality: let the paths

  $$\begin{eqnarray*}\begin{array}{rcl}{p}_{1} & = & \left({v}_{0}=a,{v}_{1},\ldots ,{v}_{i}=m,{v}_{i+1}=n,\ldots {v}_{{\ell }}=b\right),\\ {p}_{2} & = & \left({u}_{0}=m,{u}_{1},\ldots ,{u}_{j}=r,{u}_{j+1}=s,\ldots {u}_{{\ell }}=n\right)\end{array}\end{eqnarray*}$$

  be the shortest rerouting paths from a to b via edge $(m,n)$ and from m to n via $(r,s)$, respectively. Here, we assume the paths to be oriented as ${v}_{i}=m$, ${v}_{i+1}=n$ and ${u}_{i}=r,{u}_{i+1}=s$, but the proof is the same if the order of these vertices in the path is reversed. In addition to that, we assume the two distances on the right-hand side of the inequality to be finite, otherwise the proof is trivial. We can extend the path p₂ to become a cycle by adding the edge $(n,m)$ to the end of the path

  $$\begin{eqnarray*}&&{c}_{2}=({u}_{0}=m,\ldots ,\ {u}_{i}=r,{u}_{i+1}=s,\ldots {u}_{{\ell }}=n,{u}_{{\ell }+1}=m).\end{eqnarray*}$$

  Now we can explicitly construct a rerouting path from a to b via $(r,s)$. Let ${u}_{j}\equiv {v}_{p}$ be the first vertex that appears in both p₁ and c₂ and let ${u}_{k}\equiv {v}_{q}$ the last such vertex. In this case, one of the following paths is a rerouting path from a to b via $(r,s)$

  $$\begin{eqnarray*}\begin{array}{rcl}{p}_{3} & = & \left({u}_{0}=a,{u}_{1},\ldots ,\ {u}_{j}={v}_{p},{v}_{p+1},\ldots ,{v}_{q-1},{u}_{k}={v}_{q},\ldots {u}_{{\ell }}=b\right)\\ \mathrm{or}\ {p}_{4} & = & \left({u}_{0}=a,{u}_{1},\ldots ,\ {u}_{j}={v}_{p},{v}_{p-1},\ldots ,{v}_{q+1},{u}_{k}={v}_{q},\ldots {u}_{{\ell }}=b\right).\end{array}\end{eqnarray*}$$

  Assume without loss of generality that p₃ is a rerouting path from a to b via $(r,s)$. In this case, we obtain

  $$\begin{eqnarray*}&&\begin{array}{l}{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(a,b),(r,s)]\leqslant \mathrm{length}((a,b))+\mathrm{length}({p}_{3})\\ \quad \leqslant \,\mathrm{length}((a,b))+\mathrm{length}({p}_{1})+\mathrm{length}({c}_{2})\\ \quad =\,\mathrm{length}((a,b))+\mathrm{length}({p}_{1})+\mathrm{length}((m,n))+\mathrm{length}({p}_{2})\\ \quad =\,{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(a,b),(m,n)]+{\mathrm{edist}}_{\mathrm{re}}^{{\rm{w}}/{\rm{u}}}[(m,n),(r,s)].\end{array}\end{eqnarray*}$$

  Note that again the length of a path is the sum of the edge weights of all edges in the path when considering a weighted graph.

□

Lemma 4. The shortest rerouting path form r to s via edge $(m,n)$ is given by the union of the edge $(m,n)$ and the two edge-independent paths $r\to m$ and $n\to s$ or $r\to n$ and $m\to s$ which minimize the total path length.

Proof. Assume that we have found a solution to the two-disjoint shortest path problem, i.e. we have found two edge-independent paths

$$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,\ {v}_{i}=m)\ \mathrm{and}\ ({u}_{0}=s,{u}_{1}\ldots {u}_{j}=s),\end{eqnarray} \tag{ D.6 }$$

which minimize the total path length. By assumption the two paths are edge-independent such that

$$\begin{eqnarray}&&({v}_{0}=r,{v}_{1},\ldots ,{v}_{i}=m,{u}_{0}=s,{u}_{1}\ldots {u}_{j}=s)\end{eqnarray} \tag{ D.7 }$$

is a valid rerouting path. Now it remains to show that this path is indeed the shortest possible. So assume the contrary, i.e. that there exists a path

$$\begin{eqnarray}&&({w}_{0}=r,{w}_{1},\ldots ,{w}_{i}=m,{w}_{i+1}=n,\ldots {w}_{k}=s)\end{eqnarray} \tag{ D.8 }$$

which is shorter than (D.7). But then the two paths

$$\begin{eqnarray}&&({w}_{0}=r,{w}_{1},\ldots ,{w}_{i}=m)\end{eqnarray} \tag{ D.9 }$$

$$\begin{eqnarray}&&({w}_{i+1}=n,{w}_{i+2}\ldots {w}_{k}=s)\end{eqnarray} \tag{ D.10 }$$

are edge independent and have a shorter total path length than the two paths D.6. Contradiction. □