Multistability and anomalies in oscillator models of lossy power grids

The analysis of dissipatively coupled oscillators is challenging and highly relevant in power grids. Standard mathematical methods are not applicable, due to the lack of network symmetry induced by dissipative couplings. Here we demonstrate a close correspondence between stable synchronous states in dissipatively coupled oscillators, and the winding partition of their state space, a geometric notion induced by the network topology. Leveraging this winding partition, we accompany this article with an algorithms to compute all synchronous solutions of complex networks of dissipatively coupled oscillators. These geometric and computational tools allow us to identify anomalous behaviors of lossy networked systems. Counterintuitively, we show that loop flows and dissipation can increase the system’s transfer capacity, and that dissipation can promote multistability. We apply our geometric framework to compute power flows on the IEEE RTS-96 test system, where we identify two high voltage solutions with distinct loop flows.

Synchronization and flow networks.The history of scientific investigation about synchronization is traditionally traced back to Huygens' observation of an "odd kind of sympathy" in the XVIIth century. 1 It is however only in the last decades that a tractable framework has been developed, [2][3][4] thanks in particular to the pioneering works of Winfree in the 1960s, 5 and Kuramoto in the 1970s-80s. 6,7Shortly thereafter, the problem of synchronization has been embedded in the framework of network systems, 8-10 first based mostly on numerical simulations, evolving progressively towards more and more analytical results. 2,4,11Even in the simplest form of coupled oscillators, the interplay between dynamics and network structures leads to rich and sometimes unexpected behaviors.
The interactions between synchronizing oscillators is naturally interpreted as a flow of information or commodity between the nodes of a network.This dual interpretation of synchronization and flows is predominant in the modeling of voltage dynamics in high voltage power grids. 12,13Indeed, in power systems, a rotating turbine in a plant accumulates kinetic energy and accelerate if all the power it produces is not transmitted to its neighboring buses.While there is a natural link between synchronization and flow balance in power grids, a similar duality underlies dynamical systems as diverse as springconnected rotating masses, 11 motion planning, 14,15 or chemical oscillations 7 to name but a few.The rate of change of an oscillator's state is then determined by the imbalance of the flows received from or sent to its neighbors.When the flows of commodity balance out at each agent, such that all agents have identical rate of change, then the relative positions of the agents are constant in time: we say that they are synchronized.In particular, this is the desired state for an AC power grid.
Lossless oscillator networks.One of the simplest models of synchronization considers a set of oscillators, each described by a phase θ ∈ S 1 [−π, π), interacting with each other through a 2π-periodic coupling, function of their phase difference, for i ∈ {1, ..., n}.The real parameters m i and d i represent the effective inertia and damping of oscillator i respectively, and ω i is the natural frequencies that the oscillator would hold without interactions.The function f ij is the coupling between nodes i and j and a ij ∈ R ≥0 is the edge weight that scales the strength of the coupling and determines the underlying network structure.These coefficients are nonzero if and only if oscillators i and j interact.The system described by equation (1) evolves on the n-torus T n = (S 1 ) n .The oscillators are frequency synchronized (or phase locked) if, at some point in time, θi ≡ ϕ, for all i.The intrinsic compact nature of each oscillator's domain and the continuity of the coupling function require f ij to be nonlinear, and periodic in the systems of interest here.Whereas linear networked systems are well-understood, 16 the nonlinear nature of the coupling between oscillators can lead to rich and intricate behaviors. 13,17he vast majority of the literature about synchronization of coupled oscillators assumes symmetric couplings, i.e., two coupled oscillators influence each other with the same strength.In the flow network interpretation, symmetric couplings correspond to lossless flows, i.e., the flow of commodity between i and j contributes with equal magnitude and opposite sign to each end of the edge between i and j.Shortly put in mathematical terms, f ij (x) = −f ji (−x).Consequences of this strong relation between f ij and f ji are in particular: (i) conservation of the total flow in the system, simplifying the calculation of the asymptotics of equation (1) and; (ii) symmetry of the Jacobian matrix of the system, guaranteeing nice and convenient spectral features.The properties of systems with lossless couplings allowed to derive a long list of results about their dynamics: conditions for existence and uniqueness of their synchronous states; 11,18,19 multistability; 13,20 and clustering 21,22 to name but a few.An approach, common to various works, is to design a fixedpoint iteration, [18][19][20] whose convergence is guaranteed under some convexity properties of the energy landscape of the system. 23,24hallenges in the lossy oscillator systems.While the lossless assumption is reasonable in many cases, it is often not realistic and can lead to inaccurate predictions (see the power flow problem in the Results Section).In the flow interpretation, the transfer of a commodity, say electric power, is subject to dissipation, e.g., due to line resistance, meaning that the amount sent from i to j is strictly larger than the amount received by j from i.In mathematical terms, dissipation are introduced in equation (1) by adding a term to the coupling function satisfying g ij (x) = g ji (−x).Note that any pair of coupling functions (from i to j and from j to i) can be decomposed as the sum of f ij +g ij with the above properties [see equations (38) and (39) in the Methods Section].
The importance of understanding the more realistic case of dissipative couplings motivated the early work by Sakaguchi and Kuramoto 25,26 and is still an active field of research.Recent numerical investigations 27,28 as well as analytical studies in regular systems [29][30][31] are beginning to shed light on a more in-depth understanding of dissipative networks.][34] Up to this day, it is unclear to what extent the properties enjoyed by lossless networks are preserved in more realistic, dissipative systems.In the global scientific aim of a faithful modeling of real systems, it is of utmost importance to decipher the impact of dissipation in standard models of networked dynamics.Indeed, conditions for existence, uniqueness, and multiplicity of synchronous states or for the emergence of clustering in lossless systems 11,13,[20][21][22] are yet to be adapted to their dissipative counterpart.Furthermore, it is now largely documented 35 that phase frustration can lead to the occurrence of solitary and chimera states, 28,[36][37][38][39][40] that are extensively studied, but still only partially understood.
Understanding dissipative systems is challenging for a number of reasons.In such systems, flow conservation is lost and the linearization of the system typically loses its symmetry.Furthermore, while equation ( 1) can be formalized as a gradient system over an energy landscape, this property immediately fails in dissipative systems such as equation (2).Therefore, technical approaches based upon energy landscapes are not applicable any longer.Incorporating dissipation in the system even requires to re-think the intuitive vectorial formulation of equation (1), in order to recognize the directionality of flows.Notice that, surprisingly, even a clear vectorial form of the dissipative dynamics is lacking in the literature.
Models of lossy power grids.In the context of power grids, taking in equation ( 1) yields precisely the lossless approximation of the swing equations, 12,13 describing the time evolution of voltage phase angles, with ω i being the power injection or consumption at bus i.In normal operation, it is desired that the system is maintained in the vicinity of a synchronous state of the swing equations, which is reached at the solution of the lossless power flow equations (see the Methods Section).
As discussed above, while the lossless approximation can be fair and useful in power grids modeling, it is never exact.Therefore, an accurate mathematical analysis of voltage dynamics and power flows requires to take dissipation into consideration.Considering resistive losses in the swing equations boils down to taking (see Methods Section) The mathematical challenges encountered when relaxing the lossless line assumption confined most of the literature to numerical investigations. 27,28bjectives and contributions.The overarching goal of this article is to develop an analysis framework to characterize the location, properties, and stability of synchronous solutions of dissipative oscillator networks.In the task of globally characterizing synchronous states of lossless oscillator networks, an instructive and effective approach has been to leverage the concepts of winding numbers and winding cells. 20Given a cycle of oscillators σ = (θ 1 , ..., θ |σ| , θ 1 ), the associated winding number q σ (θ) counts the number of times the oscillators' angles wrap around the origin when following σ (a rigorous definition is given in the Results Section).A winding cell is a subset of the n-torus T n whose points share the same winding number around each cycle.Winding cells form a partition of the n-torus (the winding partition) and directly result from the network structure of the system.The concepts of winding numbers and winding cells are illustrated in Fig. 1.
In this article, we rigorously draw the link between the winding cells and the occurrence of a series of surprising FIG. 1. Winding cell partition and equilibria for a 3node spring network.Projected winding cell partition of the 3-torus (left) and representation of three equilibria of a spring network (right).Points on the 3-torus are projected on the 2-dimensional space of angular differences θ2 −θ1 and θ3 − θ2.A winding number q ∈ Z is associated to each equilibrium, counting the number of times its angles wind counterclockwise around the origin.The three dots (left) are equilibria of the spring network, labeled with their winding numbers q.The phase synchronous state has all angles identical and therefore q = 0. Equilibria with q = ±1 are the so-called splay states or twisted states.The colored areas in the left panel represent the set of points with same winding number, i.e., the winding cells, forming a partition of the 3-torus.
behaviors of dissipative oscillator systems, that have escaped analysis up to this day.Specifically, we show that, in some cases, increased dissipation can lead to more robust and more stable systems.We argue that the winding partition provides a clear phase portrait for the analysis of such behaviors.
Motivated by these first observations, we proceed to the second contribution of this article.Namely, we provide an analytical, statistical, and computational understanding of the solutions of dissipative flow problems.In particular, we show that, exactly as in lossless networks, there is at most one solution of the dissipative flow problem in each winding cell.This at most uniqueness property is rigorously proven for a small amount of dissipation and verified numerically for a wide range thereof.For acyclic networks we provide an algorithm computing the unique solution, if it exists.For general, cyclic graphs, we provide an iteration map that, under technical assumptions, converges to the unique solution in a given winding cell.An implementation of both algorithms is provided freely online. 41e conclude by providing a compelling example of a realistic power grid where multiple solutions of the power flows are accurately captured by the distinct winding cells.Namely, we find two power flow solutions on the IEEE RTS-96 test system, 42 belonging to two different winding cells.
Overall, in this article, we illustrate our findings with the Kuramoto-Sakaguchi model (see the Methods Section).This model is particularly appealing in our con-text: it is a natural extension of the (lossless) Kuramoto model; it is a first order version of the lossy swing equations; and the amount of dissipation in the coupling can easily be tuned with a continuous parameter, namely the phase frustration φ ∈ R. Nevertheless, our analytical results are valid for a much broader class of coupling functions f ij + g ij , that will be of interest to the dynamical systems and network science communities.
Remark.In addition to the demonstrations provided, the framework proposed here is naturally suited to the analytical study of networked dynamical system with directed interactions.For sake of conciseness and clarity, we limit our focus to dissipative interactions over undirected edges, but the framework covers naturally any type of directed interactions.We discuss these generalizations to a greater extent in the Discussion Section.

