Topologically protected loop flows in high voltage AC power grids

Power grids are networks of electrical lines whose purpose is to deliver electric power from producers to consumers. The ensuing power flows do not usually follow specified paths, instead they divide among all possible paths following Kirchhoff's laws. Circulating loop flows around closed, geographically constrained loops have been observed in the North American high voltage power grid, sometimes reaching as much as 1GW [1, 2], delivering no power but dissipating it ohmically. To try and prevent them, grid operators have issued new market regulations and recommended integrating phase angle regulating transformers into the grid [3]. The network conditions under which circulating power flows emerge, why they are so robust, how much power they dissipate and whether specific network topologies, if any, could prevent them in the first place are issues of paramount importance which have not been addressed to date. At a conceptual level, the definition of circulating loop flows is furthermore ambiguous, being arbitrarily based on an ill-defined separation of power flows into direct, parallel-path and circulating loop flows. Our goal in this manuscript is to understand better the nature of these circulating power flows.

A deep and unexpected analogy between high voltage electric power transmission and macroscopic quantum states such as superfluids and superconductors has been overlooked so far. The operational state of an AC power grid is determined by the complex voltage at each bus, ${V}_{l}=| {V}_{l}| \exp [i{\theta }_{l}]$ which must be single-valued. Summing over voltage angle differences around any loop in the network must therefore give an integer multiple of $2\pi$. This defines topological winding numbers q_α ,

$$\begin{eqnarray}&&{q}_{\alpha }={(2\pi )}^{-1}\displaystyle \sum _{l=1}^{{n}_{\alpha }}| {\theta }_{l+1}-{\theta }_{l}| \,\in \,{\mathbb{Z}},\end{eqnarray} \tag{ 1 }$$

where the sum runs over all n_α nodes around any (the $\alpha \mathrm{th}$) loop in the network, $| {\theta }_{l+1}-{\theta }_{l}|$ gives the angle difference along the $l\mathrm{th}$ line in this loop, counted modulo $2\pi$, and node indices are taken modulo n_α , i.e. ${n}_{\alpha }+1\to 1$. The topological meaning of q_α is obvious, as it counts the number of times the complex voltage winds around the origin in the complex plane as one goes around the $\alpha \mathrm{th}$ loop. The condition ${q}_{\alpha }\in {\mathbb{Z}}$ is the same as the condition that leads to quantization of circulation around superfluid vortices [4, 5] or to flux quantization through a superconducting ring [6].

The analogy with superconductivity is complete under the lossless line approximation, where AC transmission lines are assumed purely susceptive (see Supplemental Material) [7]. Then, the active power flowing between two nodes l and m is given by ${P}_{{lm}}={B}_{{lm}}| {V}_{l}| | {V}_{m}| \sin ({\theta }_{l}-{\theta }_{m})$, with the elements B_lm of the susceptance matrix. This is the DC Josephson current that would flow between two superconductors with phases ${\theta }_{l}$ and ${\theta }_{m}$, coupled by a short tunnel junction of transparency ${T}_{{lm}}={\hslash }{B}_{{lm}}| {V}_{l}| | {V}_{m}| /8e$ [8]. It has been shown within the lossless line approximation that, in complex networks, various solutions to the power flow problem exist, which differ only by circulating loop currents [9, 10]. Neglecting ohmic dissipation, it is therefore expected, and has been reported in simple networks [10–12], that AC power grids may carry circulating loop flows. Topological winding numbers, equation (1), lead to the discretization of circulating loop flows, much in the same way as superfluid circulation is quantized around a vortex [4, 5]. Therefore we refer to circulating loop flows as vortex flows from now on. The existence of integer winding numbers has the important consequence that the vortex flows are topologically protected, much in the same way as persistent currents in superconducting loops [13]. The integer q_α in equation (1) measures the number of times the complex voltage ${V}_{l}=| {V}_{l}| \exp [i{\theta }_{l}]$ rotates in the complex plane as one goes around a loop in the network. Changing a vortex flow requires to change the number q_α of such rotations, i.e. to untwist V_l , which cannot be done smoothly without driving $| {V}_{l}| \to 0$ somewhere. It is thus hard to get rid of a vortex flow without topological changes in operating electrical networks.

