# Multistability and anomalies in oscillator models of lossy power grids

Robin Delabays,1,∗ Saber Jafarpour,1,2 and Francesco Bullo1 1Center for
Control, Dynamical Systems, and Computation, UC Santa Barbara, Santa Barbara
CA 93106-5070, USA.
2School of Electrical and Computer Engineering, Georgia Institute of
Technology, Atlanta GA 30332-0250, USA.
∗Corresponding author (robindelabays@ucsb.edu).

###### Abstract

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; 7 Shortly thereafter, the problem of synchronization has
been embedded in the framework of network systems, 8; 9; 10 first based mostly
on numerical simulations, evolving progressively towards more and more
analytical results. 4; 2; 11 Even 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; 13 Indeed, 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 spring-connected
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
$\theta\in\mathbb{S}^{1}\simeq[-\pi,\pi)$, interacting with each other through
a $2\pi$-periodic coupling, function of their phase difference,

$\displaystyle m_{i}\ddot{\theta}_{i}+d_{i}\dot{\theta}_{i}$
$\displaystyle=\omega_{i}-\sum_{j=1}^{n}a_{ij}f_{ij}(\theta_{i}-\theta_{j})\,,$
(1)

for $i\in\\{1,...,n\\}$. The real parameters $m_{i}$ and $d_{i}$ represent the
effective inertia and damping of oscillator $i$ respectively, and $\omega_{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}\in\mathbb{R}_{\geq 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
$\mathbb{T}^{n}=(\mathbb{S}^{1})^{n}$. The oscillators are frequency
synchronized (or phase locked) if, at some point in time,
$\dot{\theta}_{i}\equiv\varphi$, 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; 17

The 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 fixed-point iteration, 18;
19; 20 whose convergence is guaranteed under some convexity properties of the
energy landscape of the system. 23; 24

Challenges 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

$\displaystyle m_{i}\ddot{\theta}_{i}+d_{i}\dot{\theta}_{i}$
$\displaystyle=\omega_{i}-\sum_{j=1}^{n}a_{ij}\left[f_{ij}(\theta_{i}-\theta_{j})+g_{ij}(\theta_{i}-\theta_{j})\right]\,,$
(2)

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 (50) and (51) 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. More
generally, the extension of standard approaches to more realistic systems is
gaining momentum in the fields of synchronization and complex networks. 32;
33; 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, 36; 37;
38; 39; 28; 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

$\displaystyle f_{ij}(x)$ $\displaystyle=\sin(x)\,,$ (3)

in equation (1) yields precisely the lossless approximation of the swing
equations, 12; 13 describing the time evolution of voltage phase angles, with
$\omega_{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)

$\displaystyle f_{ij}(x)$ $\displaystyle=\cos(\phi_{ij})\sin(x)\,,$
$\displaystyle g_{ij}(x)$ $\displaystyle=\sin(\phi_{ij})[1-\cos(x)]\,.$ (4)

The mathematical challenges encountered when relaxing the lossless line
assumption confined most of the literature to numerical investigations. 28; 27

Objectives 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. 20 Given a cycle of oscillators
$\sigma=(\theta_{1},...,\theta_{|\sigma|},\theta_{1})$, the associated winding
number $q_{\sigma}(\bm{\theta})$ counts the number of times the oscillators’
angles wrap around the origin when following $\sigma$ (a rigorous definition
is given in the Results Section). A winding cell is a subset of the $n$-torus
$\mathbb{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.

Figure 1: Winding cell partition and equilibria for a 3-node 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
$\theta_{2}-\theta_{1}$ and $\theta_{3}-\theta_{2}$. A winding number
$q\in\mathbb{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=\pm 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.

In this article, we rigorously draw the link between the winding cells and the
occurrence of a series of surprising 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. 41

We 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 context: 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 $\phi\in\mathbb{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.

##  I 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 $\bm{\theta}$ in the $n$-torus
$\mathbb{T}^{n}$, each component being a point $\theta_{i}$ of the circle
$\mathbb{S}^{1}$. Comparing points on $\mathbb{S}^{1}$ requires to define
angular differences, which is somewhat arbitrary. In this article, we use the
counterclockwise difference

$\displaystyle d_{\rm cc}(\theta_{1},\theta_{2})$ $\displaystyle={\rm
mod}(\theta_{1}-\theta_{2}+\pi,2\pi)-\pi\in[-\pi,\pi)\,.$ (5)

Intuitively, the counterclockwise difference is a projection of the angular
difference on the interval $[-\pi,\pi)$.

Figure 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 $\Omega(+1;\sigma)$, the winding cell of winding
number $q=+1$, and the solid volume is the $3\pi/4$-cohesive set, i.e., the
subset of $\Omega(+1;\sigma)$ where the counterclockwise differences do not
exceed $3\pi/4$. b The transparent volume is $\Omega(0;\sigma)$, the winding
cell of winding number $q=0$, and the solid volume is the $\pi/2$-cohesive
set. c The transparent volume is $\Omega(-1;\sigma)$, the winding cell of
winding number $q=-1$, and the solid volume is the $3\pi/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.

Given a cycle $\sigma=(i_{1},...,i_{|\sigma|},i_{1})$ in a graph ${G}_{\rm u}$
(see Methods for details), one can calculate the winding number around cycle
$\sigma$ associated with the state $\bm{\theta}\in\mathbb{T}^{n}$,

$\displaystyle q_{\sigma}(\bm{\theta})$
$\displaystyle=(2\pi)^{-1}\sum_{j=1}^{|\sigma|}d_{\rm
cc}(\theta_{i_{j}},\theta_{i_{j+1}})\in\mathbb{Z}\,.$ (6)

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 $\bm{\theta}$ wind around the origin when following the cycle $\sigma$.
Then, given a cycle basis $\Sigma=\\{\sigma_{1},...,\sigma_{c}\\}$ of the
graph, we naturally define the winding vector associated to a state
$\bm{\theta}\in\mathbb{T}^{n}$,

$\displaystyle{\bf q}_{\Sigma}(\bm{\theta})$
$\displaystyle=[q_{\sigma_{1}}(\bm{\theta}),...,q_{\sigma_{c}}(\bm{\theta})]^{\top}\in\mathbb{Z}^{c}\,.$
(7)

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 ${\bf u}\in\mathbb{Z}^{c}$ can
be uniquely associated to each state in $\mathbb{T}^{n}$. Therefore we can
define the winding cell associated with winding vector ${\bf u}$,

$\displaystyle\Omega({\bf u};\Sigma)$
$\displaystyle=\left\\{\bm{\theta}\in\mathbb{T}^{n}\colon{\bf
q}_{\Sigma}(\bm{\theta})={\bf u}\right\\}\,.$ (8)

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}_{\rm
u}$, forming a finite partition of $\mathbb{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. LABEL:fig:anomalies for the
Kuramoto-Sakaguchi model. The first one is a direct extension of a phenomenon
already noted for lossless systems. 20; 49 The 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. LABEL:fig:anomaliesc and LABEL:fig:anomaliesf, 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. LABEL:fig:anomaliesc. 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.
LABEL:fig:anomaliesf, a system with low dissipation ($\phi\in[0,0.3]$) and low
injection ($p\approx 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 between 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 1).

Let ${G}_{\rm u}$ be the undirected graph describing the interactions in
equation (2). Each edge $e=\\{i,j\\}$ of ${G}_{\rm u}$ is endowed with two
coupling functions,

$\displaystyle h_{ij}(x)$ $\displaystyle=f_{ij}(x)+g_{ij}(x)\,,$
$\displaystyle h_{ji}(x)$ $\displaystyle=f_{ji}(x)+g_{ji}(x)\,,$ (9)

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_{\bar{e}}=h_{ji}$, with
$\bar{e}$ denoting edge $e$ with reversed orientation.

In our framework, equation (2) can be written in vectorial form as

$\displaystyle M\ddot{\bm{\theta}}+D\dot{\bm{\theta}}$
$\displaystyle=\bm{\omega}-B_{\rm o}\left[{\bf f}(B_{\rm
u}^{\top}\bm{\theta})+{\bf g}(B_{\rm u}^{\top}\bm{\theta})\right]$
$\displaystyle=\bm{\omega}-B_{\rm o}{\bf h}(B_{\rm u}^{\top}\bm{\theta})\,,$
(10)

using the diagonal inertia and damping matrices $M$ and $D$, the incidence
matrix $B_{\rm u}$, and the out-incidence matrix $B_{\rm o}$ induced by the
chosen orientation for the incidence matrix [Methods Section, equations (30)
and (37)]. The coupling function ${\bf
h}\colon\mathbb{R}^{m}\to\mathbb{R}^{2m}$ relates a vector of angular
differences over the $m$ undirected edges to the flows that are distinct for
each edge orientation, hence ${\bf h}$ has $2m$ components,

$\displaystyle[{\bf h}({\bf y})]_{e}$ $\displaystyle=h_{e}(y_{e})\,,$
$\displaystyle[{\bf h}({\bf y})]_{e+m}$ $\displaystyle=h_{\bar{e}}(-y_{e})\,.$
(11)

The edge indices $e\in\\{1,...,m\\}$ follow the orientation induced by the
incidence matrix $B_{\rm 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 ${\bm{\Delta}}\in\mathbb{R}^{m}$, rather than in
terms of angle variables $\bm{\theta}\in\mathbb{R}^{n}$. Constructing the
vector of angular differences ${\bm{\Delta}}$ from a vector of angles
$\bm{\theta}$ is straightforward, using the transpose of the incidence matrix,
${\bm{\Delta}}=B_{\rm u}^{\top}\bm{\theta}$. The other direction however is
not that direct. Indeed, from a difference vector ${\bm{\Delta}}$, one can
recover the associated angle vector $\bm{\theta}$ 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
$\theta_{i}-\theta_{j}=\Delta_{e}+2\pi k$, $k\in\mathbb{Z}$. The integer
multiple of $2\pi$ does not matter because the angles are compact variables
over $\mathbb{S}^{1}$. Mathematically speaking, these cycle constraints can be
formalized using the cycle-edge incidence matrix $C_{\Sigma}$ associated to a
cycle basis $\Sigma=(\sigma_{1},...,\sigma_{c})$ [formally defined in the
Methods Section, equation (35)]

$\displaystyle C_{\Sigma}{\bm{\Delta}}$ $\displaystyle=2\pi{\bf u}\,,$ (12)

for some winding vector ${\bf u}\in\mathbb{Z}^{c}$. In summary, an angular
difference vector ${\bm{\Delta}}$ needs to satisfy equation (12) in order to
be consistent with an angle vector $\bm{\theta}$.

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. 50 Following Gershgorin Circles Theorem, 51 stability
is guaranteed if
$h_{e}^{\prime}(\Delta_{e}),h_{\bar{e}}^{\prime}(-\Delta_{e})>0$ for all $e$.
Formally, we assume that, in a neighborhood of the origin, both $h_{e}$ and
$h_{\bar{e}}$ are strictly increasing, and for each edge $e$ of ${G}_{\rm u}$,
we require $|\Delta_{e}|\leq\gamma_{e}$, so derivatives are positive. The
vector of angular differences is then restricted to the hypercube

$\displaystyle R({\bm{\gamma}})$ $\displaystyle=\bigcap_{e\in{E}_{\rm
u}}[-\gamma_{e},\gamma_{e}]\subset\mathbb{R}^{m}\,.$ (13)

The set of points $\bm{\theta}\in\mathbb{T}^{n}$ whose angular differences
along the edges of ${G}_{\rm u}$ are in $R({\bm{\gamma}})$ is referred to as a
$\bm{\gamma}$-cohesive set. The solid volumes in Fig. 2 show the intersections
of the various winding cells of a $3$-cycle and $R({\bm{\gamma}})$ for
different values of $\gamma_{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}_{\rm u}$ with $n$ nodes, $m$ edges, and cycle
basis $\Sigma$, a vector of natural frequencies
$\bm{\omega}\in\mathbb{R}^{n}$, and appropriate coupling functions
$h_{e},h_{\bar{e}}$, associated to each edge $e$, find a solution
${\bm{\Delta}}\in R({\bm{\gamma}})$ of

$\displaystyle B_{\rm o}{\bf
h}({\bm{\Delta}})-\bm{\omega}=\varphi\bm{1}_{n}\,,$ (14a) $\displaystyle
C_{\Sigma}{\bm{\Delta}}=2\pi{\bf u}\,,$ (14b)

for some synchronous frequency $\varphi\in\mathbb{R}$ and winding vector ${\bf
u}\in\mathbb{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}(\theta_{i}-\theta_{j})=-h_{\bar{e}}(\theta_{j}-\theta_{i})$, hence our
need of the out-incidence matrix $B_{\rm 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}\neq h_{\bar{e}}$.

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 $n$-torus.

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 $n$-torus 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. Since then, there
has been a plethora of evidence, both analytical and numerical, that the power
flow equations allow the coexistence of different solutions. 52; 53; 54; 13;
55 Even some ”real-world” events advocate in this direction. 56; 57 However, 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; 28 In 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; 62 Our 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. 42; 41 In
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.

Figure 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.

##  II 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 (49)] can always
be made contracting by taking a sufficiently small value of $\epsilon>0$. We
have therefore strong numerical evidence that the aforementioned iteration map
can be made contracting at any point of the $\bm{\gamma}$-cohesive set, with
$\gamma_{e}$ taken such that each coupling function $h_{e}$ is strictly
increasing on $[-\gamma_{e},\gamma_{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 (57)
clearly pinpoints how dissipation works against synchronization. When
considering lossless couplings, one realizes that the inequality in (57) 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_{\epsilon}$,
equation (49)]. It is therefore straightforward to actually implement the
proofs as routines, which we provide online. 41

The anomalies illustrated in Fig. LABEL:fig:anomalies 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. LABEL:fig:anomalies
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.

##  III 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.

Symbol | Name/Description
---|---
${V}$ | Set of vertices.
${G}_{\rm u}$, ${E}_{\rm u}$ | Undirected graph, undirected edge set.
$B_{\rm u}$, $L_{\rm u}$ | Incidence, and Laplacian matrices of
| an undirected graph [equation (30)].
$C_{\Sigma}$ | cycle-edge incidence matrix of the set
| of cycles $\Sigma$ [equation (35)].
${G}_{\rm d}$, ${E}_{\rm d}$ | Directed graph, set of directed edges.
${G}_{\rm b}$, ${E}_{\rm b}$ | Bidirected counterpart of the undirected
| graph ${G}_{\rm u}$ and its set of directed edges.
$s_{e}$, $t_{e}$ | Source and target of edge $e$.
$\bar{e}$ | Edge $e$ with opposite direction.
$A_{\rm d}$, $B_{\rm d}$, $L_{\rm d}$ | Adjacency, incidence, and Laplacian
| matrices of a digraph [equations (17), (24)].
$B_{\rm b}$ | Incidence matrix of a bidirected graph.
$B_{\rm o}$, $B_{\rm i}$ | Out- and in-incidence matrices of a
| digraph.
$\bm{\theta}=(\theta_{1},...,\theta_{n})^{\top}$ | Vector of phase angles.
$\bm{\omega}=(\omega_{1},...,\omega_{n})^{\top}$ | Vector of natural frequencies.
$a_{ij}$, $\phi_{ij}$ | Coupling strength and phase
| frustration between nodes $i$ and $j$.
$\gamma_{e}$ | Bound on the angular difference over
| the edge $e$.
$h_{e}$, $h_{\bar{e}}$ | Coupling functions over the edge $e$.
$f_{e}$, $g_{e}$ | Odd and even parts of the coupling
| over edge $e$ [equation (50), (51)].
$R({\bm{\gamma}})$ | Domain of bounded angular
| differences [equation (13)].
${\bf q}_{\Sigma}$ | Winding map for the cycles in $\Sigma$
| [equation (7)].
$\Omega({\bf u};\Sigma)$ | Winding cell with winding vector ${\bf u}$ in
| the graph ${G}_{\rm u}$ [equation (8)].
$S_{\epsilon}$ | Flow mismatch iteration [equation (49)].
Table 1: List of symbols.

Directed graphs. A directed graph (or digraph) ${G}_{\rm d}$ is the pair
$({V},{E}_{\rm d})$ composed of a set of vertices (or nodes)
${V}=\\{1,...,n\\}$ and a set of directed edges ${E}_{\rm
d}\subset{V}\times{V}$, which are ordered pairs of vertices. For an edge
$e=(i,j)\in{E}_{\rm 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 $\bar{e}=(t_{e},s_{e})$. The existence of edges is
encoded in the graph’s adjacency matrix

$\displaystyle(A_{\rm d})_{ij}$
$\displaystyle=\left\\{\begin{array}[]{ll}1\,,&\text{if }(i,j)\in{E}_{\rm
d}\,,\\\ 0\,,&\text{otherwise.}\end{array}\right.$ (17)

The out-degrees (resp. in-degrees) matrix is obtained as $D_{\rm o}={\rm
diag}(A_{\rm d}\bm{1})$ (resp. $D_{\rm i}={\rm diag}(A_{\rm
d}^{\top}\bm{1})$). We define the Laplacian matrix of ${G}_{\rm d}$ as $L_{\rm
d}=D_{\rm o}-A_{\rm d}$. For a digraph with $n$ vertices and $m$ directed
edges, we define the $n\times m$ out-incidence and in-incidence matrices

$\displaystyle(B_{\rm o})_{ie}$
$\displaystyle=\left\\{\begin{array}[]{ll}1\,,&\text{if }e=(i,j)\text{ for
some }j\,,\\\ 0\,,&\text{otherwise,}\end{array}\right.$ (20)
$\displaystyle(B_{\rm i})_{ie}$
$\displaystyle=\left\\{\begin{array}[]{ll}1\,,&\text{if }e=(j,i)\text{ for
some }j\,,\\\ 0\,,&\text{otherwise,}\end{array}\right.$ (23)

which form the standard incidence matrix

$\displaystyle B_{\rm d}$ $\displaystyle=B_{\rm o}-B_{\rm i}\,.$ (24)

We notice the following relations, the fourth being unknown as far as we can
tell.

######  Proposition 1.

The adjacency matrix $A_{\rm d}$, the out- and in-degree matrices $D_{\rm o}$
and $D_{\rm i}$, and the Laplacian matrix $L_{\rm d}$ of a directed graph can
be written in terms of its out- and in-incidence matrices $B_{\rm o}$ and
$B_{\rm i}$:

$\displaystyle D_{\rm o}$ $\displaystyle=B_{\rm o}B_{\rm o}^{\top}\,,$
$\displaystyle D_{\rm i}$ $\displaystyle=B_{\rm i}B_{\rm i}^{\top}\,,$
$\displaystyle A_{\rm d}$ $\displaystyle=B_{\rm o}B_{\rm i}^{\top}\,,$
$\displaystyle L_{\rm d}$ $\displaystyle=B_{\rm o}B^{\top}\,.$ (25)

###### Proof.

The proofs for the adjacency matrix $A_{\rm d}$ and for the degree matrices
$D_{\rm o}$ and $D_{\rm i}$ can be found in Ref. 63 (Lemmas 3.1 and 4.1). The
proof for the Laplacian matrix directly follows,

$\displaystyle L_{\rm d}$ $\displaystyle=D_{\rm o}-A_{\rm d}=B_{\rm o}B_{\rm
o}^{\top}-B_{\rm o}B_{\rm i}^{\top}$ $\displaystyle=B_{\rm o}(B_{\rm o}-B_{\rm
i})^{\top}=B_{\rm o}B^{\top}\,.$ (26)

∎

######  Remark.

The same proof is straightforwardly adapted to weighted directed graphs.

Undirected graphs. An (undirected) graph ${G}_{\rm u}$ is a pair
$({V},{E}_{\rm 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}_{\rm u}\subset\left\\{\\{i,j\\}\colon i,j\in{V}\right\\}$. A
cycle of ${G}_{\rm u}$ is an ordered sequence of vertices
$\sigma=(i_{0},i_{1},...,i_{\ell}=i_{0})$, such that
$\\{i_{j},i_{j+1}\\}\in{E}_{\rm u}$ and $i_{j}\neq i_{k}$ for any
$j,k\in\\{1,...,\ell\\}$.

Let us now choose an arbitrary orientation [$(i,j)$ or $(j,i)$] for each
undirected edge $\\{i,j\\}\in{E}_{\rm u}$. We can then define the incidence
matrix of ${G}_{\rm u}$,

$\displaystyle(B_{\rm u})_{ie}$
$\displaystyle=\left\\{\begin{array}[]{ll}1\,,&\text{if }e=(i,j)\text{ for
some }j\,,\\\ -1\,,&\text{if }e=(j,i)\text{ for some }j\,,\\\
0\,,&\text{otherwise.}\end{array}\right.$ (30)

The Laplacian matrix of ${G}_{\rm u}$ can be computed as $L_{\rm u}=B_{\rm
u}B_{\rm u}^{\top}$. Note that the incidence matrix is not unique and depends
on the choice of edge orientations, whereas the Laplacian does not. Given a
cycle $\sigma$, we define the cycle vector ${\bf
v}_{\sigma}\in\\{-1,0,+1\\}^{m}$, indexed by edges, as

$\displaystyle({\bf v}_{\sigma})_{e}$
$\displaystyle=\left\\{\begin{array}[]{ll}+1\,,&\text{if
}e=(i_{k-1},i_{k})\text{ for some }k\,,\\\ -1\,,&\text{if
}e=(i_{k},i_{k-1})\text{ for some }k\,,\\\
0\,,&\text{otherwise.}\end{array}\right.$ (34)

The cycle space of ${G}_{\rm u}$ is the span of the cycle vectors of all
cycles of ${G}_{\rm u}$, which is equivalently defined as the kernel of the
incidence matrix $B_{\rm u}$. A set of cycles
$\Sigma=\\{\sigma_{1},...,\sigma_{c}\\}$ is a cycle basis of ${G}_{\rm u}$ if
and only if the set of corresponding cycle vectors forms a basis of the cycle
space.

Finally, given a cycle basis $\Sigma$ of the graph ${G}_{\rm u}$, we define
the cycle-edge incidence matrix,

$\displaystyle C_{\Sigma}$ $\displaystyle=\left({\bf
v}_{\sigma_{1}},\cdots,{\bf
v}_{\sigma_{c}}\right)^{\top}\in\mathbb{R}^{c\times m}\,.$ (35)

Bidirected graphs. Dissipative couplings intrinsically require to distinguish
the two orientations of each edge. Given an undirected graph ${G}_{\rm
u}=({V},{E}_{\rm u})$, its bidirected counterpart is the directed graph
${G}_{\rm b}=({V},{E}_{\rm b})$ with the same vertex set ${V}$ and where each
undirected edge $\\{i,j\\}\in{E}_{\rm u}$ is doubled in the set of directed
edges $(i,j),(j,i)\in{E}_{\rm b}$. A bidirected graph is a directed graph
induced by an undirected graph.

If the undirected graph ${G}_{\rm u}$ has incidence matrix $B_{\rm
u}\in\mathbb{R}^{n\times m}$ [equation (30)], then, with an appropriate
indexing of the directed edges, the incidence matrix of ${G}_{\rm b}$ can be
written as $B_{\rm b}=(B_{\rm u},-B_{\rm u})\in\mathbb{R}^{n\times 2m}$.
Interestingly, we note that the Laplacian matrices of ${G}_{\rm u}$ and
${G}_{\rm b}$ are the same, namely (see Prop. 1),

$\displaystyle L_{\rm u}$ $\displaystyle=B_{\rm u}B_{\rm u}^{\top}=L_{\rm
b}=B_{\rm o}B_{\rm b}^{\top}\,,$ (36)

where the out-incidence matrix $B_{\rm o}=[B_{\rm b}]_{+}$ is the positive
part of $B_{\rm b}$. Notice that here,

$\displaystyle B_{\rm o}$ $\displaystyle=\left([B_{\rm u}]_{+},~{}[B_{\rm
u}]_{-}\right)\,.$ (37)

The Kuramoto-Sakaguchi model. We illustrate the results of this article with
the generalized Kuramoto-Sakaguchi model, 25; 8; 26

$\displaystyle\dot{\theta}_{i}$
$\displaystyle=\omega_{i}-\sum_{j=1}^{n}a_{ij}\left[\sin(\theta_{i}-\theta_{j}-\phi_{ij})+\sin(\phi_{ij})\right]\,,$
(38)

for $i\in\\{1,...,n\\}$, where $\theta_{i}\in\mathbb{S}^{1}$ and
$\omega_{i}\in\mathbb{R}$ are respectively the phase angle and the natural
frequency of the $i$-th oscillator, $\phi_{ij}\in(-\pi/2,\pi/2)$ is the phase
frustration between oscillators, and in this case, $a_{ij}\in\mathbb{R}_{\geq
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}\equiv 0$, $d_{i}\equiv 1$, $f_{ij}(x)=a_{ij}\cos(\phi_{ij})\sin(x)$,
and $g_{ij}(x)=a_{ij}\sin(\phi_{ij})[1-\cos(x)]$. It is a natural extension of
the original Kuramoto model, which is recovered for $\phi_{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 $\phi_{ij}=\phi_{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.

Figure 5: Coupling functions for the Kuramoto-Sakaguhchi model. (a) Comparison
between coupling functions for Kuramoto [dashed dark blue, $h_{e}(x)=\sin(x)$]
and Kuramoto-Sakaguchi [plain cyan, $h_{e}(x)=\sin(x-\phi)+\sin(\phi)$], with
$\phi=0.5$. The light green curve illustrates the coupling on the same edge,
but with opposite orientation [$h_{\bar{e}}(-x)=\sin(-x-\phi)+\sin(\phi)$].
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 (56).

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.

######  Example (2-node system, vectorial form).

Consider a system of two coupled Kuramoto-Sakaguchi oscillators with unit
coupling and identical frustration, whose dynamics is given by

$\displaystyle\begin{split}\dot{\theta}_{1}&=\omega_{1}-\left[\sin(\theta_{1}-\theta_{2}-\phi)+\sin(\phi)\right]\,,\\\
\dot{\theta}_{2}&=\omega_{2}-\left[\sin(\theta_{2}-\theta_{1}-\phi)+\sin(\phi)\right]\,,\end{split}$
(39)

with $\phi\in(-\pi/2,\pi/2)$.

In the Kuramoto model ($\phi=0$), it is standard to write the dynamics in
vectorial form as

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=\bm{\omega}-B_{\rm u}{\cal
A}\sin\left(B_{\rm u}^{\top}\bm{\theta}\right)\,,$ (40)

where $B_{\rm u}\in\mathbb{R}^{n\times m}$ is the incidence matrix of the
(undirected) coupling graph of the system, and ${\cal A}\in\mathbb{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

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=\bm{\omega}-B_{\rm u}{\cal
A}\left[\sin\left(B_{\rm
u}^{\top}\bm{\theta}-\phi\bm{1}_{m}\right)+\sin(\phi)\bm{1}_{m}\right]\,,$
(41)

does not yields the desired equations (39).

In order to write equation (38) in vectorial form, we need to distinguish the
two orientation of each edge and consider ${G}_{\rm b}$, the bidirected
counterpart of ${G}_{\rm u}$. We need to introduce the out-incidence matrix
$B_{\rm o}$ and the incidence matrix $B_{\rm b}$ of the bidirected coupling
graph, as well as ${\cal A}_{\rm b}\in\mathbb{R}^{2m}$, defined to be a block
diagonal matrix whose two diagonal blocks are equal to ${\cal A}$. The proof
of the following proposition follows from direct computation.

######  Proposition 2.

The Kuramoto-Sakaguchi model in equation (38) is written in vectorial form as

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=\bm{\omega}-B_{\rm o}{\cal
A}_{\rm b}\left[\sin\left(B_{\rm
b}^{\top}\bm{\theta}-\phi\bm{1}_{2m}\right)+\sin(\phi)\bm{1}_{2m}\right]\,.$
(42)

######  Example (6-node cycle, sync. frequency).

Let us consider six Kuramoto-Sakaguchi oscillators, coupled in a cycle, with
identical, vanishing natural frequency, i.e.,

$\displaystyle\dot{\theta}_{i}$
$\displaystyle=-\sin(\theta_{i}-\theta_{i-1}-\phi)-\sin(\theta_{i}-\theta_{i+1}-\phi)+2\sin(\phi)\,,$
(43)

for $i\in\\{1,...,6\\}$, where we used periodic indexing. One
straightforwardly verifies that $\bm{\theta}_{0}=(0,...,0)^{\top}$ is an
equilibrium of equation (43) (and then a synchronous state).

One can also verify that the splay state
$\bm{\theta}_{1}=(0,-\pi/3,-2\pi/3,\pi,2\pi/3,\pi/3)^{\top}$ is also a
synchronous state. Indeed, in this case, equation (43) gives

$\displaystyle\begin{split}\dot{\theta}_{i}&=-\sin(-\pi/3-\phi)-\sin(\pi/3-\phi)-2\sin(\phi)\\\
&=2\cos(\pi/3)\sin(\phi)-2\sin(\phi)=-\sin(\phi)\,,\end{split}$ (44)

independently of $i\in\\{1,...,6\\}$. The state $\bm{\theta}_{1}$ is then
synchronous, but it is an equilibrium only for the Kuramoto model ($\phi=0$).

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 (42). Note that a
similar vectorial formulation of the Kuramoto-Sakaguchi model has recently
been proposed in Ref. 64, which, while more general than equation (42), 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,

$\displaystyle\sum_{i}\dot{\theta}_{i}$
$\displaystyle=\sum_{i}\omega_{i}-\sum_{i,j}a_{ij}\sin(\phi)\left[1-2\cos(\theta_{i}-\theta_{j})\right]\,,$
(45)

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 $\phi=0$, 13; 24 are not
directly applicable when $\phi\neq 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. 65; 13
The power flow equations relate the the balance of active ($P_{i}$) and
reactive powers ($Q_{i}$) to the voltage amplitude ($V_{i}$) and phase
($\theta_{i}$) at each node $i\in\\{1,...,n\\}$,

$\displaystyle P_{i}$ $\displaystyle=\sum_{j=1}^{n}V_{i}V_{j}[{\cal
B}_{ij}\sin(\theta_{i}-\theta_{j})+{\cal
G}_{ij}\cos(\theta_{i}-\theta_{j})]\,,$ (46) $\displaystyle Q_{i}$
$\displaystyle=\sum_{j=1}^{n}V_{i}V_{j}[{\cal
G}_{ij}\sin(\theta_{i}-\theta_{j})-{\cal
B}_{ij}\cos(\theta_{i}-\theta_{j})]\,,$ (47)

with ${\cal G}_{ij}$ and ${\cal B}_{ij}$ being lines conductance and
susceptance respectively. Defining

$\displaystyle a_{ij}$ $\displaystyle=V_{i}V_{j}\sqrt{{\cal B}_{ij}^{2}+{\cal
G}_{ij}^{2}}\,,$ $\displaystyle\phi_{ij}$ $\displaystyle=\arctan(-{\cal
G}_{ij}/{\cal B}_{ij})\,,$ (48)

one verifies that solutions of equation (46) are steady states of equation
(38).

Equations (46) and (47) 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. 42 The 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_{\epsilon}$ [equation (49)], with $\epsilon=0.01$. All data are available
online. 41

The Dissipative Flow Network Problem on acyclic graphs. In the case where
${G}_{\rm 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}_{\rm u}$. Then there is at most one ${\bm{\Delta}}\in
R({\bm{\gamma}})$ 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. 41

######  Remark.

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. LABEL:fig:anomalies 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_{\epsilon}$ over the space
of angular differences $\mathbb{R}^{m}$, whose fixed points are exactly the
solutions of equation (14a). Namely, let

$\displaystyle\begin{split}S_{\epsilon}\colon\mathbb{R}^{m}&\to\mathbb{R}^{m}\\\
{\bm{\Delta}}&\mapsto{\bm{\Delta}}-\epsilon B_{\rm u}^{\top}L_{\rm
u}^{\dagger}\left(B_{\rm o}{\bf
h}({\bm{\Delta}})-\bm{\omega}\right)\,,\end{split}$ (49)

where $\epsilon>0$ is a small step size and $L_{\rm u}^{\dagger}$ is the
pseudoinverse of the graph Laplacian matrix. The flow mismatch iteration
$S_{\epsilon}$ updates the vector of angular differences according to the
mismatch between the input/output of commodities $\bm{\omega}$ and the
distribution of flows that corresponds to the current angular differences. It
has two major properties:

  1. (I)

the vector ${\bm{\Delta}}^{*}\in R({\bm{\gamma}})$ is a fixed point of
$S_{\epsilon}$ if and only if it is a solution of equation (14a);

  2. (II)

the map $S_{\epsilon}$ leaves each winding cell invariant, because
$C_{\Sigma}B_{\rm u}^{\top}=\bm{0}$. It means that fixing the winding vector
of the initial conditions, imposes the winding vector of the fixed point of
$S_{\epsilon}$, 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 ${\bf
u}\in\mathbb{Z}^{m-n+1}$, we are guaranteed that if the initial conditions
${\bm{\Delta}}_{0}$ satisfy equation (14b), then each following iteration
${\bm{\Delta}}_{k+1}=S_{\epsilon}({\bm{\Delta}}_{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({\bm{\gamma}})$.

######  Theorem 4.

If the flow mismatch iteration $S_{\epsilon}$ is contracting, then there is at
most one synchronous state of equation (2) in each winding cell.

In what follows, we provide sufficient conditions for contractivity of
$S_{\epsilon}$. We rely on the decomposition of the coupling functions as
$h_{e}(x)=f_{e}(x)+g_{e}(x)$, which implies,

$\displaystyle f_{e}(x)$ $\displaystyle=[h_{e}(x)-h_{\bar{e}}(-x)]/2\,,$ (50)
$\displaystyle g_{e}(x)$ $\displaystyle=[h_{e}(x)+h_{\bar{e}}(-x)]/2\,.$ (51)

One readily verifies the identities

$\displaystyle f_{e}^{\prime}(x)$ $\displaystyle=f_{\bar{e}}^{\prime}(-x)\,,$
$\displaystyle g_{e}^{\prime}(x)$ $\displaystyle=-g_{\bar{e}}^{\prime}(-x)\,.$
(52)

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 ${\bf x}\in R({\bm{\gamma}})$:

  1. (a)

the odd weighted Laplacian matrix, which is the Laplacian matrix of ${G}_{\rm
u}$ weighed by the derivatives of the odd parts

$\displaystyle L_{f}({\bf x})$ $\displaystyle=B_{\rm u}\cdot{\rm
diag}[f_{e}^{\prime}(x_{e})]\cdot B_{\rm u}^{\top}\,.$ (53)

We emphasize that the choice of orientation for each edge $e$ does not matter
in the definition of $L_{f}$. Also, the graph ${G}_{\rm u}$ being connected
and the above weights being positive, it is standard result of algebraic graph
theory that $\lambda_{2}$, the smallest nonzero eigenvalue of $L_{f}$ (a.k.a.,
the algebraic connectivity), is positive;

  2. (b)

the even weighted degree matrix, which is the diagonal matrix weighted by the
absolute derivatives of the even parts,

$\displaystyle[D_{g}({\bf x})]_{ii}$ $\displaystyle=\sum_{e\in
E_{i}}|g_{e}^{\prime}(x_{e})|\,,$ (54)

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

$\displaystyle h_{e}(x)$
$\displaystyle=h_{\bar{e}}(x)=a_{e}[\sin(x-\phi)+\sin(\phi)]\,,$ (55)

and trigonometric identities yield

$\displaystyle\begin{split}f_{e}(x)&=a_{e}\cos(\phi)\sin(x)\,,\\\
g_{e}(x)&=a_{e}\sin(\phi)[1-\cos(x)]\,,\end{split}$ (56)

which we illustrate in Fig. 5b. We clearly see here the relation between
$g_{e}$ and the dissipativity or frustration of the coupling. When $\phi=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_{\epsilon}$, and of the dissipation, that works against contractivity of
$S_{\epsilon}$. 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\in\\{1,...,n\\}$,

$\displaystyle\sup_{{\bf x}\in R({\bm{\gamma}})}(D_{g})_{ii}$
$\displaystyle<\inf_{{\bf x}\in R({\bm{\gamma}})}\lambda_{2}(L_{f})\,,$ (57)

then there exists a sufficiently small step size $\epsilon>0$ such that the
flow mismatch iteration $S_{\epsilon}$ [equation (49)] is contracting.

The left-hand-side of equation (57) quantifies the amount of dissipation that
is ”seen” at each node of the network, which vanishes for lossless couplings.
The right-hand-side accounts both for the strength of the coupling between
each pair of oscillators, through the weights, and for the connectedness of
the graph, $\lambda_{2}$ being the algebraic connectivity. 67 Under our
assumptions [${G}_{\rm u}$ is connected, couplings are strictly increasing on
$R({\bm{\gamma}})$], the right-hand-side of equation (57) is necessarily
positive. For couplings with sufficiently low dissipation, equation (57) is
then satisfied, which, combined with Theorem 4, yields the following
corollary.

######  Corollary 6.

If equation (57) 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
$\lambda_{2}(L_{f})$ in the Supplementary Information (Prop. S2) that are
state-independent and may ease the verification of equation (57). The bounds
are adapted from standard results of algebraic graph theory.

## Data and code availability

The data used and the codes created for this article have been deposited on
the DFNSolver 41 repository: http://dx.doi.org/10.5281/zenodo.5899408.

## References

  * Huygens (1993) C. Huygens,  _Oeuvres complètes de Christiaan Huygens_, edited by M. Nijhoff (Société hollandaise des sciences, 1993).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep. 469, 93–153 (2008).
  * Strogatz (2000) S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” Physica D 143, 1–20 (2000).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, “The Kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137 (2005).
  * Winfree (1967) A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” J. Theoret. Biol. 16, 15–42 (1967).
  * Kuramoto (1975) Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in  _Lecture Notes in Physics, vol. 39. International Symposium on Mathematical Problems in Theoretical Physics_, edited by Huzihiro Araki (Springer, Berlin, 1975) pp. 420–422.
  * Kuramoto (1984) Y. Kuramoto, “Cooperative dynamics of oscillator community a study based on lattice of rings,” Prog. Theor. Phys. Suppl. 79, 223–240 (1984).
  * Sakaguchi  _et al._ (1987) H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, “Local and global self-entrainments in oscillator lattices,” Prog. Theor. Phys. 77, 1005–1010 (1987).
  * Strogatz and Mirollo (1988) S. H. Strogatz and R. E. Mirollo, “Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random intrinsic frequencies,” Physica D 31, 143–168 (1988).
  * Daido (1988) H. Daido, “Lower critical dimension for populations of oscillators with randomly distributed frequencies: A renormalization-group analysis,” Phys. Rev. Lett. 61, 231–234 (1988).
  * Dörfler and Bullo (2014) F. Dörfler and F. Bullo, “Synchronization in complex networks of phase oscillators: A survey,” Automatica 50, 1539–1564 (2014).
  * Bergen and Hill (1981) A. R. Bergen and D. J. Hill, “A structure preserving model for power system stability analysis,” IEEE Trans. Power App. Syst. , 25–35 (1981).
  * Dörfler  _et al._ (2013) F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” Proc. Natl. Acad. Sci. USA 110, 2005–2010 (2013).
  * Paley  _et al._ (2007) D. A. Paley, N. E. Leonard, R. Sepulchre, D. Grunbaum, and J. K. Parrish, “Oscillator models and collective motion,” IEEE Control Syst. 27, 89–105 (2007).
  * Leonard  _et al._ (2012) N. E. Leonard, T. Shen, B. Nabet, L. Scardovi, I. D. Couzin, and S. A. Levin, “Decision versus compromise for animal groups in motion,” Proc. Natl. Acad. Sci. USA 109, 227–232 (2012).
  * Bullo (2022) F. Bullo,  _Lectures on Network Systems, v1.6_ (Kindle Direct Publishing, 2022).
  * Levis  _et al._ (2017) D. Levis, I. Pagonabarraga, and A. Díaz-Guilera, “Synchronization in Dynamical Networks of Locally Coupled Self-Propelled Oscillators,” Phys. Rev. X 7, 011028 (2017).
  * Simpson-Porco (2018a) J. W. Simpson-Porco, “A theory of solvability for lossless power flow equations – Part I: Fixed-point power flow,” IEEE Trans. Control Netw. Syst. 5, 1361–1372 (2018a).
  * Simpson-Porco (2018b) J. W. Simpson-Porco, “A theory of solvability for lossless power flow equations – Part II: Conditions for radial networks,” IEEE Trans. Control Netw. Syst. 5, 1373–1385 (2018b).
  * Jafarpour  _et al._ (2022) S. Jafarpour, E. Y. Huang, K. D. Smith, and F. Bullo, “Flow and elastic networks on the $n$-torus: Geometry, analysis, and computation,” SIAM Review 64, 59–104 (2022).
  * Aeyels and De Smet (2008) D. Aeyels and F. De Smet, “A mathematical model for the dynamics of clustering,” Physica D 237, 2517–2530 (2008).
  * Xia and Cao (2011) W. Xia and M. Cao, “Clustering in diffusively coupled networks,” Automatica 47, 2395–2405 (2011).
  * Dvijotham and Chertkov (2015) K. Dvijotham and M. Chertkov, “Convexity of structure preserving energy functions in power transmission: Novel results and applications,” in  _American Control Conference (ACC)_ (IEEE, Chicago, IL, USA, 2015) pp. 5035–5042.
  * Ling  _et al._ (2019) S. Ling, R. Xu, and A. S. Bandeira, “On the Landscape of Synchronization Networks: A Perspective from Nonconvex Optimization,” SIAM J. Optim. 29, 1879–1907 (2019).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, “A soluble active rotater model showing phase transitions via mutual entertainment,” Prog. Theor. Phys. 76, 576–581 (1986).
  * Sakaguchi  _et al._ (1988) H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, “Mutual entrainment in oscillator lattices with nonvariational type interaction,” Prog. Theor. Phys. 79, 1069–1079 (1988).
  * Balestra  _et al._ (2019) C. Balestra, F. Kaiser, D. Manik, and D. Witthaut, “Multistability in lossy power grids and oscillator networks,” Chaos 29, 123119 (2019).
  * Hellmann  _et al._ (2020) F. Hellmann, P. Schultz, P. Jaros, R. Levchenko, T. Kapitaniak, J. Kurths, and Y. Maistrenko, “Network-induced multistability through lossy coupling and exotic solitary states,” Nat. Commun. 11, 592 (2020).
  * Bronski  _et al._ (2018) J. C. Bronski, T. Carty, and L. DeVille, “Configurational stability for the Kuramoto-Sakaguchi model,” Chaos 28, 103109 (2018).
  * Berner  _et al._ (2021) R. Berner, S. Yanchuk, Y. Maistrenko, and E. Schöll, “Generalized splay states in phase oscillator networks,” Chaos 31, 073128 (2021).
  * Mihara and Medrano-T (2019) A. Mihara and R. O. Medrano-T, “Stability in the Kuramoto–Sakaguchi model for finite networks of identical oscillators,” Nonlin. Dyn. 98, 539–550 (2019).
  * Zhang and Motter (2020) Y. Zhang and A. E. Motter, “Symmetry-Independent Stability Analysis of Synchronization Patterns,” SIAM Review 62, 817–836 (2020).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, “Networks beyond pairwise interactions: Structure and dynamics,” Phys. Rep. 874, 1–92 (2020).
  * Zhang  _et al._ (2021) Y. Zhang, J. L. Ocampo-Espindola, I. Z. Kiss, and A. E. Motter, “Random heterogeneity outperforms design in network synchronization,” Proc. Natl. Acad. Sci. USA 118, e2024299118 (2021).
  * Matheny  _et al._ (2019) M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. Hudoba de Badyn, M. Pósfai, L. Duenas-Osorio, M. Mesbahi, J. P. Crutchfield, M. C. Cross, R. M. D’Souza, and M. L. Roukes, “Exotic states in a simple network of nanoelectromechanical oscillators,” Science 363, eaav7932 (2019).
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, “Chimera states for coupled oscillators,” Phys. Rev. Lett. 93, 174102 (2004).
  * Abrams  _et al._ (2008) D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, “Solvable model for chimera states of coupled oscillators,” Phys. Rev. Lett. 101, 084103 (2008).
  * Wolfrum and Omel’chenko (2011) M. Wolfrum and O. E. Omel’chenko, “Chimera states are chaotic transients,” Phys. Rev. E 84, 015201(R) (2011).
  * Martens  _et al._ (2013) E. A. Martens, S. Thutupalli, A. Fourriere, and O. Hallatschek, “Chimera states in mechanical oscillator networks,” Proc. Natl. Acad. Sci. USA 110, 10563–10567 (2013).
  * Zhang  _et al._ (2020) Y. Zhang, Z. G. Nicolaou, J. D. Hart, R. Roy, and A. E. Motter, “Critical Switching in Globally Attractive Chimeras,” Phys. Rev. X 10, 011044 (2020).
  * Delabays (2022) R. Delabays, “DFNSolver: Dissipative Flow Networks Solver (v1.1), Zenodo,”  (2022).
  * Grigg  _et al._ (1999) C. Grigg, P. Wong, P. Albrecht, R. Allan, M. Bhavaraju, R. Billinton, Q. Chen, C. Fong, S. Haddad, S. Kuruganty, W. Li, R. Mukerji, D. Patton, N. Rau, D. Reppen, A. Schneider, M. Shahidehpour, and C. Singh, “The IEEE Reliability Test System-1996. A report prepared by the Reliability Test System Task Force of the Application of Probability Methods Subcommittee,” IEEE Trans. Power Syst. 14, 1010–1020 (1999).
  * Delabays  _et al._ (2016) R. Delabays, T. Coletta, and Ph. Jacquod, “Multistability of phase-locking and topological winding numbers in locally coupled Kuramoto models on single-loop networks,” J. Math. Phys. 57, 032701 (2016).
  * Manik  _et al._ (2017) D. Manik, M. Timme, and D. Witthaut, “Cycle flows and multistability in oscillatory networks,” Chaos 27, 083123 (2017).
  * Feynman (1955) R. P. Feynman, “Chapter II Application of Quantum Mechanics to Liquid Helium,” in  _Progress in Low Temperature Physics_, Vol. 1 (Elsevier, 1955) pp. 17–53.
  * Abrikosov (1957) A. A. Abrikosov, “On the magnetic properties of superconductors of the second group,” Sov. Phys. JETP 5, 1174 (1957).
  * Korsak (1972) A. J. Korsak, “On the question of uniqueness of stable load-flow solutions,” IEEE Trans. Power App. Syst. 91, 1093–1100 (1972).
  * Ermentrout (1985) G. B. Ermentrout, “The behavior of rings of coupled oscillators,” J. Math. Biol. 23, 55–74 (1985).
  * Coletta  _et al._ (2016) T. Coletta, R. Delabays, I. Adagideli, and P. Jacquod, “Topologically protected loop flows in high voltage AC power grids,” New J. Phys. 18, 103042 (2016).
  * Khalil (2002) H. K. Khalil,  _Nonlinear systems_ , 3rd ed. (Prentice Hall, 2002).
  * Horn and Johnson (1994) R. A. Horn and C. R. Johnson,  _Matrix Analysis_ (Cambridge University Press, New York, 1994).
  * Tamura  _et al._ (1983) Y. Tamura, H. Mori, and S. Iwamoto, “Relationship between voltage instability and multiple load flow solutions in electric power systems,” IEEE Trans. Power App. Syst. 102, 1115–1125 (1983).
  * Klos and Wojcicka (1991) A. Klos and J. Wojcicka, “Physical aspects of the nonuniqueness of load flow solutions,” Int. J. Elect. Power Energy Syst. 13, 268 – 276 (1991).
  * Ma and Thorp (1993) W. Ma and J. S. Thorp, “An efficient algorithm to locate all the load flow solutions,” IEEE Trans. Power Syst. 8, 1077–1083 (1993).
  * Mehta  _et al._ (2016) D. Mehta, D. K. Molzahn, and K. Turitsyn, “Recent advances in computational methods for the power flow equations,” in  _American Control Conference (ACC)_ (IEEE, Boston, MA, USA, 2016) pp. 1753–1765.
  * Casazza (1998) J. Casazza, “Blackouts: Is the risk increasing?” Electr. World 212, 62 (1998).
  * Whitley (2008) S. G. Whitley,  _Lake Erie Loop Flow Mitigation_ (Technical Report; New York Independent System Operator, 2008).
  * Simpson-Porco  _et al._ (2016) J. W. Simpson-Porco, F. Dörfler, and F. Bullo, “Voltage collapse in complex power grids,” Nat. Commun. 7, 10790 (2016).
  * Schmietendorf  _et al._ (2019) K. Schmietendorf, O. Kamps, M. Wolff, P. G. Lind, P. Maass, and J. Peinke, “Bridging between load-flow and Kuramoto-like power grid models: A flexible approach to integrating electrical storage units,” Chaos 29, 103151 (2019).
  * Böttcher  _et al._ (2022) P. C. Böttcher, D. Witthaut, and L. Rydin Gorjão, “Dynamic stability of electric power grids: Tracking the interplay of the network structure, transmission losses, and voltage dynamics,” Chaos 32, 053117 (2022).
  * Dvijotham  _et al._ (2015) K. Dvijotham, M. Chertkov, and S. Low, “A differential analysis of the power flow equations,” in  _IEEE Conference on Decision and Control (CDC)_ (IEEE, Osaka, Japan, 2015) pp. 23–30.
  * Park  _et al._ (2021) S. Park, R. Y. Zhang, J. Lavaei, and R. Baldick, “Uniqueness of Power Flow Solutions Using Monotonicity and Network Topology,” IEEE Trans. Control Netw. Syst. 8, 319–330 (2021).
  * Balbuena  _et al._ (2003) C. Balbuena, D. Ferrero, X. Marcote, and I. Pelayo, “Algebraic properties of a digraph and its line digraph,” J. Interconnect. Netw. 4, 377–393 (2003).
  * Arnaudon  _et al._ (2021) A. Arnaudon, R. L. Peach, G. Petri, and P. Expert, “Connecting Hodge and Sakaguchi-Kuramoto: a mathematical framework for coupled oscillators on simplicial complexes,” arXiv preprint: 2111.11073 (2021).
  * Machowski  _et al._ (2008) J. Machowski, J. W. Bialek, and J. R. Bumby,  _Power System Dynamics_ , 2nd ed. (Wiley, Chichester, U.K, 2008).
  * Bergen and Vittal (2000) A. R. Bergen and V. Vittal,  _Power Systems Analysis_ (Prentice Hall, 2000).
  * Fiedler (1973) M. Fiedler, “Algebraic connectivity of graphs,” Czech. Math. J. 23, 298–305 (1973).

#### Acknowledgements.

RD was supported by the Swiss National Science Foundation, under grant number
P400P2_194359. This work was partly supported by AFOSR grant FA9550-22-1-0059.

#### Authors contributions.

All three authors designed the research; RD and SJ derived the mathematical
results; RD did the simulations; RD and FB analyzed the simulations and wrote
the paper.

#### Competing interests.

The authors declare that they have no competing interests.

Supplementary Information for

Multistability and Anomalies in Oscillator Models of Lossy Power Grids

Robin Delabays,* Saber Jafarpour, and Francesco Bullo

*Corresponding author. Email<EMAIL_ADDRESS>
### Content

Section S1

The proof of Theorem 3;

Section S2

The proof of Theorem 5;

Section S3

Additional bounds on the algebraic connectivity.

### Section S1. Section S1. Proof of Theorem 3

Without loss of generality, let us renumber the nodes such that node $1$ is a
leaf of the graph ${G}_{\rm u}$, and such that, for all $i\in\\{2,...,n-1\\}$,
node $i$ is a leaf of the subgraph $G_{i}\subset{G}_{\rm u}$ where nodes with
index up to $i-1$ are pruned. Furthermore, let us denote by $e_{i}$ the unique
edge that connects node $i\in\\{1,...,n-1\\}$ to the set of nodes
$\\{i+1,...,n\\}$, and orient it such that $i$ is the source of $e$, i.e.,
$s_{e_{i}}=i$. The reversed edge is denoted $e_{i+m}=\bar{e}_{i}$. All
directed edges are then indexed. This can be done iteratively and we provide
an implementation of this renumbering in Ref. 1.

For convenience, we recall and slightly rephrase equation (15a),

$\displaystyle\varphi$
$\displaystyle=\omega_{j}-\sum_{\begin{subarray}{c}e\colon\\\
s_{e}=j\end{subarray}}h_{e}(\Delta_{e})\,,$ $\displaystyle j$
$\displaystyle\in\\{1,...,n\\},~{}\varphi\in\mathbb{R}\,,$ (S58a)
$\displaystyle|\Delta_{e}|$ $\displaystyle\leq\gamma_{e}\,,$ $\displaystyle e$
$\displaystyle\in{E}_{\rm u}\,.$ (S58b)

A solution ${\bm{\Delta}}$ of equation (15a) is also a solution of equations
(S58). Showing that there is at most one solution to equations (S58) then
implies that there is at most one solution to equation (15a).

Assume that ${\bm{\Delta}}$ and ${\bm{\Delta}}^{\prime}$ are two solutions of
equation (15a), with respective synchronous frequencies $\varphi$ and
$\varphi^{\prime}$. Without loss of generality, assume that
$\varphi\leq\varphi^{\prime}$. We compare now the components of the two
solutions recursively.

Initial step, ${i=1}$. According to our choice of indexing, equation (S58a)
gives

$\displaystyle h_{e_{1}}(\Delta_{e_{1}})$
$\displaystyle=\omega_{1}-\varphi\geq\omega_{1}-\varphi^{\prime}=h_{e_{1}}(\Delta_{e_{1}}^{\prime})\,,$
(S59)

and by monotonicity of the coupling functions,

$\displaystyle\Delta_{e_{1}}$ $\displaystyle\geq\Delta_{e_{1}}^{\prime}\,.$
(S60)

Recursion step, ${2\leq i\leq n-1}$. By the previous steps, we have
$\Delta_{e_{j}}\geq\Delta_{e_{j}}^{\prime}$ for $j\in\\{1,...,i-1\\}$.
According to our choice of indexing, there is a single out-going edge from
node $i$ whose angle differences in the two solutions have not been compared
in the previous steps, namely $e_{i}$ (see Fig. S6). Again by equation (S58a)
and monotonicity of the coupling functions, we get

$\displaystyle h_{e_{i}}(\Delta_{e_{i}})$
$\displaystyle=\omega_{i}-\sum_{\begin{subarray}{c}e:s_{e}<i,\\\
t_{e}=i\end{subarray}}h_{\bar{e}}(-\Delta_{e})-\varphi\geq\omega_{i}-\sum_{\begin{subarray}{c}e:s_{e}<i,\\\
t_{e}=i\end{subarray}}h_{\bar{e}}(-\Delta_{e}^{\prime})-\varphi^{\prime}=h_{e_{i}}(\Delta_{e_{i}}^{\prime})\,,$
(S61)

and

$\displaystyle\Delta_{e_{i}}$ $\displaystyle\geq\Delta_{e_{i}}^{\prime}\,.$
(S62)

Final step, ${i=n}$. All the previous steps together with equation (S58a) and
monotonicity of the coupling functions give

$\displaystyle\varphi$
$\displaystyle=\omega_{n}-\sum_{\begin{subarray}{c}e:\\\
t_{e}=n\end{subarray}}h_{\bar{e}}(-\Delta_{e})\geq\omega_{n}-\sum_{\begin{subarray}{c}e:\\\
t_{e}=n\end{subarray}}h_{\bar{e}}(-\Delta_{e}^{\prime})=\varphi^{\prime}\geq\varphi\,,$
(S63)

which implies $\varphi=\varphi^{\prime}$. The inequalities in equations (S59),
(S60), (S61), (S62), and (S63) are then equalities and the two solutions are
identical. Note that we crucially used that the coupling functions are
strictly increasing. ∎

Figure S6: Illustration of the node and edge indexing. In this example, our
construction implies $j_{1},j_{2}<i\leq m$.

### Section S2. Section S2. Proof of Theorem 5

For technical purposes, we need to define the extended flow function ${\bf
h}_{\gamma}\colon\mathbb{R}^{m}\to\mathbb{R}^{2m}$, whose components are
straight lines extensions of the components of ${\bf h}$. For each edge $e$,

$\displaystyle\left[{\bf h}_{\gamma}({\bf y})\right]_{e}$
$\displaystyle=\left\\{\begin{array}[]{ll}h_{e}^{\prime}(-\gamma_{e})(y_{e}+\gamma_{e})+h_{e}(-\gamma_{e})\,,&y_{e}<-\gamma_{e}\,,\\\
h_{e}(y_{e})\,,&|y_{e}|\leq\gamma_{e}\,,\\\
h_{e}^{\prime}(\gamma_{e})(y_{e}-\gamma_{e})+h_{e}(\gamma_{e})\,,&y_{e}>\gamma_{e}\,,\end{array}\right.$
(S67) $\displaystyle\left[{\bf h}_{\gamma}({\bf y})\right]_{e+m}$
$\displaystyle=\left\\{\begin{array}[]{ll}-h_{\bar{e}}^{\prime}(\gamma_{e})(y_{e}+\gamma_{e})+h_{\bar{e}}(\gamma_{e})\,,&y_{e}<-\gamma_{e}\,,\\\
h_{\bar{e}}(-y_{e})\,,&|y_{e}|\leq\gamma_{e}\,,\\\
-h_{\bar{e}}^{\prime}(-\gamma_{e})(y_{e}-\gamma_{e})+h_{\bar{e}}(-\gamma_{e})\,,&y_{e}>\gamma_{e}\,,\end{array}\right.$
(S71)

which is well-defined and continuously differentiable.

Let $\bm{\xi},\bm{\eta}\in\mathbb{R}^{m}$, such that
$C_{\Sigma}\bm{\xi}=C_{\Sigma}\bm{\eta}={\bf u}\in\mathbb{Z}^{c}$, and define
${\bf y}=\bm{\xi}-\bm{\eta}\in{\rm Ker}(C_{\Sigma})$. We construct the two
diagonal $m\times m$ matrices $\Lambda_{1}(\bm{\xi},\bm{\eta})$,
$\Lambda_{2}(\bm{\xi},\bm{\eta})$ as

$\displaystyle(\Lambda_{1})_{e}$
$\displaystyle=(\xi_{e}-\eta_{e})^{-1}\int_{\eta_{e}}^{\xi_{e}}({\bf
h}_{\gamma})^{\prime}_{e}(t)dt\,,$ (S72) $\displaystyle(\Lambda_{2})_{e}$
$\displaystyle=(\xi_{e}-\eta_{e})^{-1}\int_{\eta_{e}}^{\xi_{e}}({\bf
h}_{\gamma})^{\prime}_{\bar{e}}(-t)dt\,.$ (S73)

We verify that

$\displaystyle{\bf h}_{\gamma}(\bm{\xi})-{\bf h}_{\gamma}(\bm{\eta})$
$\displaystyle=\begin{pmatrix}\Lambda_{1}\\\
-\Lambda_{2}\end{pmatrix}(\bm{\xi}-\bm{\eta})=\Lambda{\bf y}\,,$ (S74)

and the two matrices $\Lambda_{1},\Lambda_{2}$ are nonnegative, because the
coupling functions are assumed strictly increasing.

Then, by definition of $S_{\epsilon}$,

$\displaystyle\left\|S_{\epsilon}(\xi)-S_{\epsilon}(\eta)\right\|_{2}^{2}={\bf
y}^{\top}\left(I_{m}-\epsilon B_{\rm u}^{\top}L_{\rm u}^{\dagger}B_{\rm
o}\Lambda\right)^{\top}\left(I_{m}-\epsilon B_{\rm u}^{\top}L_{\rm
u}^{\dagger}B_{\rm o}\Lambda\right){\bf y}=\|{\bf y}\|_{2}^{2}-\epsilon{\bf
y}^{\top}M{\bf y}+O(\epsilon^{2})\,,$ (S75)

where

$\displaystyle M$ $\displaystyle=B_{\rm u}^{\top}L_{\rm u}^{\dagger}B_{\rm
o}\Lambda+\Lambda^{\top}B_{\rm o}^{\top}(L_{\rm u}^{\dagger})^{\top}B_{\rm
u}\,.$ (S76)

The remainder of the proof will show that, under the assumptions of the
theorem, ${\bf y}^{\top}M{\bf y}>0$. Then for $\epsilon>0$ small enough, the
second order term in equation (S75) is dominated by the first order term,
which is negative. Therefore, equation (S75) is strictly smaller than $\|{\bf
y}\|_{2}^{2}$, and $S_{\epsilon}$ is contracting.

We remark that $B_{\rm o}\Lambda$ has exactly the same sparsity and sign
pattern as $B_{\rm u}$, with the difference being that the nonzero terms are
not all the same, not even for the terms corresponding to the ends of a given
edge. Namely,

$\displaystyle\left(B_{\rm o}\Lambda\right)_{ie}$
$\displaystyle=\left\\{\begin{array}[]{ll}(\Lambda_{1})_{ee}\,,&\text{if
}e=(i,j)\text{ for some }j\,,\\\ -(\Lambda_{2})_{ee}\,,&\text{if
}e=(j,i)\text{ for some }j\,,\\\ 0\,,&\text{otherwise.}\end{array}\right.$
(S80)

We illustrate this property for a simple network in the example below.

We now define

$\displaystyle\Lambda_{\rm p}$
$\displaystyle=\frac{\Lambda_{1}+\Lambda_{2}}{2}\,,$
$\displaystyle\Lambda_{\rm m}$
$\displaystyle=\frac{\Lambda_{1}-\Lambda_{2}}{2}\,,$ (S81)

and recall that the out-incidence matrix is given by $B_{\rm o}=([B_{\rm
u}]_{+},[B_{\rm u}]_{-})$ in equation (25), where the square brackets denote
the positive and negative parts (see the Methods Section for details). Then we
can rewrite

$\displaystyle B_{\rm o}\Lambda$ $\displaystyle=[B_{\rm
u}]_{+}\Lambda_{1}-[B_{\rm u}]_{-}\Lambda_{2}=[B_{\rm u}]_{+}(\Lambda_{\rm
p}+\Lambda_{\rm m})-[B_{\rm u}]_{-}(\Lambda_{\rm p}-\Lambda_{\rm m})=B_{\rm
u}\Lambda_{\rm p}+|B_{\rm u}|\Lambda_{\rm m}\,,$ (S82)

with the absolute value taken elementwise.

######  Example.

We illustrate equations (S80) and (S82) for the case of a triangular network,
shown in Fig. S7. The matrices involved are

$\displaystyle B_{\rm o}$ $\displaystyle=\begin{pmatrix}1&0&0&0&0&1\\\
0&1&0&1&0&0\\\ 0&0&1&0&1&0\end{pmatrix}\,,$ $\displaystyle\Lambda$
$\displaystyle=\begin{pmatrix}\lambda_{1,1}&0&0\\\ 0&\lambda_{1,2}&0\\\
0&0&\lambda_{1,3}\\\ -\lambda_{2,1}&0&0\\\ 0&-\lambda_{2,2}&0\\\
0&0&-\lambda_{2,3}\end{pmatrix}\,,$ $\displaystyle B_{\rm o}\Lambda$
$\displaystyle=\begin{pmatrix}\lambda_{1,1}&0&-\lambda_{2,3}\\\
-\lambda_{2,1}&\lambda_{1,2}&0\\\
0&-\lambda_{2,2}&\lambda_{1,3}\end{pmatrix}\,,$ (S83)

illustrating equation (S80). Now defining $\lambda_{{\rm
p},i}=(\lambda_{1,i}+\lambda_{2,i})/2$ and $\lambda_{{\rm
m},i}=(\lambda_{1,i}-\lambda_{2,i})/2$, we can decompose $B_{\rm o}\Lambda$ as

$\displaystyle B_{\rm o}\Lambda$ $\displaystyle=\begin{pmatrix}\lambda_{{\rm
p},1}&0&-\lambda_{{\rm p},3}\\\ -\lambda_{{\rm p},1}&\lambda_{{\rm p},2}&0\\\
0&-\lambda_{{\rm p},2}&\lambda_{{\rm
p},3}\end{pmatrix}+\begin{pmatrix}\lambda_{{\rm m},1}&0&\lambda_{{\rm m},3}\\\
\lambda_{{\rm m},1}&\lambda_{{\rm m},2}&0\\\ 0&\lambda_{{\rm
m},2}&\lambda_{{\rm m},3}\end{pmatrix}\,,$

which is precisely equation (S82).

Figure S7: Triangle network. Illustration of the example’s triangle network.

Based on equations (S76) and (S82), on the fact that ${\bf y}=B_{\rm
u}^{\top}L^{\dagger}B_{\rm u}{\bf y}$ (see Prop. S7), and using ${\bf
z}=L_{\rm u}^{\dagger}B_{\rm u}{\bf y}$, we compute

$\displaystyle{\bf y}^{\top}M{\bf y}$ $\displaystyle={\bf
z}^{\top}\left(B_{\rm o}\Lambda B_{\rm u}^{\top}+B_{\rm u}\Lambda B_{\rm
o}^{\top}\right){\bf z}={\bf z}^{\top}\left(2B_{\rm u}\Lambda_{\rm p}B_{\rm
u}^{\top}+|B_{\rm u}|\Lambda_{\rm m}B_{\rm u}^{\top}+B_{\rm u}\Lambda_{\rm
m}|B_{\rm u}|^{\top}\right){\bf z}=2{\bf z}^{\top}\left(L_{\rm p}+D_{\rm
m}\right){\bf z}\,,$ (S84)

where

$\displaystyle L_{\rm p}$ $\displaystyle=B_{\rm u}\Lambda_{\rm p}B_{\rm
u}^{\top}$ (S85) $\displaystyle D_{\rm m}$ $\displaystyle=\left(|B_{\rm
u}|\Lambda_{\rm m}B_{\rm u}^{\top}+B_{\rm u}\Lambda_{\rm m}|B_{\rm
u}|^{\top}\right)/2\,.$ (S86)

First, notice that $L_{\rm p}$ is a weighted Laplacian matrix with the graph
structure of ${G}_{\rm u}$. The weight of edge $e$ is given by

$\displaystyle(\Lambda_{\rm p})_{ee}$
$\displaystyle=(\xi_{e}-\eta_{e})^{-1}\int_{\eta_{e}}^{\xi_{e}}f_{e}^{\prime}(t)dt\,,$
(S87)

which, by the Mean Value Theorem, is bounded by

$\displaystyle\inf_{x}f_{e}^{\prime}(x)$ $\displaystyle\leq(\Lambda_{\rm
p})_{e}\leq\sup_{x}f_{e}^{\prime}(x)\,,$ (S88)

with the infimum and supremum taken over the admissible values of $x$.
Therefore, by definition of the odd weighted Laplacian matrix [equation (41)],
we can bound the second eigenvalue of $L_{\rm p}$ by Weyl’s inequalities, 2

$\displaystyle\lambda_{2}(L_{\rm p})$ $\displaystyle\geq\inf_{\bf
x}\lambda_{2}\left(L_{f}({\bf x})\right)\,.$ (S89)

Finally, using that ${\bf z}=L^{\dagger}B_{\rm u}{\bf y}$ is orthogonal to the
null space of $L_{\rm p}$,

$\displaystyle{\bf z}^{\top}L_{\rm p}{\bf z}$
$\displaystyle\geq\lambda_{2}(L_{\rm p})\|{\bf z}\|_{2}^{2}\geq\inf_{\bf
x}\lambda_{2}(L_{f}({\bf x}))\|{\bf z}\|_{2}^{2}\,.$ (S90)

Second, let us break down the matrix $D_{\rm m}$. Using the positive and
negative parts of the incidence matrix, $B_{\rm u}=[B_{\rm u}]_{+}-[B_{\rm
u}]_{-}$ and $|B_{\rm u}|=[B_{\rm u}]_{+}+[B_{\rm u}]_{-}$, direct computation
shows

$\displaystyle D_{\rm m}$ $\displaystyle=[B_{\rm u}]_{+}\Lambda_{\rm m}[B_{\rm
u}]_{+}^{\top}-[B_{\rm u}]_{-}\Lambda_{\rm m}[B_{\rm u}]_{-}^{\top}\,.$ (S91)

According to Prop. 1, we know that $D_{\rm m}$ is diagonal, and its elements
can be computed as

$\displaystyle\left(D_{\rm m}\right)_{ii}$
$\displaystyle=\sum_{e}\left([B_{\rm u}]_{+}\right)_{ie}^{2}(\Lambda_{\rm
m})_{ee}-\sum_{e}\left([B_{\rm u}]_{-}\right)_{ie}^{2}(\Lambda_{\rm
m})_{ee}=\sum_{e}\left([B_{\rm u}]_{+}-[B_{\rm
u}]_{-}\right)_{ie}(\Lambda_{\rm m})_{ee}=\sum_{e\in E_{i}}\pm(\Lambda_{\rm
m})_{ee}\,,$ (S92)

where $E_{i}$ is the set of edges incident to node $i$ in ${G}_{\rm u}$. Note
that the elements of $\Lambda_{\rm m}$ are given by

$\displaystyle(\Lambda_{\rm m})_{e}$
$\displaystyle=(\xi_{e}-\eta_{e})^{-1}\int_{\eta_{e}}^{\xi_{e}}g_{e}^{\prime}(t)dt\,,$
(S93)

and, by the Mean Value Theorem, are bounded by

$\displaystyle\left|(\Lambda_{\rm m})_{e}\right|$
$\displaystyle\leq\sup_{x}|g_{e}^{\prime}(x)|\,.$ (S94)

Therefore, by definition of the even weighted degree matrix [equation (42)],
we can bound the elements of $D_{\rm m}$,

$\displaystyle\min_{i}(D_{\rm m})_{ii}$ $\displaystyle\geq-\sup_{{\bf
x},i}(D_{g}({\bf x}))_{ii}\,,$ (S95)

which gives the bound

$\displaystyle{\bf z}^{\top}D_{\rm m}{\bf z}$
$\displaystyle\geq\min_{i}\left(D_{\rm m}\right)_{ii}\|{\bf
z}\|_{2}^{2}\geq-\sup_{{\bf x},i}\left(D_{g}({\bf x})\right)_{ii}\|{\bf
z}\|_{2}^{2}\,.$ (S96)

To conclude, we introduce equations (S90) and (S96) into equation (S84),
yielding

$\displaystyle{\bf y}^{\top}M{\bf y}$ $\displaystyle\geq\left[\inf_{\bf
x}\lambda_{2}\left(L_{f}({\bf x})\right)-\sup_{{\bf x},i}\left(D_{g}({\bf
x})\right)\right]\|{\bf z}\|_{2}^{2}\,,$ (S97)

which is strictly positive under the assumptions of the theorem. Going back to
equation (S75), we have shown that the first order term is strictly negative,
and therefore, for $\epsilon>0$ sufficiently small, the flow mismatch
iteration $S_{\epsilon}$ is contracting, which concludes the proof. ∎

######  Proposition S7.

Let ${G}$ be a graph and define its incidence matrix $B$, Laplacian matrix
$L$, and the cycle-edge incidence matrix $C_{\Sigma}$ (see the Methods Section
for definitions). Then

$\displaystyle{\rm Ker}(C_{\Sigma})$ $\displaystyle={\rm
Ker}(I_{m}-B^{\top}L^{\dagger}B)\,.$ (S98)

###### Proof.

Let ${\bf x}\in{\rm Ker}(I_{m}-B^{\top}L^{\dagger}B)$, then ${\bf
x}=B^{\top}L^{\dagger}B{\bf x}$. We compute

$\displaystyle C_{\Sigma}{\bf x}$
$\displaystyle=C_{\Sigma}B^{\top}L^{\dagger}B{\bf x}=0\,,$ (S99)

because $C_{\Sigma}B^{\top}=0$. Thus ${\rm
Ker}(I_{m}-B^{\top}L^{\dagger}B)\subset{\rm Ker}(C_{\Sigma})$.

The rows of $C_{\Sigma}$ are linearly independent by definition. Then its
kernel has dimension $m-(m-n+1)=n-1$. The matrix $I_{m}-B^{\top}L^{\dagger}B$
is the orthogonal projection onto the kernel of $B$, therefore its rank is the
nullity of $B$. By the Rank-Nullity Theorem,

$\displaystyle{\rm null}(I_{m}-B^{\top}L^{\dagger}B)$ $\displaystyle=m-{\rm
null}(B)=n-1\,.$ (S100)

The two kernel have the same dimension.

The set ${\rm Ker}(I_{m}-B^{\top}L^{\dagger}B)$ is then a subspace of ${\rm
Ker}(C_{\Sigma})$ and has the same dimension, they are then identical. ∎

### Section S3. Section S3. Bounds on the algebraic connectivity

The following bounds are adapted from standard results of algebraic graph
theory. We summarize the bounds of interest and the ad hoc quantities in Table
S2.

######  Proposition S8.

With the definition of $L_{f}$ given in equation (41), we have the following
bounds on its Fiedler eigenvalue:

  1. (i)

$\lambda_{2}\geq 2c_{\rm e}[1-\cos(\pi/n)]$ (Ref. 3, paragraph 4.3);

  2. (ii)

$\lambda_{2}\geq(nd_{\rm w})^{-1}$ (adapted from Ref. 4, Lemma 1.9);

  3. (iii)

$\lambda_{2}\geq 4\min_{e}\sigma_{\min,e}/nd_{{G}_{\rm u}}$ (adapted from Ref.
5, Theorem 4.2).

All relevant quantities are defined in Tab. S2.

###### Proof.

(i). Defining the weighted edge connectivity, 3

$\displaystyle c({\bf x})$ $\displaystyle=\min_{\cal
S}\sum_{e}f_{e}^{\prime}(x_{e})\,,$ (S101)

where the minimum is taken over subsets of edges ${\cal S}\subset{E}_{\rm u}$
that split the graph ${G}_{\rm u}$, we can adapt the proof of Ref. 3, in
paragraph 4.3, yielding

$\displaystyle\lambda_{2}(L_{f})$ $\displaystyle\geq 2c({\bf
x})\left[1-\cos(\pi/n)\right]\geq 2c_{\rm e}\left[1-\cos(\pi/n)\right]\,,$
(S102)

independently of ${\bf x}$, where $c_{\rm e}$ is defined in Tab. S2.

(ii). We adapt here the proof of Ref. 4, Lemma 1.9. Let ${\bf v}$ be the
eigenvector of $L_{f}$ associated with $\lambda_{2}$ and assume that
$|v_{i}|=\max_{k}|v_{k}|$ (recall that all these quantities depend on ${\bf
x}$). Because $L_{f}$ is a symmetric Laplacian matrix, ${\bf 1}^{\top}{\bf
v}=0$ and there is an index $j$ such that $v_{i}v_{j}<0$. We denote with
$P_{ij}$ the shortest (weighted) path from $i$ to $j$. Now,

$\displaystyle\lambda_{2}$ $\displaystyle=\frac{{\bf v}^{\top}L_{f}{\bf
v}}{{\bf v}^{\top}{\bf
v}}=\frac{\sum_{e=(k,\ell)}f_{e}^{\prime}(x_{e})(v_{k}-v_{\ell})^{2}}{\sum_{k}v_{k}^{2}}\geq\sum_{e=(k,\ell)\in
P_{ij}}\frac{f_{e}^{\prime}(x_{e})(v_{k}-v_{\ell})^{2}}{nv_{i}^{2}}\,.$ (S103)

Defining the odd weighted diameter

$\displaystyle D_{\rm w}({\bf x})$
$\displaystyle=\max_{i,j}\min_{P_{ij}}\sum_{e\in
P_{ij}}[f_{e}^{\prime}(x_{e})]^{-1}\,,$ (S104)

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 6 (direct consequence of the Cauchy-Schwarz inequality) to the
numerator of equation (S103) and get

$\displaystyle\lambda_{2}$ $\displaystyle\geq\frac{[D_{\rm w}({\bf
x})]^{-1}(v_{i}-v_{j})^{2}}{nv_{i}^{2}}\geq\frac{1}{nD_{\rm w}({\bf
x})}\geq\frac{1}{nd_{\rm w}}\,,$ (S105)

where $d_{\rm w}$ is defined in Tab. S2.

(iii). Alternatively, both sides of the identity

$\displaystyle{\bf v}^{\top}L_{f}{\bf v}$ $\displaystyle=\lambda_{2}{\bf
v}^{\top}{\bf v}\,,$ (S106)

can be bounded as follows,

$\displaystyle{\bf v}^{\top}L_{f}{\bf v}$
$\displaystyle=\sum_{e=(i,j)\in{E}_{\rm
u}}f_{e}^{\prime}(x_{e})(v_{i}-v_{j})^{2}\geq\min_{e}\sigma_{\min,e}\,,$
(S107)

and

$\displaystyle\begin{split}2n{\bf v}^{\top}{\bf
v}&=\sum_{i}\sum_{j}(v_{i}-v_{j})^{2}\leq\sum_{i}\sum_{j}|P_{ij}|\sum_{(k,\ell)\in
P_{ij}}(v_{k}-v_{\ell})^{2}\leq\sum_{(k,\ell)\in{E}_{\rm
u}}(v_{k}-v_{\ell})^{2}d_{{G}_{\rm u}}\sum_{i}\sum_{j}\chi_{ij}(k,\ell)\\\
&\leq\sum_{(k,\ell)\in{E}_{\rm u}}(v_{k}-v_{\ell})^{2}d_{{G}_{\rm
u}}n^{2}/2\,,\end{split}$ (S108)

where we used that ${\bf v}^{\top}\bm{1}=0$ and the Cauchy-Schwartz inequality
at the first line, $P_{ij}$ is a chosen (unweighted) shortest path between $i$
and $j$, $\chi_{ij}(e)$ is its indicator function, and we used Lemma 4.1 in
Ref. 5 at the last inequality. Plugging these bounds together yields

$\displaystyle 2n\min_{e}\sigma_{\min,e}$
$\displaystyle\leq\lambda_{2}d_{{G}_{\rm u}}n^{2}/2\,,$ (S109)

which concludes the proof. ∎

Bounds | Reference
---|---
$\displaystyle\max_{i}(D_{g})_{ii}\leq\max_{i}\sum_{e\in E_{i}}\alpha_{\max,e}$ |  Direct computation ($E_{i}$ is the set of edges incident to node $i$).
$\displaystyle\lambda_{2}\geq 2c_{\rm e}[1-\cos(\pi/n)]$ |  Prop. S8, adapted from Ref. 3, paragraph 4.3).
$\displaystyle\lambda_{2}\geq(nd_{\rm w})^{-1}$ |  Prop. S8, adapted from Lemma 1.9 in Ref. 4.
$\displaystyle\lambda_{2}\geq 4\min_{e}\sigma_{\min,e}/nd_{{G}_{\rm u}}$ |  Prop. S8, adapted from Theorem 4.2 in Ref. 5.
Definitions | Name
$\displaystyle\alpha_{\max,e}=\sup_{x}|g_{e}^{\prime}(x)|\phantom{\sum}$ | Maximal even slope.
$\displaystyle\sigma_{\min,e}=\inf_{x}f_{e}^{\prime}(x)\phantom{\sum}$ | Minimal odd slope.
$\displaystyle c_{\rm e}=\min_{E\subset{\cal S}}\sum_{e\in E}\sigma_{\min,e}$ |  Minimized odd weighted edge connectivity (${\cal S}$ is the set of splitting edge sets).
$\displaystyle d_{\rm w}=\max_{i,j}\min_{P_{ij}}\sum_{e\in P_{ij}}\sigma_{\min,e}^{-1}$ |  Maximized odd weighted diameter ($P_{ij}$ denotes a path between $i$ and $j$).
$\displaystyle d_{{G}_{\rm u}}=\max_{i,j}\min_{P_{ij}}|P_{ij}|$ |  Graph diameter ($P_{ij}$ denotes a path between $i$ and $j$).
Table S2: List of bounds on the components of equation (18) and definition of
the ad hoc quantities. Note that all quantities are independent of the state
of the system and can be determined beforehand.

## References

  * Delabays (2022) R. Delabays, “DFNSolver: Dissipative Flow Networks Solver (v1.1), Zenodo,”  (2022).
  * Horn and Johnson (1994) R. A. Horn and C. R. Johnson,  _Matrix Analysis_ (Cambridge University Press, New York, 1994).
  * Fiedler (1973) M. Fiedler, “Algebraic connectivity of graphs,” Czech. Math. J. 23, 298–305 (1973).
  * Chung (1997) F. R. K. Chung,  _Spectral Graph Theory_ , CBMS Regional Conference Series in Mathematics ; Number 92 (Published for the Conference Board of the Mathematical Sciences by the American Mathematical Society, Providence, Rhode Island, 1997).
  * Mohar (1991) B. Mohar, “Eigenvalues, diameter, and mean distance in graphs,” Graphs and Combin. 7, 53–64 (1991).
  * Sedrakyan and Sedrakyan (2018) H. Sedrakyan and N. Sedrakyan,  _Algebraic Inequalities_ , Problem Books in Mathematics (Springer International Publishing, Cham, 2018).