RESULTS
After a formal definition of the winding partition of the n-torus, we provide a careful description of a series of unexpected behaviors of dissipative oscillator networks.This section culminates with a presentation of our rigorous mathematical results.We conclude by giving an example of two solutions of the power flow equations, coexisting on the IEEE RTS-96 test system.A detailed formalism can be found in the Methods Section and proofs are deferred to the Supplementary Information.
Algebraic graph theory on the torus.Our framework is inspired by Ref. 20.The states of the system of equation ( 2) are points θ in the n-torus T n , each component being a point θ i of the circle S 1 .Comparing points on S 1 requires to define angular differences, which is somewhat arbitrary.In this article, we use the counterclockwise difference Intuitively, the counterclockwise difference is a projection of the angular difference on the interval [−π, π).
Given a cycle σ = (i 1 , ..., i |σ| , i 1 ) in a graph G u (see Methods for details), one can calculate the winding number around cycle σ associated with the state θ ∈ T n , Three states with different winding numbers are illustrated in Fig. 1 for the 3-cycle.Intuitively, the winding number counts the number of times the angles in θ wind around the origin when following the cycle σ.Then, given a cycle basis Σ = {σ 1 , ..., σ c } of the graph, we naturally define the winding vector associated to a state θ ∈ T n , a The transparent volume is Ω(+1; σ), the winding cell of winding number q = +1, and the solid volume is the 3π/4-cohesive set, i.e., the subset of Ω(+1; σ) where the counterclockwise differences do not exceed 3π/4.b The transparent volume is Ω(0; σ), the winding cell of winding number q = 0, and the solid volume is the π/2-cohesive set.c The transparent volume is Ω(−1; σ), the winding cell of winding number q = −1, and the solid volume is the 3π/4-cohesive set.d Union of the cohesive sets of the previous panels.Each color corresponds to a winding cell.A key result of this article is that there is at most one solution to equation ( 2) in a certain cohesive subset of each winding cell, i.e., in each solid volume.
Nonzero winding numbers are typically associated to loop flows, 20,43,44 i.e., a commodity flow of constant magnitude around a cycle of the network.Such loop flows occupy line capacity, but do not deliver commodity anywhere.
Remark.The winding number is a natural extension to complex networks, of the quantification of vortex flows in regular lattices, that arise in statistical physics [e.g., superfluids 45 or superconductors 46 ].As far as we can tell, the notion of winding numbers in systems of coupled oscillators can be traced back to Refs.47 (referee discussion) and 48.
For a graph with c cycles, a winding vector u ∈ Z c can be uniquely associated to each state in T n .Therefore we can define the winding cell associated with winding vector u, The counterclockwise difference is bounded, and so are the winding numbers.There is then a finite number of winding cells for a given graph G u , forming a finite partition of T n .See Fig. 2 for an illustration of winding cells in a cycle of n = 3 oscillators.Anomalous behaviors of dissipative systems.Three unexpected behaviors of lossy oscillator networks are illustrated in Fig. 3 for the Kuramoto-Sakaguchi model.The first one is a direct extension of a phenomenon already noted for lossless systems. 20,49The two other have not been reported to the best of our knowledge.
Anomaly 1, loop flows increase capacity: One would expect the presence of a loop flow (i.e., nonzero winding number) to reduce the transmission capacity of the system, because lines are occupied by the aformentioned loop flow.In Figs.3c and 3f, we see that the solutions with larger winding numbers tolerate larger commodity transfer.Even though such observations have been documented in the past for lossless systems, 20,49 it remains somewhat counterintuitive.
Anomaly 2, dissipation increases capacity: An initial reasoning would suggest that increasing dissipation would reduce the robustness and reliability of a system.Indeed, if part of the transmitted commodity is lost on the way, then more of it needs to be injected and the system is operated closer to criticality.However, the relation between dissipation and robustness is not that simple, as we illustrate in Fig. 3c.Indeed, for a nonzero winding number, the ability of the system to synchronize can evolve non-monotonously with respect to the dissipation (see solution at q = −1).Such phenomenon is quite unexpected and, to the best of our knowledge, has not been reported so far.Anomaly 3, dissipation promotes multistability: Different solutions differ by a collection of loop flows, 43 i.e., for some solutions, the lines are more loaded than for others.Similarly as in the previous anomaly, one would expect that increased dissipation would prevent the occurrence of loop flows and therefore of multiple solutions.However, according to Fig. 3f, a system with low dissipation (φ ∈ [0, 0.3]) and low injection (p ≈ 0) can have fewer solutions than more loaded and dissipative systems.Indeed, one would assume that lower loads and lower frustration leads to a larger margin of freedom in the system.Apparently, this is not necessarily the case and this can be attributed to the underlying network structure.
The anomalous behaviors identified above are typically related to the coexistence of different solutions.As we show in this article, there is a strong and direct link be- Solutions at: q = (0,0) q = (-1,+1) q = (+1,-1) 1 sol.Angles corresponding to the two solutions of panel a.One clearly sees that, for the solution at q = +1, the angles wrap around the circle, but not in the solution at q = 0. c Boundaries (colored curves) of the existence regions for solutions at different winding numbers, in the parameter space of phase frustration φ and injection magnitude p.The solutions in panels a and b were obtained for (φ, p) = (0.3, 1.0).The darkness of each area in the parameter space represents the number of existing synchronous states.It is surprising that (i) the solution at q = +1, i.e., with larger winding number, can carry a larger flow than the solution at q = 0, and (ii) for the solution at q = −1, the maximal tolerated commodity injection is not monotone in the frustration.d, e, f: Same (a, b, c) respectively, for the 2-cycle network in panel d.f Surprisingly, this network has fewer solutions for light load and frustration, (φ, p) ≈ (0.0, 0.0), rather than for larger parameter values, e.g., (φ, p) = (0.3, 1.0).
tween different solutions and the winding partition of the n-torus.
Problem setup and solution: Synchronous states with dissipative couplings.We now formalize the problem of flow distribution in dissipative networks and present our main formal results.We provide a summary of the main notation symbols in the Methods Section (Table I).
Let G u be the undirected graph describing the interactions in equation (2).Each edge e = {i, j} of G u is endowed with two coupling functions, one for each orientation.Without loss of generality, we incorporate the edge weight a ij in the coupling functions f ij , g ij , and h ij .For each edge, we choose an arbitrary orientation, say (i, j).We will refer to the coupling functions as h e = h ij and h ē = h ji , with ē denoting edge e with reversed orientation.
In our framework, equation ( 2) can be written in vectorial form as using the diagonal inertia and damping matrices M and D, the incidence matrix B u , and the out-incidence matrix B o induced by the chosen orientation for the incidence matrix [Methods Section, equations ( 21) and ( 25)].The coupling function h : R m → R 2m relates a vector of angular differences over the m undirected edges to the flows that are distinct for each edge orientation, hence h has The edge indices e ∈ {1, ..., m} follow the orientation induced by the incidence matrix B u .We refer to the discussion about the Kuramoto-Sakaguchi model in the Methods Section for an instructive example of the construction of equation (10).
From now on, it will be convenient to formulate the problem in terms of angular difference variables ∆ ∈ R m , rather than in terms of angle variables θ ∈ R n .Constructing the vector of angular differences ∆ from a vector of angles θ is straightforward, using the transpose of the incidence matrix, ∆ = B u θ.The other direction however is not that direct.Indeed, from a difference vector ∆, one can recover the associated angle vector θ over a spanning tree of the graph.Now, the angular difference vector is consistent with the graph structure only if some cycle constraints.Namely, over the remaining edges of the graph e = {i, j}, that are not in the spanning tree, the constraint is The integer multiple of 2π does not matter because the angles are compact variables over S 1 .Mathematically speaking, these cycle constraints can be formalized using the cycle-edge incidence matrix C Σ associated to a cycle basis Σ = (σ 1 , ..., σ c ) [formally defined in the Methods Section, equation (23)] for some winding vector u ∈ Z c .In summary, an angular difference vector ∆ needs to satisfy equation ( 12) in order to be consistent with an angle vector θ.
From now on, we will search for stable synchronous states of equation (10), i.e., we require the Jacobian matrix of the system to have its spectrum in the left complex half-plane. 50Following Gershgorin Circles Theorem, 51 stability is guaranteed if h e (∆ e ), h ē(−∆ e ) > 0 for all e.Formally, we assume that, in a neighborhood of the origin, both h e and h ē are strictly increasing, and for each edge e of G u , we require |∆ e | ≤ γ e , so derivatives are positive.The vector of angular differences is then restricted to the hypercube The set of points θ ∈ T n whose angular differences along the edges of G u are in R(γ) is referred to as a γ-cohesive set.The solid volumes in Fig. 2 show the intersections of the various winding cells of a 3-cycle and R(γ) for different values of γ e .
Gathering the above observations, we formulate the following problem, whose solutions are in one-to-one correspondence with synchronous states of equation (10).
Problem statement (Dissipative Flow Network).Given a connected graph G u with n nodes, m edges, and cycle basis Σ, a vector of natural frequencies ω ∈ R n , and appropriate coupling functions h e , h ē, associated to each edge e, find a solution ∆ ∈ R(γ) of for some synchronous frequency ϕ ∈ R and winding vector u ∈ Z m−n+1 .
In contrast with previous works on lossless systems, the flow map h e is not odd, meaning that we do not impose the constraint h e (θ i − θ j ) = −h ē(θ j − θ i ), hence our need of the out-incidence matrix B o instead in equation (14a).Note that, even though in our example of the Kuramoto-Sakaguchi model all coupling functions are identical, in full generality, we allow h e = h ē.
In the Methods Section, we provide a series of rigorous results, proving the following claim.
Claim.There is at most one solution to the Dissipative Flow Network problem in each winding cell of the ntorus.
More precisely, we distinguish the cases of acyclic graphs and of general graphs with cycles.For acyclic graphs, the winding partition is trivial and the whole ntorus is a single winding cell.By proving that there is at most one solution to the Dissipative Flow Network problem for such graphs, the claim is proven (see Theorem 3).
In the case of cyclic graph, we design a iterative map that we prove to be contracting within winding cells, under reasonable conditions (see Corollary 6).Solutions of the Dissipative Flow Network problem are fixed points of this iteration map and therefore, there is at most one solution in each winding cell.
The formal results are formulated in the Methods Section and proof are deferred to the Supplementary Information.
Multistability in power grids.The last decades have seen a large scale effort of the complex systems community to provide an analytical description of the power flow equations and of their solutions (see the Methods Section).In 1972 already, Korsak 47 showed that, mathematically speaking, the power flow equations tolerate multiple solutions on cyclic networks.3][54][55] Even some "realworld" events advocate in this direction. 56,57However, a large proportion of the work mentioned above relies on the lossless line assumption, namely, neglecting dissipation, voltage amplitude dynamics, and reactive power flows.Recently, there has been a common effort in trying to pursue a more realistic mathematical analysis of power grids, by incorporating reactive power flows, 18,19 voltage amplitude dynamics, 58,59 and dissipation. 27,28In particular, the recent linear stability analysis of the extendend swing equations proposed in Ref. 60 show that the conditions for the local stability of a syncrhonous state in lossy systems are very similar to those for lossless sytems.Despite all this work, there is still neither a clear extension of the winding partition to the full active-reactive power flows, nor a global phase portrait for lossy oscillator systems, even though there are some notable related preliminary works. 61,62Our results are an advance in the aforementioned collective effort.
To put our results in perspective with the resolution of the power flow equations, we solved both the Dissipative Flow Network Problem and the power flow equations on an adapted version of the IEEE RTS-96 test case. 41,42In Fig. 4, we compare synchronous states of the Kuramoto-Sakaguchi model (panels b and c, inner circle), with the corresponding solutions of the full power flow equations (panels b and c, outer annulus).We elaborate on the resolution of the full power flow equations in the Methods Section.First of all, one clearly sees that the main qualitative features (e.g., winding number, cohesiveness, clustering) of the power flow solutions are captured by the corresponding synchronous states of the Kuramoto-Sakaguchi model.Furthermore, it is remarkable that two solutions to the full power flow equations coexist, satisfying all voltage amplitude constraints as well as voltage angle stability.This example shows that loop flows and the winding partition are fundamental features of power flow solutions.Our work is a contribution to the joint and long lasting effort in the quest of an accurate mathematical analysis of power grids, which is a landmark in the area of power grids analysis.