AC power grids however differ from superfluids and superconducting systems in at least two significant ways in that (i) whereas vortices are generated by external magnetic fields (in a superconductor) or sample rotation (in a superfluid), how to create vortex flows in AC power grids is not quite understood, and (ii) superfluids and superconductors are nondissipative quantum fluids, whereas AC power lines dissipate ohmically part of the power they transmit. The lossless line approximation is in fact only partially justified in very high voltage AC power grids, where lines have a conductance that is at least ten times smaller than their susceptance [7]. Still, ohmic losses typically reach 5%–10 % of the total transported power. It is therefore important to find out whether the above analogy between high voltage AC power grids and macroscopic quantum states is at all physically relevant. The two main purposes of this manuscript are therefore (i) to investigate how vortex flows can be created in electric power grids and (ii) to investigate how resilient they are to the presence of ohmic dissipation. Our investigations of creation mechanisms amplify on the work of Janssens and Kamagate [14] who succinctly discussed one of the three mechanisms we identify below. Investigating basins of stability for different solutions via the Lyapunov function allows us furthermore to shed analytical light on the line reclosing mechanism they proposed [14] and give precise bounds on when vortex flows are created in this way. We furthermore find that vortex flows are resilient to reasonable amounts of ohmic dissipation typical of high voltage power grids and that vortex-carrying operating states ohmically dissipate significantly more electric power than vortex-free states.

Loop flows in electric power grids have been investigated before in a number of theoretical and numerical works. We list some of the most important published works we know about. Korsak investigated a simple network where different, linearly stable solutions exist that differ by some circulating loop current [15]. Tavora and Smith related the existence of different stable fixed points of the power flow problem to the presence of integer winding numbers [16], reflecting the $2\pi$-periodicity of the complex voltage around any loop in the network and its single-valuedness. The characterization of circulating loop flows with topological winding numbers has been pushed further by Janssens and Kamagate [14], who also investigated how to generate such loop flows and found one of the three creation mechanisms we discuss below. More recently, [9, 10] showed that within the lossless line approximation, different power flow solutions must be related to one another by circulating loop flows.

While vortex flows and the analogy we just pointed out between macroscopic quantum states and high voltage AC power grids are intellectually interesting in their own right, we stress that they are physically and technically relevant. Circulating loop flows in the GW range have been observed in power grids [1, 2], which delivered no power but consumed it ohmically. Furthemore it can be expected that with the changes in operational conditions of the power grid brought about by the energy transition—substituting delocalized productions with smaller primary power reserve for large power plants—such circulating loop flows may occur more frequently. We stress, however, that, while power flow solutions are similar to vortex-carrying quantum mechanical states, there is no quantumness in the power flow problem and that high voltage AC power grids are not superconducting.

The paper is organized as follows. In section 2 we argue that different solutions to the power flow problem in meshed networks are related to one another via vortex flows. Sections 3, 4 and 5 next discuss sequentially the three mechanisms we identified for generating vortex flows. Section 5 in particular investigates vortex flow creation via reclosing of a line from a basin of attraction point of view, which allows to quantitatively understand how vortex flows are born. Section 6 makes the case that vortex flows are robust against the relatively modest amount of ohmic dissipation present in high voltage power networks. Having gained much understanding of vortex flows in simple models in these early sections, we next discuss vortex creation with and without dissipation on a complex network with the topology of the UK grid. This is done in section 7. Conclusions and future perspectives are briefly discussed in section 8.

We start with the lossless line approximation and neglect voltage variations, $| {V}_{l}| ={V}_{0}$, $\forall l$. Active power flows are then governed by the following set of equations (see Supplemental Material)

$$\begin{eqnarray}&&{P}_{l}=\displaystyle \sum _{m}{\tilde{B}}_{{lm}}\sin ({\theta }_{l}-{\theta }_{m}).\end{eqnarray} \tag{ 2 }$$

Here, P_l is the active power injected (${P}_{l}\gt 0$) or consumed (${P}_{l}\lt 0$) at node l, ${\theta }_{l}$ is the complex voltage angle and ${\tilde{B}}_{{lm}}={B}_{{lm}}{V}_{0}^{2}$. At this level, one has an exact balance between production and consumption, ${\sum }_{l}{P}_{l}=0$. The theorem of [9, 10] states that different solutions to equation (2) on a meshed network differ only by circulating flows around loops in the network. We go beyond the lossless line approximation to see how much validity this theorem keeps in the presence of ohmic losses. Being interested in high voltage grids we neglect voltage fluctuations all through this manuscript, as they correspond to few percents of the rated voltage V₀. Ohmic losses are introduced by rewriting equation (2) as

$$\begin{eqnarray}&&{P}_{l}=\displaystyle \sum _{m}({\tilde{B}}_{{lm}}\sin ({\theta }_{l}-{\theta }_{m})+{\tilde{G}}_{{lm}}[1-\cos ({\theta }_{l}-{\theta }_{m})]),\end{eqnarray} \tag{ 3 }$$

where ${\tilde{G}}_{{lm}}={G}_{{lm}}{V}_{0}^{2}$ with elements G_lm of the conductance matrix. Because of ohmic dissipation, total production now exceeds total consumption,

$$\begin{eqnarray}&&{\rm{\Delta }}P=\displaystyle \sum _{l}{P}_{l}=\displaystyle \sum _{l,m}{\tilde{G}}_{{lm}}[1-\cos ({\theta }_{l}-{\theta }_{m})]\gt 0.\end{eqnarray} \tag{ 4 }$$

Equation (4) implies in particular that different solutions to equation (3) dissipate different amounts ${\rm{\Delta }}P$ of active power and therefore require different power injections $\{{P}_{l}\}\to \{{P}_{l}+\delta {P}_{l}\}$ to compensate for ohmic losses. The set $\{\delta {P}_{l}\}$ is not uniquely defined, nevertheless, choices exist for which the theorem of [9, 10] remains valid (see Supplemental Material). This is in agreement with Baillieul and Byrnes [17] who stated that 'models for lossless power networks provide valuable insight and understanding for systems with small transfer conductances' and further suggests that circulating loop flows are robust against a moderate amount of ohmic dissipation. Below we numerically confirm this conjecture. Earlier works however suggested that the number of solutions decreases at fixed susceptance when the conductance increases [17, 18], so that it is expected that circulating loop flows are eventually suppressed when the conductance exceeds some network-dependent threshold.

Focusing next on the operational conditions under which circulating power flows occur, we find three different mechanisms for creating them, (i) by loss of stability of the solution carrying no circulating loop flow, (ii) by tripping of a line traversing a large loop in the network and (iii) by reclosing a loop that tripped or was open earlier [14, 19]. We discuss these three mechanisms sequentially, first without ohmic dissipation.

Consider first a single-ring, lossless network as illustrated in figure 1(d). The power flow around the ring is governed by equation (2) with ${\tilde{B}}_{{lm}}={B}_{0}$. From [10], the system carries at most nine solutions differing by loop flows, when ${B}_{0}\to \infty$. Figure 1(b) shows seven of these solutions, each characterized by a winding number $q=-3,-2,...3$. The stability of each solution is determined by the swing equations [7], which govern the dynamics of the voltage angles ${\theta }_{l}$ under changes in operating conditions. In this work we neglect inertia terms in swing equations, since their presence affects neither the nature of the stationary states, nor at which parameter values they become unstable [19, 20]. In a frame rotating at the grid frequency of 50 or 60 Hz, the swing equations read (see Supplemental Material) [7, 21–23]

$$\begin{eqnarray}&&\dot{{\theta }_{l}}={P}_{l}-\displaystyle \sum _{m=1}^{n}{B}_{0}\,\sin ({\theta }_{l}-{\theta }_{m}),\qquad l=1,\ldots n.\end{eqnarray} \tag{ 5 }$$

Stationary solutions to equation (5) are obviously solutions to equation (2) and their linear stability depends on the stability matrix M obtained after linearizing equation (5) (see Supplemental Material) [24]. A solution of equation (2) is stable if M is negative semidefinite. We therefore study the stability of each solution via the largest nonvanishing eigenvalue ${\lambda }_{2}$ of M. Figure 1(a) shows, together with figure 1(b) how solutions disappear as they lose their stability, ${\lambda }_{2}\to 0$. The solution with q = 0 has the smallest ${\lambda }_{2}$ at small P. Remarkably enough, the q = 0 solution loses its stability at $P/{B}_{0}\approx 1.6$, before the $q=-1$ solution, which remains stable until $P/{B}_{0}\approx 1.85$. Starting from the q = 0 solution and increasing P beyond 1.6B₀, we observe a loss of stability followed by a short transient after which the operating state has been transferred to the $q=-1$ state. This transient is illustrated in figures 1(c) and (d), which show that mostly one voltage angle, corresponding to the consumer node rotates while all other angles move very little (a movie of this transient can be found in the Supplemental Material). The rotation of this angle changes q which oscillates between q = 0 and $q=-1$, eventually stabilizing at $q=-1$.

Zoom In Zoom Out Reset image size

Reducing next P starting from the $q=-1$ solution at $P/{B}_{0}\gt 1.6$, one remains on the $q=-1$ solution. This hysteretic behavior is indicated by arrows in figures 1(a) and (b) and illustrates the topological protection brought about by the integer winding number q. We have found that this behavior is generic for single-ring networks (see Supplemental Material). The process by which the winding number changes is similar to quantum phase slips in small rings of Josephson junctions [25]. To emphasize this similarity, while stressing the different physical ingredients at work, we call dynamical phase slip this first mechanism of creation of circulating loop flows.