DISCUSSION
Theorem 3 and Corollary 6 (see Methods Section) rigorously prove that, in each winding cell of the n-torus, there is at most a unique synchronous solution for dissipative networks of oscillators.In acyclic networks, the whole n-torus is trivially the unique winding cell, and there is therefore at most one solution to the Dissipative Flow Network problem (Theorem 3), independently of the amount of dissipation.For systems over more general, cyclic graphs, the winding partition provides a natural decomposition of the n-torus in subsets containing at most one solution.These results are a straight generalization of Ref. 20 to dissipative systems.
Even though the relation established in Corollary 6 is formally valid for relatively small amounts of dissipation, numerical experiments did not lead to any counterexample.Indeed, we empirically observed for a large range of network structures, frustration parameters, and initial conditions, that the iteration map defined in Theorem 5 [Methods Section, equation (37)] can always be made contracting by taking a sufficiently small value of > 0. We have therefore strong numerical evidence that the aforementioned iteration map can be made contracting at any point of the γ-cohesive set, with γ e taken such that each coupling function h e is strictly increasing on [−γ e , γ e ].We conjecture that Corollary 6 is actually very conservative in general and that the at most uniqueness property therein is valid for a much broader range of dissipation-to-coupling ratio.Furthermore, the comparison between coupling and dissipation in equation (45)  clearly pinpoints how dissipation works against synchronization.When considering lossless couplings, one realizes that the inequality in ( 45) is always satisfied, recovering the results of Ref. 20.
Both the proofs of Theorem 3 and Corollary 6 are algorithmic by nature.Namely, the proof of Theorem 3 considers recursively the flows on the edges of the acyclic graph, and Corollary 6 relies on an iteration map [S , equation (37)].It is therefore straightforward to actually implement the proofs as routines, which we provide online. 41he anomalies illustrated in Fig. 3 emphasize that the introduction of dissipation in the coupling between oscillators has a nontrivial and surprising impact of the dynamics.The fact that both loop flows and dissipation can increase the transmission capacity of a system (Anomalies 1 and 2) is arguably counterintuitive.We remark that both Anomalies 1 and 2 occur for solutions with nonzero winding numbers (q = +1 and q = −1 respectively).It is also quite unexpected that more loaded and dissipative systems can possess more flow network solutions for a given network structure (Anomaly 3).Again, this last anomaly involves solutions in different winding cells.All anomalies identified in Fig. 3 are strongly linked to solutions with nontrivial winding numbers.A general and thorough description of the different operating states of dissipative networks of oscillators then requires to tackle these systems through the prism of the winding partition.On top of that, through our realistic example on the IEEE RTS-96 test system, we show that the winding partition will be relevant in the analysis of multiple solutions to the full power flow equations.
We trust that the notion of winding partition has the potential to contribute elucidating many open problems in the fascinating phenomenon of synchronization in complex networks.We reiterate that even though we restricted our discussion to bidirected interactions for sake of clarity, the whole framework developed in this article naturally applies to any system with directed interactions.Namely, our formalism is a first step towards a unified analysis of synchronization in any network of coupled oscillators, no matter the nature of the interactions.