We next investigate the second mechanism for vortex flow generation, by tripping of a line. High voltage AC power grids have a meshed structure, where multiple paths connect production and consumption centers. This ensures that a single line failure does not preclude the supply of electric power. Consider the model shown in the inset of figure 2(a), where a producer is connected to a consumer via three different paths. Assume then that a line on the middle path trips. The power initially transmitted via that path is redistributed and if P is relatively large, the angle differences ${{\rm{\Delta }}}_{{\rm{L}},{\rm{R}}}$ on each remaining path increase significantly. When one of these two paths, say the left one, goes through many more lines than the other one, ${N}_{{\rm{L}}}\gg {N}_{{\rm{R}}}$, it is then possible that ${N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}-{N}_{{\rm{R}}}{{\rm{\Delta }}}_{{\rm{R}}}=2\pi q$ with $q\gt 0$, even if the system carried no vorticity initially. This simple example shows how one line tripping in an asymmetrical double-loop system can generate a vortex flow.

Figure 2 illustrates how a $q\ne 0$ state emerges from a q = 0 state after a line tripping. The initial state is a stationary state of the double-loop system with zero winding number on both loops. The red line in the inset of figure 2(a) is then cut, which induces a transient driving the system to a stationary state of the resulting single-loop system. Figure 2(a) shows the obtained winding number. One sees that for small P, the final state has q = 0, while for larger P, a q = 1 state is reached. We have found that this behavior is generic of sufficiently asymmetric double-loop systems. We discuss cases with ohmic dissipation below.

Zoom In Zoom Out Reset image size

One may wonder what is the fate of the q = 1 state created by line tripping when the line is reclosed. Line reclosing is a topological change that has the potential to induce integer changes in the winding number q so that a vortex-free state with q = 0 can be expected after line reclosing. We show in figure 3 that in the present case, line reclosing does not change the winding number and that the vortex flow persists, and that the winding number remains the same, q = 1. Inspecting the angles in the final state, we find that the vortex in the final state is supported by the larger, left loop.

Zoom In Zoom Out Reset image size

We have just showed that line tripping can lead to a vortex flow that is robust against the reclosing of the tripped line. In this paragraph, we finally consider vortex flow creation via reclosing of a line. We consider again the single-loop model sketched in the inset of figure 1(d). We consider here the case with only one consumer and one producer with some produced (consumed) power P ($-P$) but checked that our conclusions remain the same and that vortex formation proceeds similarly in more complicated single-loop models (see Supplemental Material). We start from the closed-loop system with an operating state with q = 0. The power P is transferred from producer to consumer both clockwise (via the right path) and counterclockwise (left). One line along the right path then trips, which forces P to be transmitted exclusively along the left path. The voltage angle difference ${{\rm{\Delta }}}_{{\rm{L}}}$ between any two nodes along the left path increases, while angle differences along the right path vanish, ${{\rm{\Delta }}}_{{\rm{R}}}=0$. This gives an angle difference ${{\rm{\Delta }}}_{0}={N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}$ between the two nodes on each side of the tripped line. Upon reclosing that line, a current flows through it whose initial direction depends on ${{\rm{\Delta }}}_{0}$, but which eventually relaxes following a dynamical process determined by equation (5). The creation of a vortex flow can however be understood without investigating the voltage angle dynamics, by instead adopting an approach based on the Lyapunov function [7, 22]. The latter determines the basin of attraction, in voltage-angle space, for the various solutions to the power flow problem [26, 27]. The Lyapunov function corresponding to equation (5) reads [28]

$$\begin{eqnarray}&&{ \mathcal V }(\{{\theta }_{i}\})=-\displaystyle \sum _{l}{P}_{l}{\theta }_{l}-\displaystyle \sum _{\langle l,m\rangle }{B}_{0}\,\cos ({\theta }_{l}-{\theta }_{m}),\end{eqnarray} \tag{ 6 }$$

where the second sum runs over connected nodes only. Minima of ${ \mathcal V }$ have ${{\rm{\nabla }}}_{\theta }{ \mathcal V }=0$, and thus correspond to stationary power flow solutions. The Lyapunov function can be rewritten as a function of angle differences ${{\rm{\Delta }}}_{l}:= {\theta }_{l}-{\theta }_{l+1}\in [-\pi ,\pi ]$ as

$$\begin{eqnarray}&&{ \mathcal V }(\{{{\rm{\Delta }}}_{i}\})=-\displaystyle \sum _{l}{P}_{l}^{* }{{\rm{\Delta }}}_{l}-\displaystyle \sum _{l}{B}_{0}\cos ({{\rm{\Delta }}}_{l}),\end{eqnarray} \tag{ 7 }$$

where ${P}_{l}^{* }:= {\sum }_{j=1}^{l}{P}_{j}$, for l = 1,...,n. The single-loop model we consider here can be split in two paths (left and right) from producer to consumer. We have

$$\begin{eqnarray}{P}_{l}^{* }=\left\{\begin{array}{ll}P & \mathrm{if}\ \mathrm{the}\ \mathrm{line}\ \mathrm{from}\ l\ \mathrm{to}\ l+1\ \mathrm{is}\ \mathrm{on}\ \mathrm{the}\ \mathrm{left}\ \mathrm{path},\\ 0 & \mathrm{if}\ \mathrm{the}\ \mathrm{line}\ \mathrm{from}\ l\ \mathrm{to}\ l+1\ \mathrm{is}\ \mathrm{on}\ \mathrm{the}\ \mathrm{right}\ \mathrm{path}.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

Furthermore, any solution of the nondissipative power flow equations has the same angle differences, ${{\rm{\Delta }}}_{{\rm{L}}}={\theta }_{l}-{\theta }_{l+1}$, along each line on the left path and ${{\rm{\Delta }}}_{{\rm{R}}}={\theta }_{l+1}-{\theta }_{l}$ on the right path.

Going around the loop, the voltage phases must be well defined. Therefore, just before line reclosing the phase difference ${{\rm{\Delta }}}_{0}$ between the two ends of the tripped line can be written as a function of ${{\rm{\Delta }}}_{{\rm{L}}}$ and ${{\rm{\Delta }}}_{{\rm{R}}}$, ${{\rm{\Delta }}}_{0}={N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}-({N}_{{\rm{R}}}-1){{\rm{\Delta }}}_{{\rm{R}}}-2\pi q$, where ${N}_{{\rm{L}}}\gt {N}_{{\rm{R}}}\geqslant 2$ are the number of edges on the left and right paths. Then we can project the Lyapunov function on the $({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})$-plane,

$$\begin{eqnarray}&&{ \mathcal V }({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})=-{N}_{{\rm{L}}}P{{\rm{\Delta }}}_{{\rm{L}}}-{N}_{{\rm{L}}}{B}_{0}\cos {{\rm{\Delta }}}_{{\rm{L}}}-({N}_{{\rm{R}}}-1){B}_{0}\cos {{\rm{\Delta }}}_{{\rm{R}}}-{B}_{0}\cos ({N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}-({N}_{{\rm{R}}}-1){{\rm{\Delta }}}_{{\rm{R}}}).\end{eqnarray} \tag{ 9 }$$

Figures 4(c) and (d) show contour plots of ${ \mathcal V }({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})$. Local minima are indicated, together with the corresponding integer winding numbers. To each minimum corresponds a basin of attraction containing the set of initial states that converge towards that minimum under equation (5). All points around a minimum belong to that basin, until one reaches a saddle point or a ridge, beyond which points belong to another basin of attraction. Cutting the right path projects ${{\rm{\Delta }}}_{{\rm{R}}}\to 0$. Right before line reclosing, the system is at ${{\rm{\Delta }}}_{{\rm{L}}}=\arcsin (P/{B}_{0})$ on the dashed lines in figures 4(c) and (d) which correspond to $P/{B}_{0}=0.259$ and $P/{B}_{0}=0.707$ respectively. The solution towards which the system converges after line reclosing depends on the basin of attraction to which the initial state belongs. For P such that ${N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}={{\rm{\Delta }}}_{0}=(2p+1)\pi$ with $p\in {\mathbb{Z}}$, the point $({{\rm{\Delta }}}_{L},{{\rm{\Delta }}}_{{\rm{R}}})=(\arcsin (P/{B}_{0}),0)$ lies right on a saddle point at the boundary between two basins of attraction, as we now proceed to show.

Zoom In Zoom Out Reset image size

The gradient of the Lyapunov function ${ \mathcal V }$ in the $({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})$-plane is given by

$$\begin{eqnarray*}&&{\rm{\nabla }}{ \mathcal V }=\left(\begin{array}{c}-{N}_{{\rm{L}}}P+{N}_{{\rm{L}}}{B}_{0}\sin {{\rm{\Delta }}}_{{\rm{L}}}+{N}_{{\rm{L}}}{B}_{0}\sin ({N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}-({N}_{{\rm{R}}}-1){{\rm{\Delta }}}_{{\rm{R}}})\\ ({N}_{{\rm{R}}}-1){B}_{0}\sin {{\rm{\Delta }}}_{{\rm{R}}}-({N}_{{\rm{R}}}-1){B}_{0}\sin ({N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}-({N}_{{\rm{R}}}-1){{\rm{\Delta }}}_{{\rm{R}}})\end{array}\right).\end{eqnarray*}$$

It is easy to check that ${\rm{\nabla }}{ \mathcal V }=0$ at $({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})=(\arcsin (P/{B}_{0}),0)$, which is thus a critical point. The nature of this critical point is determined by the two eigenvalues of the Hessian ${{ \mathcal H }}_{{ \mathcal V }}$ of ${ \mathcal V }$. At our critical point, we obtain

$$\begin{eqnarray*}{{ \mathcal H }}_{{ \mathcal V }}(\arcsin (P/{B}_{0}),0)={B}_{0}\left(\begin{array}{cc}{N}_{{\rm{L}}}\cos {{\rm{\Delta }}}_{{\rm{L}}}-{N}_{{\rm{L}}}^{2} & {N}_{{\rm{L}}}({N}_{{\rm{R}}}-1)\\ {N}_{{\rm{L}}}({N}_{{\rm{R}}}-1) & ({N}_{{\rm{R}}}-1)-{({N}_{{\rm{R}}}-1)}^{2}\end{array}\right)=: \left(\begin{array}{cc}a & b\\ b & c\end{array}\right).\end{eqnarray*}$$

The two eigenvalues of ${{ \mathcal H }}_{{ \mathcal V }}$ are then the two roots of

$$\begin{eqnarray*}&&\chi (\lambda )={\lambda }^{2}-(a+c)\lambda +{ac}-{b}^{2}\Rightarrow {\lambda }^{\pm }=(a+c\pm \sqrt{({(a+c)}^{2}-4{ac}-{b}^{2})})/2.\end{eqnarray*}$$

Then ${\lambda }^{+}$ is always positive and ${\lambda }^{-}$ is negative if and only if ${ac}-{b}^{2}\lt 0$. Replacing a, b and c we have

$$\begin{eqnarray}&&{ac}-{b}^{2}=-{B}_{0}^{2}{N}_{{\rm{L}}}({N}_{{\rm{R}}}-1)({N}_{{\rm{L}}}+({N}_{{\rm{R}}}-2)\cos {{\rm{\Delta }}}_{{\rm{L}}}),\end{eqnarray} \tag{ 10 }$$

which is necessarily negative since first, at the moment of the reclosing, ${{\rm{\Delta }}}_{{\rm{L}}}=\arcsin (P/{B}_{0})$ implying that $\cos {{\rm{\Delta }}}_{{\rm{L}}}\gt 0$, and second ${N}_{{\rm{R}}}\geqslant 2$. We conclude that $({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})=(\arcsin (P/{B}_{0}),0)$ for ${N}_{{\rm{L}}}\arcsin (P/{B}_{0})=(2p+1)\pi$ is a saddle point of the projected Lyapunov function. It can actually be shown that it is a saddle point of the full Lyapunov function. One concludes that vortex generation by this mechanism occurs for ${{\rm{\Delta }}}_{0}\gt \pi$, and that the final winding number increases by one each time ${{\rm{\Delta }}}_{0}$ crosses an odd integer multiple of π. The same line of argument with ${N}_{{\rm{L}}}{{\rm{\Delta }}}_{{\rm{L}}}\leftrightarrow {N}_{{\rm{R}}}{{\rm{\Delta }}}_{{\rm{R}}}$ applies when the line to be cut is on the left path.

The above argument is based on the projected Lyapunov function. It neglects the fact that, after line reclosing, the transient dynamics leaves the $({{\rm{\Delta }}}_{{\rm{L}}},{{\rm{\Delta }}}_{{\rm{R}}})$-plane until a new stationary state is reached. We therefore check its validity numerically. Figure 4 shows what final winding number is obtained upon line reclosing depending on the location D = 1,...,n of the open line (counted counterclockwise, starting from the main producer) and the rescaled power $P/{B}_{0}$ (figure 4(a)) and ${{\rm{\Delta }}}_{0}$ (figure 4(b)). Figure 4(b) confirms that the final winding number changes by one each time ${{\rm{\Delta }}}_{0}$ crosses an odd integer multiple of π, except when the injected power gets close to its maximal allowed value, $P\to {B}_{0}$. We attribute this change of behavior to a more complicated transient in this case. Figures 4(a) and (b) further show that ${{\rm{\Delta }}}_{0}\approx 2\pi$ around $P/{B}_{0}\approx 0.5$. Taken modulo $2\pi$, this means that the angle difference at the ends of the tripped line is small so that line reclosing is technically feasible.

Figures 4(c) and (d) shed a new light on the work of Araposthatis et al [29], who discussed the existence of different power flow solutions in separated stable domains in voltage angle space in simple models. Our method allows to visualize different such domains and in particular to infer the precise location of saddle points separating them. Recent works have advocated a new line of research in dynamical systems including coupled oscillators models [26, 27] and AC power networks [30], investigating the size of basins of attraction. These works were restricted to numerical statistical studies. Our projective approach allows to visualize basins of attractions and the separatrices in between. Quite remarkably, it allows us to understand quantitatively how vortices emerge and when winding numbers change. In both line tripping and line reclosing mechanisms, vortices are created by a topological change in the network, which twists the voltage angles around a loop. We therefore collectively refer to these two mechanisms as topological phase twist.

We next investigate the persistence of circulating loop flows in the presence of ohmic dissipation. Voltage angle differences between connected nodes in operational states of AC power grids seldomly reach more than few tens of degrees, beyond which the power line's thermal limit is exceeded. Therefore, the dynamical phase slip mechanism (i) is of little relevance for power grid operation, because it occurs when one angle difference $| {\theta }_{i}-{\theta }_{j}| \gtrsim \pi /2$ [10]. We therefore focus on the topological phase twist mechanisms (ii) and (iii).

The green and red lines in figure 2 illustrate how ohmic dissipation affects vortex flow creation by tripping of a line. Because ohmic dissipation requires an excess of production to compensate for losses, power production is $P+{\rm{\Delta }}P$, larger than the power demand P. One sees that the presence of a finite conductance, ${\tilde{G}}_{{lm}}\equiv {G}_{0}\ne 0$ in equation (3), reduces the range in $P/{B}_{0}$ at which transitions between q = 0 and q = 1 occur upon line tripping, but that the overall behavior remains the same as long as ${G}_{0}/{B}_{0}$ is not too large. Figure 2(c) shows furthermore how much more power is consumed in the presence of vortex flows, with a huge stepwise increase in ohmic losses by almost 50% of the losses ${\rm{\Delta }}P$ in the presence of a still moderate conductance ${G}_{0}/{B}_{0}=0.03$. Figure 2 finally shows topological protection by the winding number, where once the vortex has been created, returning the operating conditions to smaller P (as indicated by arrows) does not bring the system back to the vortex-free state. Instead, the operational state remains at q = 1, with losses well above those for q = 0.

Figure 5 makes it clear that vortex generation under the line tripping and reclosing mechanism (iii) proceeds in the same way in the presence of ohmic dissipation. We start from a q = 0 stationary state, cut a line and let the system relax to a stationary state, after which we close the line again. What final state is obtained depends on $P/{B}_{0}$. When P is small, the system relaxes back to the initial state with q = 0 (yellow area in figure 5(a)), however at larger P, the transient dynamically moves the system towards the q = 1 stationary state (green area in figure 5(a)). Figure 5(b) shows that the transition to q = 1 occurs precisely when the voltage angle difference between the nodes surrounding the faulted line reaches π, even in the presence of ohmic dissipation, in complete agreement with the basin of attraction theory discussed above. We finally note that line reclosing is technically feasible around $P/{B}_{0}\approx 0.5$ where the angle difference between the two ends of the tripped line is small (modulo $2\pi$). This would lead to a q = 1 vortex flow state.

Zoom In Zoom Out Reset image size

Dissipation renders operational conditions different for different stationary states. In particular, states with vortex flows generically have larger ohmic losses, because they have larger angle differences in equation (4). They therefore require an additional power injection ${\rm{\Delta }}{P}_{q}$ depending on the vorticity q. One would think that changing the operational conditions by e.g. reducing ${\rm{\Delta }}{P}_{q\ne 0}\to {\rm{\Delta }}{P}_{q=0}$ makes the q-vortex flow disappear. Figures 5(d) and (e) show however that topological protection by the winding number remains active, despite the presence of dissipation. In panel (d) we decrease the additional power injection ${\rm{\Delta }}P$ to the amount of losses incurred in the q = 0 stationary state when the system is in the q = 1 state, while in panel (e) we increase ${\rm{\Delta }}P$ to its value in the q = 1 state, starting from the q = 0 state, at $P/{B}_{0}=1$. In both cases, the winding number remains the same. Synchronization is furthermore not destroyed, however angles rotate at a modified frequency, ${\theta }_{i}(t)={\theta }_{i}^{(0)}\pm \delta {\rm{\Omega }}t$ in the rotating frame, with $\delta {\rm{\Omega }}\simeq {N}^{-1}({\rm{\Delta }}{P}_{q=1}-{\rm{\Delta }}{P}_{q=0})$. Adapting ${\rm{\Delta }}P$ to what is required by another q-state changes the synchronous frequency but leaves q unchanged.

We finally export the knowledge obtained from investigating simple systems to a network model with the topology of the UK high voltage AC power grid [19, 31]. Figure 6 illustrates vortex flow creation, enhanced ohmic losses and topological protection of stationary states with vortex flows. The system is initially stabilized in a stationary state without vorticity on any of its loops. It is later perturbed by a line tripping, at the position indicated in figure 6(e). The state is then left to stabilize towards a new state, after which it is again perturbed by the reclosing of the line. Finally, the power injected is modified to try and move the system back to the stationary state with q = 0, without success. Figures 6(a) and (b) show the behavior of the voltage angles and angle differences between connected nodes. It is seen that line tripping at t = 1600 essentially makes a single voltage angle difference significantly change, similarly to the single-loop model considered in figure 5. The same angle difference is the only one to move sensibly upon line reclosing at t = 5410. Figure 6(c) shows that the line fault creates a vortex flow. The latter is affected by adapting the power at t = 6970 to the losses of the initial state with q = 0 only insofar as all its angles start to rotate at a modified frequency $\delta {\rm{\Omega }}\simeq {N}^{-1}({\rm{\Delta }}{P}_{q=0}-{\rm{\Delta }}{P}_{q=-1})$. However this does not affect its winding number—topological protection is at work also in this case of a complex meshed grid with significant ohmic dissipation, ${G}_{{lm}}/{B}_{{lm}}=0.1$. Figure 6(d) shows that the vortex flow doubles ohmic losses, despite the fact that it affects only a small fraction of the grid, as is seen in figure 6(e) which shows differences in flows between the $q=-1$ and q = 0 states. The reduction in losses in figure 6(d) after t = 6970 leads to the synchronization of the grid at a frequency different from the rated frequency of 50 Hz. Such a change in frequency would be intolerable in a real power grid, and would quickly lead to either controlled or uncontrolled line trippings, cascades of failures, possibly leading to blackouts [32–34].

Zoom In Zoom Out Reset image size

Out of the three mechanisms for creating vortex flows we discussed, the topological phase twist mechanisms (ii) and (iii) are relevant to AC power grid operation as they can occur at relatively small voltage angle differences between connected nodes. The reliability of AC power grids is constantly evaluated via $N-1$ feasibility and transient stability analysis, where the existence of a stationary solution as well as the convergence towards that solution is checked after any one of the N major components (lines, transformers etc.) is removed from the network. We believe that this analysis should be complemented by checks of the presence of vortex flows, since the tripping or the reclosing of a line has the potential to generate them, resulting in reduced stability and higher, persistent ohmic losses, that are very difficult to get rid of.

The operating conditions of power grids are expected to change drastically as the energy transition steadily substitutes smaller, delocalized production for large power plants. As but one consequence, power generators have less and less mechanical inertia, thus less primary power reserve. To compensate for these changes, power electronics devices, phase angle regulating transformers and other devices whose task it is to effectively modify admittances and voltage phase angles are often incorporated into the power grid. The two topological phase twist mechanisms discussed above are in a way extreme, in that they rely on line tripping or reclosing. The conditions under which less stringent actions such as reducing line admittances would generate vortex flows should be investigated. Work along those lines is in progress and preliminary results seem to indicate that changes in power grids brought about by the energy transition have the potential to generate vortex flows more frequently.

We also discussed vortex flow creation via dynamical phase slip, despite its lack of relevance for AC power grids. High voltage power transmission is however closely connected to problems of coupled oscillators via the celebrated Kuramoto model [35–38]. Connections between the latter model and Josephson junction arrays were noted in [39]. Stationary states in systems of coupled oscillators with different winding numbers were discussed in [40, 41] and in a biological context in [42]. Our theory can be exported to those situations to explain how such states are created in the first place and how they disappear.

We finally comment on further analogies with vortex physics in superconductors. First, the point has already been made above that there is no counterpart to external magnetic fields in electrical grids that could generate vortices, and at present it is not known to us whether specific sequences of injection/consumption changes could lead to the creation of vortex flows. But if such a sequence exists, figure 4 shows that it should move the system over an 'energy' barrier, i.e. over a ridge or a saddle point of the Lyapunov function. This, in a sense, is to be related to vortex formation in the Ginzburg-Landau model for type II superconductors where an energy barrier is passed at a critical magnetic field beyond which the vortex state is energetically favorable. Second, vortex creation in superfluids occurs either via creation of a pair of vortex-antivortex or via vortex nucleation at the boundary of the system. Electrical power grids being rather small, we have found that vortex nucleation occurs at the network boundary, but would need to perform numerically intensive investigations on much bigger networks (the North American or the Pan-European networks for instance) to see if/when vortex-antivortex pairs are created. While we cannot rule out single vortex creation in the bulk from topological line modifications such as line tripping or reclosing, we suspect this occurs only very rarely, if at all, as such a process would require winding numbers to change simultaneously for a macroscopic number of loops encircling the vortex. Finally, it is known that vortices in superconductors move in the presence of a transport current, which leads to dissipation. So far we have seen vortex flows move only in numerical investigations on regular lattices. That we never saw it in complex networks may be due to either the finiteness of the system size considered, or to the absence of large loops neighboring those where vortex flows sit or to vortex pinning in these complex, effectively disordered networks. With our current knowledge, whether vortex flows can move around in large electric power grids is an open question.

This work has been supported by the Swiss National Science Foundation under an AP Energy Grant.