METHODS
We first provide the necessary grounds of directed and undirected graph theory, as well as a link between them.We point to Ref. 16 for an extended discussion about graph and digraph theory.We then discuss the Kuramoto-Sakaguchi model and its link with the power flow equations.We conclude this section with the formulation of our theoretical results.Directed graphs.A directed graph (or digraph) G d is the pair (V, E d ) composed of a set of vertices (or nodes) V = {1, ..., n} and a set of directed edges E d ⊂ V × V , which are ordered pairs of vertices.For an edge e = (i, j) ∈ E d , i is the source of e, denoted s e , and j is its target, denoted t e , i.e., e = (s e , t e ).We denote the edge with opposite direction as ē = (t e , s e ).The existence of edges is encoded in the graph's adjacency matrix The out-degrees (resp.in-degrees) matrix is obtained as For a digraph with n vertices and m directed edges, we define the n×m out-incidence and in-incidence matrices which form the standard incidence matrix We notice the following relations, the fourth being unknown as far as we can tell.
Proposition 1.The adjacency matrix A d , the out-and in-degree matrices D o and D i , and the Laplacian matrix L d of a directed graph can be written in terms of its outand in-incidence matrices B o and B i : Proof.The proofs for the adjacency matrix A d and for the degree matrices D o and D i can be found in Ref. 63 (Lemmas 3.1 and 4.1).The proof for the Laplacian matrix directly follows, Remark.The same proof is straightforwardly adapted to weighted directed graphs.
Undirected graphs.An (undirected) graph G u is a pair (V, E u ) composed of a set of n vertices (or nodes) V = {1, ..., n} and a set of m edges, which are unordered pairs of vertices, E u ⊂ {{i, j} : i, j ∈ V }.A cycle of G u is an ordered sequence of vertices σ = (i 0 , i 1 , ..., i = i 0 ), such that {i j , i j+1 } ∈ E u and i j = i k for any j, k ∈ {1, ..., }.
Let us now choose an arbitrary orientation [(i, j) or (j, i)] for each undirected edge {i, j} ∈ E u .We can then define the incidence matrix of G u , The Laplacian matrix of G u can be computed as L u = B u B u .Note that the incidence matrix is not unique and depends on the choice of edge orientations, whereas the Laplacian does not.Given a cycle σ, we define the cycle vector v σ ∈ {−1, 0, +1} m , indexed by edges, as The cycle space of G u is the span of the cycle vectors of all cycles of G u , which is equivalently defined as the kernel of the incidence matrix B u .A set of cycles Σ = {σ 1 , ..., σ c } is a cycle basis of G u if and only if the set of corresponding cycle vectors forms a basis of the cycle space.
Finally, given a cycle basis Σ of the graph G u , we define the cycle-edge incidence matrix, Bidirected graphs.Dissipative couplings intrinsically require to distinguish the two orientations of each edge.Given an undirected graph G u = (V, E u ), its bidirected counterpart is the directed graph G b = (V, E b ) with the same vertex set V and where each undirected edge {i, j} ∈ E u is doubled in the set of directed edges (i, j), (j, i) ∈ E b .A bidirected graph is a directed graph induced by an undirected graph.
If the undirected graph G u has incidence matrix B u ∈ R n×m [equation ( 21)], then, with an appropriate indexing of the directed edges, the incidence matrix of G b can be written as B b = (B u , −B u ) ∈ R n×2m .Interestingly, we note that the Laplacian matrices of G u and G b are the same, namely (see Prop. 1), where the out-incidence matrix Notice that here, The Kuramoto-Sakaguchi model.We illustrate the results of this article with the generalized Kuramoto-Sakaguchi model, 8,25,26 for i ∈ {1, ..., n}, where θ i ∈ S 1 and ω i ∈ R are respectively the phase angle and the natural frequency of the i-th oscillator, φ ij ∈ (−π/2, π/2) is the phase frustration between oscillators, and in this case, a ij ∈ R ≥0 is the coupling strength between oscillators i and j.The Kuramoto-Sakaguchi model directly translates to the framework of equation ( 2), with m i ≡ 0, It is a natural extension of the original Kuramoto model, which is recovered for φ ij = 0.The coupling function of the Kuramoto-Sakaguchi model is illustrated in Fig. 5a.
Remark.While in its original formulation, the Kuramoto-Sakaguchi model assumes homogeneous, all-to-all couplings, here we take the couplings to be given by an underlying network structure.For the sake of simplicity, in our examples we consider a ij = a ji and φ ij = φ ji for all connected nodes i and j.Nevertheless, we keep in mind that these assumptions are not necessary for our results and that our framework is appropriate for much more general cases.with φ = 0.5.The light green curve illustrates the coupling on the same edge, but with opposite orientation [hē(−x) = sin(−x − φ) + sin(φ)].The thick parts (cyan and green) emphasize the region where the curve is increasing (resp.decreasing).The shaded gray area shows the interval where the coupling in both orientations is strictly monotone.(b) Odd (cyan) and even (green) parts of the Kuramoto-Sakaguchi coupling function, as defined in equation (44).
In order to illustrate some fundamental complications that arise in the Kuramoto-Sakaguchi model, compared to the original Kuramoto model, we detail two simple examples below.
In the Kuramoto model (φ = 0), it is standard to write the dynamics in vectorial form as where B u ∈ R n×m is the incidence matrix of the (undirected) coupling graph of the system, and A ∈ R m is the diagonal matrix of the edge weights a ij .In the Kuramoto-Sakaguchi model, this vectorial form is not that simple.Direct computation shows that naively writing does not yields the desired equations (27).
In order to write equation (26) in vectorial form, we need to distinguish the two orientation of each edge and consider G b , the bidirected counterpart of G u .We need to introduce the out-incidence matrix B o and the incidence matrix B b of the bidirected coupling graph, as well as A b ∈ R 2m , defined to be a block diagonal matrix whose two diagonal blocks are equal to A. The proof of the following proposition follows from direct computation.
There are at least two main messages that can be taken from these examples.First, by extending our framework to directed graphs, we are able to write the Kuramoto-Sakaguchi model in vectorial from, in equations (30).Note that a similar vectorial formulation of the Kuramoto-Sakaguchi model has recently been proposed in Ref. 64, which, while more general than equation (30), does not provide as much insight in the underlying network structure.
Second, unlike the Kuramoto model, the average frequency of the system is not preserved along arbitrary trajectories.Also, if multiple synchronous states exist, then they have, in general, different synchronous frequencies.
These claims are backed up by showing that the average frequency of the system depends on angular differences, which is time varying over the trajectories of the system, and not identical for different synchronous states.
On top of that, we reiterate that, contrary to the Kuramoto model, the Kuramoto-Sakaguchi model is not the gradient of any function (even locally).Therefore, the energy landscape approaches, valid for φ = 0, 13,24 are not directly applicable when φ = 0.
The power flow equations.Under the assumption that voltage amplitudes are fixed, synchronous states of the Kuramoto-Sagauchi model are in direct correspondence with the solutions of the active power flow equations. 13,65The power flow equations relate the the balance of active (P i ) and reactive powers (Q i ) to the voltage amplitude (V i ) and phase (θ i ) at each node i ∈ {1, ..., n}, with G ij and B ij being lines conductance and susceptance respectively.Defining one verifies that solutions of equation (34) are steady states of equation (26).Equations ( 34) and ( 35) are usually solved by iterative methods.In Fig. 4b and c, outer annulus, we used a Newton-Raphson scheme 66 with different, carefully chosen initial conditions to solve the full power flow equations on our version of the IEEE RTS-96 test case. 42he squares are PV buses, the circles are PQ buses, and the slack bus is node 23.The synchronous states of the Kuramoto-Sakaguchi models were computed by the flow mismatch iteration S [equation (37)], with = 0.01.All data are available online. 41he Dissipative Flow Network Problem on acyclic graphs.In the case where G u is acyclic, we show that there is at most a unique solution to the Dissipative Flow Network Problem.Here there are obviously no cycle constraint and thus equation (14b) is trivially satisfied.
Theorem 3. Consider the Dissipative Flow Network Problem on a connected acyclic undirected graph G u .Then there is at most one ∆ ∈ R(γ) that satisfies equation (14a).
The proof of Theorem 3 proceeds recursively and we provide it in the Supplementary Information.An implementation of an algorithm deciding the existence of the unique solution is provided online. 41mark.Theorem 3 is the dissipative version of Theorem 2.2 in Ref. 20.The spirit of Theorem 3 is somewhat similar to Ref. 27, even though therein, the authors restrict their investigation to the Kuramoto-Sakaguchi model and cannot extend their approach to more general couplings.
The Dissipative Flow Network Problem on general graphs.The presence of cycles in the network can induce the existence of multiple solutions to the Dissipative Flow Network Problem [see Fig. 3 or Ref. 27].We rigorously show here that winding vectors characterize these solutions for sufficiently moderate dissipation.
To do so, we define the flow mismatch iteration S over the space of angular differences R m , whose fixed points are exactly the solutions of equation (14a).Namely, let where > 0 is a small step size and L † u is the pseudoinverse of the graph Laplacian matrix.The flow mismatch iteration S updates the vector of angular differences according to the mismatch between the input/output of commodities ω and the distribution of flows that corresponds to the current angular differences.It has two major properties: (I) the vector ∆ * ∈ R(γ) is a fixed point of S if and only if it is a solution of equation (14a); (II) the map S leaves each winding cell invariant, because C Σ B u = 0.It means that fixing the winding vector of the initial conditions, imposes the winding vector of the fixed point of S , if ever it exists.
One of the main lessons from Ref. 20 is that different solutions of the flow network problem on the n-torus are better understood when put in the context of their winding cell.Accordingly, and thanks to the property (II) above, we split the Dissipative Flow Network Problem in each winding cell of the n-torus induced by the network structure.Fixing a winding vector u ∈ Z m−n+1 , we are guaranteed that if the initial conditions ∆ 0 satisfy equation (14b), then each following iteration ∆ k+1 = S (∆ k ) will satisfy it as well.
We summarize the above observations in the following theorem, whose proof is a direct consequence of the compactness of R(γ).In what follows, we provide sufficient conditions for contractivity of S .We rely on the decomposition of the coupling functions as h e (x) = f e (x) + g e (x), which implies, One readily verifies the identities Remind that g e quantifies to what extent the coupling is dissipative.In the particular case where the coupling is lossless, then g e = 0. Equipped with this decomposition of the couplings, we define two state dependent matrices, for x ∈ R(γ): (a) the odd weighted Laplacian matrix, which is the Laplacian matrix of G u weighed by the derivatives of the odd parts We emphasize that the choice of orientation for each edge e does not matter in the definition of L f .Also, the graph G u being connected and the above weights being positive, it is standard result of algebraic graph theory that λ 2 , the smallest nonzero eigenvalue of L f (a.k.a., the algebraic connectivity), is positive; (b) the even weighted degree matrix, which is the diagonal matrix weighted by the absolute derivatives of the even parts, where E i is the set of (undirected) edges incident to node i.The diagonal terms of D g quantify the dissipativity of the couplings.In particular, for lossless couplings, D g = 0.
Example.In the case of the Kuramoto-Sakaguchi model, the coupling functions are and trigonometric identities yield which we illustrate in Fig. 5b.We clearly see here the relation between g e and the dissipativity or frustration of the coupling.When φ = 0, we recover the original Kuramoto model, where the coupling is lossless, and g e = 0.
We are now ready to formulate the main theorem of this work.It clearly separates the impact of network connectivity, that promote the contractivity of S , and of the dissipation, that works against contractivity of S .We defer the proof to the Supplementary Information.
Theorem 5. Given a Dissipative Flow Network Problem, define the odd weighted Laplacian L f and the even weighted degree matrix D g .If, for all i ∈ {1, ..., n}, then there exists a sufficiently small step size > 0 such that the flow mismatch iteration S [equation (37)] is contracting.
The left-hand-side of equation ( 45) quantifies the amount of dissipation that is "seen" at each node of the network, which vanishes for lossless couplings.The righthand-side accounts both for the strength of the coupling between each pair of oscillators, through the weights, and for the connectedness of the graph, λ 2 being the algebraic connectivity. 67Under our assumptions [G u is connected, couplings are strictly increasing on R(γ)], the right-handside of equation ( 45) is necessarily positive.For couplings with sufficiently low dissipation, equation ( 45) is then satisfied, which, combined with Theorem 4, yields the following corollary.Corollary 6.If equation ( 45) is satisfied, then there is at most a unique synchronous state of equation ( 2) in each winding cell.The number of synchronous states is then bounded from above by the number of winding cells.
Theorem 5 and Corollary 6 give a rigorous, even though conservative, sufficient condition for at most uniqueness of synchronous states in each winding cell.However, computing the eigenvalues of the odd weighted Laplacian can be time consuming.We therefore propose some lower bounds on λ 2 (L f ) in the Supplementary Information (Prop.S2) that are state-independent and may ease the verification of equation (45).The bounds are adapted from standard results of algebraic graph theory.(ii).We adapt here the proof of Ref. S4, Lemma 1.9.Let v be the eigenvector of L f associated with λ2 and assume that |vi| = max k |v k | (recall that all these quantities depend on x).Because L f is a symmetric Laplacian matrix, 1 v = 0 and there is an index j such that vivj < 0. We denote with Pij the shortest (weighted) path from i to j.Now, where the maximum is taken over all pairs of vertices and the minimum is taken over all simple paths joining i and j, we can apply the Sedrakyan inequality S6 (direct consequence of the Cauchy-Schwarz inequality) to the numerator of equation (S82) and get where dw is defined in Tab.S2.
(iii).Alternatively, both sides of the identity

FIG. 2 .
FIG.2.Winding cells and cohesive sets in a 3-node system.Winding cells and their cohesive subsets, for a cycle of length n = 3.The 3-dimensional plots shows the unfolded 3-torus, where each dimension parametrizes one of the three angles and the winding cells become polytopes.The sides of the cube have then to be considered as identified (left-right, top-bottom, front-back).a The transparent volume is Ω(+1; σ), the winding cell of winding number q = +1, and the solid volume is the 3π/4-cohesive set, i.e., the subset of Ω(+1; σ) where the counterclockwise differences do not exceed 3π/4.b The transparent volume is Ω(0; σ), the winding cell of winding number q = 0, and the solid volume is the π/2-cohesive set.c The transparent volume is Ω(−1; σ), the winding cell of winding number q = −1, and the solid volume is the 3π/4-cohesive set.d Union of the cohesive sets of the previous panels.Each color corresponds to a winding cell.A key result of this article is that there is at most one solution to equation (2) in a certain cohesive subset of each winding cell, i.e., in each solid volume.

FIG. 3 .
FIG. 3. Anomalous behaviors in the Kuramoto-Sakaguchi model.Illustration of the anomalous behaviors identified for the Kuramoto-Sakaguchi model on a cycle network (a, b, c) and a 2-cycle network (d, e, f), with unit coupling weights.a Qualitative distributions of flows over a cycle of n = 18 oscillators, with commodity injection +p (resp.withdrawal −p) at node 16 (resp.4).The arrows of two different colors visualize different flow solutions, with different winding numbers.bAngles corresponding to the two solutions of panel a.One clearly sees that, for the solution at q = +1, the angles wrap around the circle, but not in the solution at q = 0. c Boundaries (colored curves) of the existence regions for solutions at different winding numbers, in the parameter space of phase frustration φ and injection magnitude p.The solutions in panels a and b were obtained for (φ, p) = (0.3, 1.0).The darkness of each area in the parameter space represents the number of existing synchronous states.It is surprising that (i) the solution at q = +1, i.e., with larger winding number, can carry a larger flow than the solution at q = 0, and (ii) for the solution at q = −1, the maximal tolerated commodity injection is not monotone in the frustration.d, e, f: Same (a, b, c) respectively, for the 2-cycle network in panel d.f Surprisingly, this network has fewer solutions for light load and frustration, (φ, p) ≈ (0.0, 0.0), rather than for larger parameter values, e.g., (φ, p) = (0.3, 1.0).

FIG. 4 .
FIG.4.Multistability in the full power flow equations for the IEEE RTS-96 test case.Comparison of the power flow solutions and Kuramoto-Sakaguchi synchronous states on the IEEE RTS-96 test case.42(a) Geographic representation of the system.Circles are loads and squares are generators.The network is composed of n = 73 nodes, m = 108 edges, and therefore c = 36 independent cycles.The long cycle with thick edges is of particular interest, because its length promotes the existence of loop flows while keeping angular differences small [see Refs.43,44 for an extended discussion].(b, c) Combined representations of: (outer annulus) the complex voltages for solutions to the full power flow equations for an adapted version of the IEEE RTS-96 test case; (inner circle) the phase angles of synchronous states of the Kuramoto-Sakaguchi model on the same system.For sake of readability, only the values of the nodes around the long cycle of panel a are represented.The outer annulus represents the tolerated margin of variation for the voltage amplitudes in the power flow equations.The power flow solution in panel b has a nonzero winding number (q = +1) and there is a reasonable correspondence (ordering, clustering) between its voltage phases and the angles of the Kuramoto-Sakaguchi synchronous state.Similarly, both the power flow solution and the Kuramoto-Sakaguchi synchronous states in panel c have zero winding number, with all angles in a relatively short arc.

Theorem 4 .
If the flow mismatch iteration S is contracting, then there is at most one synchronous state of equation (2) in each winding cell.