Multistability in lossy power grids and oscillator networks
Chiara Balestra, Franz Kaiser, Debsankha Manik, Dirk Witthaut
MMultistability in lossy power grids and oscillator networks
Chiara Balestra,
1, 2, a) Franz Kaiser,
1, 3, b) Debsankha Manik, c) and Dirk Witthaut
1, 3, d) Institute for Theoretical Physics, University of Cologne, Köln, 50937, Germany Department of Mathematics "Giuseppe Peano", Università degli Studi di Torino, 10123 Torino,Italy Forschungszentrum Jülich, Institute for Energy and Climate Research (IEK-STE), 52428 Jülich,Germany Max Planck Institute for Dynamics and Self-Organization, Am Faßberg 17, 37077 Göttingen,Germany (Dated: 9 January 2020)
Networks of phase oscillators are studied in various contexts, in particular in the modeling of the electric powergrid. A functional grid corresponds to a stable steady state, such that any bifurcation can have catastrophicconsequences up to a blackout. But also the existence of multiple steady states is undesirable, as it can leadto transitions or circulatory flows. Despite the high practical importance there is still no general theory ofthe existence and uniqueness of steady states in such systems. Analytic results are mostly limited to gridswithout Ohmic losses. In this article, we introduce a method to systematically construct the solutions of thereal power load-flow equations in the presence of Ohmic losses and explicitly compute them for tree and ringnetworks. We investigate different mechanisms leading to multistability and discuss the impact of Ohmiclosses on the existence of solutions.
This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. Thefollowing article appeared in Chaos , 123119 (2019) and may be found at https://doi.org/10.1063/1.5122739 . The stable operation of the electric power gridrelies on a precisely synchronized state of all gen-erators and machines. All machines rotate at ex-actly the same frequency with fixed phase differ-ences leading to steady power flows throughoutthe grid. Whether such a steady state exists fora given network is of eminent practical impor-tance. The loss of a steady state typically leadsto power outages up to a complete blackout. Butalso the existence of multiple steady states is un-desirable, as it can lead to sudden transitions, cir-culating flows and eventually also to power out-ages. Steady states are typically calculated nu-merically, but this approach gives only limitedinsight into the existence and (non-)uniquenessof steady states. Analytic results are availableonly for special network configurations, in par-ticular for grids with negligible Ohmic losses orradial networks without any loops. In this ar-ticle, we introduce a method to systematicallyconstruct the solutions of the real power load-flow equations in the presence of Ohmic losses.We calculate the steady states explicitly for ele-mentary networks demonstrating different mech-anisms leading to multistability. Our results alsoapply to models of coupled oscillators which arewidely used in theoretical physics and mathemat-ical biology. a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected]
I. INTRODUCTION
The electric power grid is one of the largest man-madesystems, and a stably operating grid is integral for theentire economy, industry, and almost all other technicalinfrastructures. The complexity of the power grid withthousands of generators, substations and transmission el-ements calls for an interdisciplinary approach to ensurestability in a transforming energy system . In partic-ular, the interrelation of structure and stability of com-plex grids has received widespread attention in recentyears, see e.g. . These endeavours have been aidedby the similarity of mathematical models across scien-tific disciplines. The fundamental models for power griddynamics such as the classical model or the structure-preserving model are mathematically equivalent tothe celebrated Kuramoto model with inertia . There-fore, results obtained on networks of Kuramoto oscil-lators can be easily translated to power grids and viceversa.A central question across disciplines is whether a sta-ble steady state exists and whether it is unique given acertain network structure. In the context of power grids,it is desirable to have a unique steady state. Grid opera-tors strive to maintain the flows across each line below acertain limit to avoid disruptions. Ensuring this is muchmore difficult if one has to take into account multiplesteady states, and hence multiple flow patterns acrossthe lines. Analytic results have been obtained for vari-ous special cases. In particular, multistability has beenruled out for lossless grids in the two limiting cases ofvery densely connected networks as well as tree-likenetworks (very sparse) . The existence of a steady stateis determined by two factors: the distribution of the realpower injections (natural frequencies for Kuramoto oscil-lators) and the strength of connecting lines. A variety ofrelated results have been obtained for tree-like distribu- a r X i v : . [ n li n . AO ] J a n tion grids in power engineering, see e.g. .The situation is more involved for networks of inter-mediate sparsity such as power transmission grids, whichcan give rise to multistability . The existence ofmultiple steady states in meshed networks can be tracedback to the existence of cycle flows that do not affect thepower balance at any node in the grid. The number ofand size of the cycles in the grid is thus an essential factorthat determines the number of steady states . Explor-ing the quantitative relationship between these topolog-ical factors and multistability, rigorous bounds on thenumber of steady states and mechanisms for a grid toswitch from one steady state to another one have beenfound .Despite the great theoretical progress, a general the-ory of the solvability of the power flow equations isstill lacking. Most analytic studies focus on losslessgrids or tree-like grids . Ana-lytic results are extremely rare for the full power flowequations with Ohmic losses in meshed networks .In this article, we present a new approach to computethe steady states of the real power flow equations in gen-eral networks in the presence of Ohmic losses, extendinga prior study of lossless grids . Our main contributionis a stepwise procedure to construct solutions. In a firststep, flows and losses are treated as independent vari-ables, turning the load flow equations into a linear set ofequations. The inherent relationship between flows andlosses is reintroduced in a second step. Choosing an ap-propriate basis for the solution space of the linear set ofequations, we can explicitly compute the coefficients thatlead to a consistent solution. Using this approach, weshow that Ohmic losses in general have two contrary ef-fects on the solvability of the real power flow equations:On the one hand, increasing losses requires higher linecapacities to be able to transport the same amount ofpower thereby potentially destabilizing the grid and thuslosing stable fixed points. On the other hand, we demon-strate for two very basic topologies that high line lossesmay also cause multistability leading to additional stablefixed points through a mechanism non-existent for thelossless case.The article is organized as follows: we first specify themathematical structure of the problem and fix the no-tation in section II. We then briefly review the losslesscase in section III to illustrate the fundamental impor-tance of cycles and cycle flows. Section IV constitutesthe main part of the paper, introducing the stepwise ap-proach. We then investigate two topologies in detail: atree and a ring network, for which we lay down the pro-cedures for computing all the steady states, in sectionsV and VI. II. STEADY STATES IN POWER GRIDS ANDOSCILLATOR NETWORKS
The load-flow equations constitute the fundamentalmodel to describe the steady state of an AC power grid.The system state is defined in terms of the magnitude andphase of the nodal voltages V j e iθ j , j ∈ { , . . . , N } , thathave to satisfy the energy conservation law. The nodesprovide or consume a certain amount of real power P j such that the real power balance reads P j = (cid:88) k (cid:16) b jk V j V k sin( θ j − θ k )+ g jk (cid:0) V j − V j V k cos( θ j − θ k ) (cid:1) (cid:17) . (1)Here, g jk is the conductance of the line ( j, k ) , while thesusceptance is given by − b jk (not + b jk !). By this defini-tion both g jk and b jk are generally positive for all trans-mission elements, with g jk = b jk = 0 if the two nodes j and k are not connected. The variation of the voltagemagnitudes V j is intimately related with the provisionand demand for reactive power. In general, generatornodes adapt the reactive power to fix the voltage to thereference level V j = V ref , while load nodes consume afixed value of reactive power. The voltage magnitude V j can depart from the reference level , but strict secu-rity rules are imposed to limit this voltage variation. Inthe present article we will focus on the real power bal-ance equation (1) to explore the existence of solutionsand possible routes to multistability. We neglect volt-age variability to reduce the complexity of the problemand refer to for a detailed discussion of this issue.Technically, this corresponds to the assumption that thereactive power can be balanced at all nodes. Using ap-propriate units, referred to as the pu system in powerengineering we can thus set V j = V ref = 1 for all nodes.The network structure plays a decisive role for the ex-istence and stability of steady states. This structure isencoded in the coupling coefficients b and g . For a giventransmission line ( j, k ) with resistance r jk and reactance x jk we have g jk − ib jk = 1 r jk + ix jk . (2)In high voltage transmission grids, Ohmic losses are typi-cally small such that g is small compared to b . In the limitof a lossless line, we obtain g jk = 0 and b jk = 1 /x jk > .In contrast, b and g are of similar magnitude in distribu-tion grids.A mathematically equivalent problem arises in theanalysis of steady states of dynamical power system mod-els. In particular, the dynamics of coupled synchronousmachines is determined by the swing equation M j d θ j dt + D j dθ j dt = P j − P el j , (3)whose steady states are again determined by Eq. (1).Furthermore, coupled oscillator models are used to de-scribe the collective motion of various systems across sci-entific disciplines. For instance, the celebrated Kuramotomodel considers a set of N limit cycle oscillators whosestate is described by their phases θ j along the cycle. Inmany important applications , the equations of mo-tions of the coupled system are given by dθ j dt = ω j + N (cid:88) k =1 K jk sin( θ k − θ j + γ jk ) , (4)where ω j is the intrinsic frequency of the j -th oscillator, K jk = K kj is the coupling strength of the link betweenoscillators j and k and γ jk = γ kj is a phase shift. Thefixed points of this model are determined by the algebraicequations dθ j /dt = 0 that are cast into the following formby using basic trigonometric identities ω j + (cid:88) k K jk sin( γ jk ) = (cid:88) k (cid:16) K jk cos( γ jk ) sin( θ ∗ j − θ ∗ k )+ K jk sin( γ jk ) (cid:2) − cos( θ ∗ j − θ ∗ k ) (cid:3) (cid:17) , (5)where (cid:126)θ ∗ = ( θ ∗ , . . . , θ ∗ N ) is a fixed point. This equa-tion is identical to the real power balance (1) if we iden-tify P in j = ω j + (cid:80) k K jk sin( γ jk ) , b jk = K jk cos( γ jk ) and g jk = K jk sin( γ jk ) . We note that in the limit of a losslessline, γ jk = 0 for all edges. In the following, we will fix aslack node s that can provide an infinite amount of power P s , which translates as an additional free parameter tothe Kuramoto model given by the frequency at the nodecorresponding to the slack node ω s . Therefore, differentfixed points, i.e., solutions to Eq. (5), can have a differ-ent frequency at the slack node ω s in this set-up, whichdiffers from the way fixed points are typically consideredin the Kuramoto model.The stability of a given fixed point (cid:126)θ ∗ is assessed byadding a small perturbation and then using linear sta-bility analysis, θ j = θ ∗ j + ξ j , j = 1 , . . . , N. (6)For the first order model, the dynamics of the perturba-tion is to linear order given by dξ j dt = N (cid:88) k =1 w jk ( ξ k − ξ j ) with the weights w jk = K jk cos( θ ∗ k − θ ∗ j + γ jk )= b jk cos( θ ∗ k − θ ∗ j ) − g jk sin( θ ∗ k − θ ∗ j ) . This relation is expressed in vectorial form as d(cid:126)ξdt = − Λ (cid:126)ξ (7) with the Laplacian matrix Λ ∈ R N × N with elements Λ jk = (cid:26) − w jk for j (cid:54) = k (cid:80) (cid:96) w j(cid:96) for j = k. (8)Before we proceed we note that Λ always has a zeroeigenvalue corresponding to a global shift of all phases θ j → θ j + c that does not affect the synchronization ofthe system. We thus discard this mode and limit thestability analysis to the subspace perpendicular to it D ⊥ = { (cid:126)y ∈ R N | (1 , , . . . , (cid:126)y (cid:62) = 0 } . (9)A steady state is linearly stable if all perturbations in D ⊥ are damped exponentially, which is the case if the realpart of all eigenvalues of Λ are strictly positive (exceptfor the zero eigenvalue corresponding to a global phaseshift).Stability analysis becomes rather simple in the losslesscase. Assuming that the network is connected and thatthe phase differences along any line are limited as | θ ∗ k − θ ∗ j | < π , (10)the matrix Λ is a proper graph Laplacian of an undi-rected graph, whose relevant eigenvalues are always pos-itive. Hence, Eq. (10) is a sufficient condition for linearstability but not a necessary one. Stable steady statesthat violate condition (10) do exist at the boundary ofthe stability region, but in most cases states with phasedifferences that are this large are unstable . Hence,we typically focus on states that do satisfy (10) and referto this as the normal operation of the grid .The stability analysis is more involved in the presenceof Ohmic losses, as Λ is no longer symmetric. Hence,it rather corresponds to the Laplacian of a directed net-work, whose stability is harder to grasp analytically. Inthis case we will evaluate the linear stability of differentsteady states by direct numerical computations.However, in the case where all off-diagonal elementsof this matrix are strictly negative, we are able to gainlimited analytical insight by the following Lemma: Lemma 1.
Let (cid:126)θ ∗ ∈ R N be an equilibrium of the Ku-ramoto model with phase lags as defined in Eq. (4) . Theequilibrium is linearly stable if all edges ( j, k ) have posi-tive weights w jk = K jk cos( θ ∗ k − θ ∗ j + γ jk ) > , ∀ ( j, k ) . A proof is given in Appendix A. Note that similar re-sults have also been reported in Ref. . III. THE LOSSLESS CASE
We briefly review the analysis of the lossless case tointroduce the fundamentals of our approach as well assome notation and methodology following Ref. . A. Constructing solutions
Consider a graph G consisting of N nodes and M edges. The lossless case is recovered from equation (1)by putting g jk = 0 and assuming V j = V ref = 1 , ∀ j . Thesteady states are then determined by the equation (cid:126)P = IB d sin( I (cid:62) (cid:126)θ ) . (11)Here, the sine function is assumed to be taken element-wise and we summarized all quantities in a vectorial form (cid:126)P = ( P , . . . , P N ) (cid:62) ∈ R N ,(cid:126)θ = ( θ , . . . , θ N ) (cid:62) ∈ R N , B d = diag( b , . . . , b M ) ∈ R M × M . The topology of the network is encoded in the node-edgeincidence matrix I ∈ R N × M with elements I j,e = +1 if node j is the tail of edge e ˆ= ( j, (cid:96) ) , − if node j is the head of edge e ˆ= ( j, (cid:96) ) , otherwise . (12)Based on this matrix, we also fix an orientation for eachof the network’s edges . Steady states exists only ifthe power injections of the entire grid are balanced, i.e., (cid:80) j P j = 0 , which we assume to hold.The main idea to construct all solutions of Eq. (11)is to shift the focus from nodal quantities to edges andcycles of the network. To do so, we define a vector (cid:126)F =( F , . . . , F M ) (cid:62) ∈ R M of flows on the network’s edges (cid:126)F = B d sin( I (cid:62) (cid:126)θ ) . (13)If a component of the flow vector is larger than zero, F e > , the flow on link e = ( k, j ) is directed from k to j and if F e < from j to k . Therefore F e physicallydenotes the flow from the tail of the edge e to the headof e . Eq. (11) then becomes (cid:126)P = I (cid:126)F . (14)Solutions of Eq. (11) may be constructed by first solv-ing Eq. (14) and then rejecting all solution candidateswhich are incompatible with Eq. (13). Solutions of (14)may be obtained based on the following observation: thekernel of the incidence matrix I corresponds exactly to cycle flows , a cycle flow referring to a constant flowalong a cycle with no in- or out-flow . The kernelhas dimension M − N + 1 , which reflects the fact thatthe cycles in a graph forms a vector space of dimension M − N + 1 , a basis set of this space is called a fun-damental cycle basis. A set of fundamental cycles B isencoded in the corresponding cycle-edge incidence matrix C B ∈ R M × ( M − N +1) with elements C Be,c = +1 if the edge e is part of the cycle c − if the reversed edge e is part of cycle c otherwise . (15) Then, all solutions of equation (14) can be written as (cid:126)F = (cid:126)F ( s ) + C B (cid:126)f , (16)where (cid:126)F ( s ) ∈ R M is a specific solution and (cid:126)f ∈ R M − N +1 gives the strength of the cycle flows along each cycle inthe chosen cycle basis. Having obtained a flow vector (cid:126)F ,we can simply construct the associated phases as follows.Start at the slack node s and set θ s = 0 . Then proceedto a neighbouring node j . Assuming that the connectingedge e ˆ=( j, s ) is oriented from node s to node j , the phasevalue reads θ j = θ s + ∆ e , (17)where the phase difference ∆ e is reconstructed from theflow F e by inverting Eq. (13), ∆ + e = arcsin( F e /b e ) or ∆ − e = π − arcsin( F e /b e ) . (18)For each edge e we have to decide whether we take the + -solution or the − -solution in Eq. (18). To keep trackof this choice, we decompose the edge set of the network E into two parts, E + = { e ∈ E | ∆ e = ∆ + e } E − = { e ∈ E | ∆ e = ∆ − e } , such that E = E + ∪ E − . Not all solutions obtained thisway are physically correct. We can obtain the physicallycorrect ones by making sure that the sum of the phase dif-ferences around any fundamental cycle yields zero or aninteger multiple of π , which is equivalent to the windingnumbers (cid:36) c = 12 π M (cid:88) e =1 C Be,c ∆ ± e , (19)summarized in the vector (cid:126)(cid:36) = ( (cid:36) , . . . , (cid:36) M − N +1 ) (cid:62) be-ing integer (cid:36) c ∈ Z . It should be noted that the choice ∆ + e corresponds to the state of normal operation discussed insection II. Hence, states with E − = ∅ are guaranteed tobe stable, while states with E − (cid:54) = ∅ are typically (butnot always) unstable . We summarize these resultsin the following proposition due to Ref. . Proposition 1.
Consider a connected lossless networkwith power injections (cid:126)P ∈ R N . Then the following twostatements are equivalent:1. (cid:126)θ is a steady state, i.e., a real solution of equation(11).2. The flows (cid:126)F ∈ R M satisfy the ‘dynamic’ conditions(14) with | F e | ≤ b e such that (cid:126)F = (cid:126)F ( s ) + C B (cid:126)f (20) and the geometric condition (19) (cid:126)(cid:36) ( (cid:126)f ) ∈ Z M − N +1 . (21) for some decomposition E = E + ∪ E − . IV. POWER GRIDS WITH OHMIC LOSSES
We now extend the approach introduced above topower grids with Ohmic losses or oscillator networkswith a general trigonometric coupling. The steadystates are determined by the real power balance equa-tion (cf. Eq. (1)) P j = N (cid:88) k =1 (cid:16) b jk sin( θ j − θ k )+ g jk [1 − cos( θ j − θ k )] (cid:17) . (22)Before we proceed to construct the solution to these equa-tions we note an important difference to the lossless case.The Ohmic losses occurring on the lines are not a prioriknown as they depend on the phases θ , . . . , θ N . Hencethe real power balance for the entire grid now reads N (cid:88) j =1 P j = P losses ( θ , . . . , θ N ) . (23)Thus for arbitrary P , . . . , P N , there will typically be nosolution. This issue is solved by assuming that one ofthe nodes, referred to as the slack node, can provide anarbitrary amount of power to balance the losses. For thesake of consistency, we label the slack as j = 1 through-out this article and set θ = 0 . We note that the choice ofa particular slack node is often arbitrary. In transmissiongrids one typically chooses a node with high generation,whereas in distribution grids one can choose the connec-tion to the higher grid level. Other approaches using adistributed slack bus also exist, see e.g. .To solve the set of equations (22) for the remainingnodes j ∈ { , . . . , N } we decompose it into different partsas before and first formulate a linear system of equa-tions. Before we start, we fix some notations by definingthe unsigned incidence matrix E ∈ R N × M with elements E j,e = | I j,e | . For each edge e ˆ=( j, k ) we now define thelosses by L e = g e [1 − cos( θ j − θ k )] . Using this notation, the power balance equations canbe decomposed into three parts. First we have the dy-namic condition, which now reads ( Ia ) P j = M (cid:88) e =1 I j,e F e + E j,e L e , ∀ j ∈ { , . . . , N } . (24)Flows and losses are limited by the line parameters whichis represented by the following conditions ( Ib ) F e b e ∈ [ − , , L e g e ∈ [0 , , ∀ e ∈ { , . . . , M } . (25)In addition to that, flows and losses are not independent,but are both functions of the phase difference θ j − θ k . Us-ing the trigonometric identity sin + cos = 1 we obtain the flow-loss condition ( II ) (cid:18) F e b e (cid:19) + (cid:18) L e g e − (cid:19) = 1 , ∀ e ∈ { , . . . , M } . (26)Finally, we have a geometric condition as in the losslesscase ( III ) (cid:36) c ( F , . . . , L , . . . ) = z with z ∈ Z , ∀ cycles c. (27)In comparison to the lossless case we have M additionaldegrees of freedom L , . . . , L M and M additional nonlin-ear conditions (26) to fix them. Furthermore, the knowl-edge of both F e and L e is sufficient to fix the phasescompletely. Eq. (18) is replaced by ∆ e = (cid:26) arcsin( F e /b e ) if L e ≤ g e π − arcsin( F e /b e ) if L e > g e . (28)Still, there are two solution branches ± per edge as in thelossless case, because the quadratic equation (26) has twosolutions in general. We summarize these findings in thefollowing proposition. Proposition 2.
Consider a connected lossy network withpower injections (cid:126)P ∈ R N . Then the following two state-ments are equivalent:1. (cid:126)θ is a steady state, i.e., a real solution of equation(22).2. The flows (cid:126)F ∈ R M and losses (cid:126)L ∈ R M satisfythe ‘dynamic’ conditions (24) with | F e | ≤ b e and ≤ L e ≤ g e , the flow loss condition (26), and thegeometric condition (cid:126)(cid:36) ∈ Z M − N +1 . (29)To find actual solutions, we thus have to solve the lin-ear set of equations (24) subject to a variety of nonlin-ear constraints (25-27). Remarkably, this can be accom-plished in an iterative fashion such that we find the fol-lowing general strategy to construct solutions:1. Construct the solution space of the linear set ofequations (24), which yields a potentially large setof solution candidates . This set is gradually re-duced in the further steps until only the actual so-lutions are left.2. Use the flow-loss condition (26) to reduce the de-grees of freedom of the system. In particular, allremaining solution candidates can be expressed interms of the cycle flow strengths and a set of in-dices ± which indicate the solution branch for eachedge.3. Finally, fix the cycle flows by the geometric condi-tions (27). FIG. 1. Labeling of nodes (blue circles) and edges (blackarrows) in a tree network used in Sec. V A. The slack node istaken as the root of the tree and labeled as j = 1 as indicatedby the letter S and the darker blue colouring. Remarkably, we will see in the following that condi-tion (25) on the line limits is automatically satisfied, ifa real solution of the flow-loss condition (26) exists, sowe do not have to explicitly consider this condition (seeLemma 2). We further note that the addition of cycleflows still does not affect the power balance, so the cycleflows remain a basic degrees of freedom in equation (24).The losses L e are fixed only in the second step using theflow-loss condition (26). Hence, the resulting losses de-pend on the strength of cycle flows. We now illustratethis approach by explicitly constructing the solutions fora tree network and a single cycle. We will show that in-cluding losses gives rise to an additional mechanism ofmultistability. V. TREE NETWORKS
We will first consider tree networks, i.e. networks with-out any closed cycles. Hence, we do not have to take intoaccount the geometric condition (27) and focus on thesolution of the flow-loss condition (26).
A. Fundamentals
We first introduce the basic notation, see Fig. 1. Theslack node is chosen to be the root of the tree and labeledas j = 1 . The remaining nodes are labeled according tothe distance to the root: first nearest neighbors, thennext-to-nearest neighbors, and so on. Every edge e =1 , . . . , M = N − points to the node e + 1 . For eachnode and edge, we must keep track of how it is connectedto the root of the tree. We thus introduce the matrix T ∈ R M × M by T e,j = +1 if edge e is on the path from node j + 1 to the root otherwise . Note that the edges are labeled in such a way that T e,k also indicates whether edge e is on the path from edge k to the root. Furthermore, we introduce the vectorialnotation (cid:126)F = ( F , . . . , F M ) (cid:62) ,(cid:126)L = ( L , . . . , L M ) (cid:62) ,(cid:126)x = ( F , . . . , F M , L , . . . , L M ) (cid:62) . The dynamic condition (24) then reads (cid:126)P = I (cid:126)x, (30)where the matrix I ∈ R ( N − × M is obtained by concate-nating the signed and unsigned incidence matrix ( I | E ) and removing the first line corresponding to the slacknode. In particular, the matrix elements are given by I j − ,e = +1 if e ≤ M and j is the tail of edge e or if e > M and j is the tail orhead of edge e − M − if e ≤ M and j is the head of edge e otherwise . (31)First, we need a specific solution (cid:126)x ( s ) of the dynamiccondition (30). For the sake of simplicity, we choose asolution with no losses, that is (cid:126)x ( s ) = ( (cid:126)F ( s )1 , . . . , (cid:126)F ( s ) M , , . . . (cid:62) , (32)where F ( s ) e = − N (cid:88) j =2 T e,j − P j . (33)Then we have to construct the general solution to thedynamic conditions, i.e., we need a basis for the N -dimensional kernel of the matrix I . The basis vectorsare constructed such that they have losses only at oneparticular line, which yields (cid:126)x ( e ) = (cid:20) (cid:126)F ( e ) (cid:126)L ( e ) (cid:21) , ∀ e ∈ { , . . . , M } F ( e ) k = 2 T e,k + δ e,k L ( e ) k = δ e,k with the Kronecker symbol δ e,k . This set of basis vectorsis illustrated in Fig. 2 for an elementary example. Wenote that these basis vectors are linearly independent asrequired, but not orthogonal. All solution candidates of a bc d FIG. 2. (a) Simple tree network with N = 4 nodes M =3 edges. Arrows indicate the orientation of edges which inturn determines the direction of flows. (b-d) Illustration ofthe basis vectors of the kernel of the matrix I . The vectors (cid:126)x e , e = 1 , , include losses at exactly one edge e , indicatedby the dotted red arrows at the terminal nodes, and the flowsneeded to compensate these losses. the dynamic and the flow-loss conditions can be writtenas (cid:126)x = (cid:126)x ( s ) + M (cid:88) e =1 α e (cid:126)x ( e ) , (34)In terms of the flows and losses this yields F e = F ( s ) e + 2 N (cid:88) k = e +1 T e,k α k + α e ,L e = α e . (35)To simplify the notation, we introduce the abbreviation F e = − N (cid:88) j =2 T e,j − P j + 2 N (cid:88) k = e +1 T e,k α k , (36)which is the flow on the line e minus the losses, F e = F e − L e = F e − α e . Now we can calculate the parameters α e by substitut-ing ansatz (35) into the flow-loss condition (26): (cid:18) F e + α e b e (cid:19) + (cid:18) α e g e − (cid:19) = 1 . (37)To solve these quadratic equations we now have to pro-ceed iteratively from e = N − to e = 1 as the quan-tity F e depends on the losses α k on the lines k = e + 1 , . . . , N − . In each step, we have to check whetherthe solutions are real, positive and respect the line limits(25). Fortunately, these conditions can be simplified toa single inequality condition as stated in the followinglemma. Lemma 2.
Eq. (37) has two real positive solutions α ± e which both satisfy the line limits (25) if and only if b e ≥ F e + 2 g e F e . (38) The two solutions coalesce in the case of equality.
We emphasize that condition (38) has to be satisfied forall edges e ∈ { , . . . , M } , which again has to be verifiediteratively. A proof of this result is given in Appendix B.Finally, we can summarize our findings as follows. Lemma 3.
All potential solutions of the dynamic condi-tions and the load-flow condition for a tree network canbe written as F e = − N (cid:88) j =2 T e,j − P j + 2 N (cid:88) k = e +1 T e,k α k (cid:124) (cid:123)(cid:122) (cid:125) =: F e + α e L e = α e , where the parameters α e , e ∈ { M, M − , . . . , } are de-termined iteratively as α ± e = g e b e ( g e + b e ) (cid:20) b e − g e b e F e − σ e (cid:112) b e − F e − g e F e (cid:21) , (39) where the sign σ e ∈ {− , +1 } indicates the solutionbranch. Hence, each potential solution is uniquely char-acterized by the sign vector (cid:126)σ = ( σ , . . . , σ M ) (cid:62) ∈{− , +1 } M . B. Example
As an example we consider a grid with N = 4 nodesand M = 3 edges as depicted in Fig. 2 (a). The node-edgeincidence matrix I and its modulus E read I = +1 0 0 − − − ⇒ E = +1 0 0+1 +1 +10 +1 00 0 +1 , and the tree matrix is given by T = +1 +1 +10 +1 00 0 +1 . The dynamic condition (24) thus reads P = − F + F + F + L + L + L P = − F − L P = − F − L . g a + 0 10 20 30 40 50 g b + 0 2 4 6 8 g c + ++ + ++ ++ ++ + + FIG. 3. Multiple solutions in a tree network with Ohmic losses. The possible values of the line losses α e are shown as a functionof the conductance g for the simple four node tree network shown in Fig. 2 and parameters b = 10 , P = − , P = − and P = − for varying g as calculated according to Eq. (44). Solid, coloured lines indicate dynamically stable solutions anddotted, coloured lines indicate unstable ones. The black dotted line indicates points with α e = g , thus determining whichbranch to choose when calculating angular differences according to Eq. (28). (a,b) Branching of α and α into two differentsolutions according to the different signs of the square root in the expression (Eq. (44)). (c) The solutions found for α and α can be used to subsequently calculate the solutions for α . The solutions depend on the signs σ e for all lines e = 1 , , suchthat we find solution branches in total. The signs indicated in the legend are ordered as ( σ , σ , σ ) . In the region shadedin grey, there are two coexisting stable solutions. A particular solution of these equations is given by (cid:126)x ( s ) = (cid:20) (cid:126)F ( s ) (cid:126)L ( s ) (cid:21) = ( − P − P − P , − P , − P , , , (cid:62) . and the kernel is spanned by the basis vectors (cid:126)x = (1 , , , , , (cid:62) ,(cid:126)x = (2 , , , , , (cid:62) ,(cid:126)x = (2 , , , , , (cid:62) , which are illustrated in Fig. 2 (b-d). Hence, the generalsolution can be written as (cid:126)x = F F F L L L = F ( S )1 + 2 α + 2 α + α F ( S )2 + α F ( S )3 + α α α α . The coefficients α e , e ∈ { , , } , are directly calculatedin the order e = 3 , , via Eq. (39) with F = F ( s )3 , F = F ( s )2 and F = F ( s )1 + 2 α ± + 2 α ± . We recall that incontrast to the cyclic case we do not have to consider thegeometric condition. The values of α ± e and hence alsothe flows and losses depend only on the signs ( σ , σ , σ ) – and of course on the system parameters.To explore the emergence of multistability in networkswith Ohmic losses, we plot the different solution branches as a function of the conductances g in Fig. 3. For thesake of simplicity we assume that all lines have the sameparameters, and keep both b and the power injectionsfixed. In the limit g → , we trivially have α e = 0 forall edges such that the functions α + e and α − e coalesce.However, this does not imply that equilibria themselvescoalesce, cf. Eq. (18). For small values of g , the line losses α e then increase approximately linearly and we find different solutions in total, corresponding to the differentchoices of the signs ( σ , σ , σ ) . For each edge, the + branch corresponds to a solution with low losses L e < g e and the − branch to a solution with high losses L e > g e .Nonlinear effects become important for higher values of g : The losses in the + branches increase super-linearly,while the − branches show a non-monotonic behaviour.For even higher values of g solutions vanish pairwise. Thesolution branches (cid:126)σ = (+ , + , +) with the lowest overalllosses and the branch (cid:126)σ = (+ , + , − ) vanishes last.We further evaluate the dynamical stability for eachsolution branch by numerically testing the eigenvalues ofthe matrix Λ defined in (8). The weights used in thisLaplacian matrix can be rewritten directly in terms ofthe flows and losses. If nodes j and k are connected viaedge e , we obtain w jk = b e g e ( g e − L e ) ± g e b e F e , where the minus sign is chosen if j is the tail and k thehead of edge e and the plus sign is chosen if k is the tailand j the head of edge e . b g FIG. 4. Number of stable fixed points (colour code) of thelossy real power flow equation 1 for the four node tree networkshown in Fig. 2(a) with power injections P = − , P = − and P = − for varying line susceptance b (abscissa) andconductance g (ordinate). Whereas a minimum line capacityis required to result in any stable fixed points in the same wayas for the lossless power flow, two effects that do not exist inthe lossless case may be observed: Increasing conductances g and thus losses requires for higher line capacities b as ex-pected. In addition to that, an additional stable fixed pointarises for higher losses thus presenting a different mechanismfor multistability. The results for the stability of the different solutionbranches are indicated by displaying the lines as eitherdashed (unstable) or solid (stable) in Fig. 3. We find thatonly the (+ , + , +) -branch is stable for low losses. This isexpected because in the lossless case there can be at mostone stable solution . The (+ , + , +) -branch continuouslymerges into this stable solution in the limit g → . Moreinterestingly, also the (+ , + , − ) -branch becomes stablefor large values of g . Hence, losses can stabilize fixedpoints.A comprehensive analysis of the existence of solutionsfor the given sample network in terms of the grid pa-rameters b and g is given in Fig. 4. Remarkably, thepresence of Ohmic losses has two antithetic effects on thesolvability of the real power load-flow equations: On theone hand, losses can prohibit the existence of solutions.Real power flows are generally higher in lossy networksas losses have to be balanced by additional flows. Hence,the minimum line strength b required for the existenceof a solution increases with g . On the other hand, lossesfacilitate multistability. While the lossless equation canhave at most one stable fixed point for tree networks, twostable fixed points can exist if losses are added.For example, for three consumer nodes with power in-jections P = − , P = − and P = − and uniformline parameters of b = 10 and g = 8 , we find a dynami-cally stable solution branch with (cid:126)σ = (+ , + , − ) with flows FIG. 5. Labeling of nodes (blue circles) and edges (blackarrows) in a cyclic network used in Sec. VI A. The slack nodeis located at j = 1 and indicated here by the letter S and acolouring in darker blue. (cid:126)F ≈ (9 . , . , . (cid:62) and losses (cid:126)L ≈ (4 . , . , . (cid:62) and another one with (cid:126)σ = (+ , + , +) with flows (cid:126)F ≈ (6 . , . , . (cid:62) and losses (cid:126)L ≈ (1 . , . , . (cid:62) . We re-call that node 1 serves as a slack node. Hence, the powerinjection P (or the natural frequency ω in the oscillatorcontext) is different for the two stable steady states. VI. CYCLIC NETWORKA. Fundamentals
We now consider a single closed cycle as depicted inFig. 5. We label all nodes by j ∈ { , . . . , N } around thecycle in the mathematically positive direction startingat the slack node j = 1 . Similarly, we label all lines e ∈ { , . . . , N } where line e corresponds to ( e, e + 1) andline e = N corresponds to ( N, .We now construct the solutions of the dynamic condi-tion (30). As before, we choose a specific solution withno losses (cf. Eq. (32)), where the flows satisfy P j = N (cid:88) e =1 I j,e F ( s ) e , ∀ j ∈ { , . . . , N } . A solutions always exists as the linear set of equationshas rank N − . A proper initial guess can be obtained,for example, by solving the DC approximation .To construct the general solution, we further need abasis for the ( N + 1) -dimensional kernel of the matrix I .As before, we use a set of basis vectors that have losses0 a bc d FIG. 6. Illustration of the basis vectors of the kernel of thematrix I for a small cyclic network with N = 3 nodes. (a-c)The vectors (cid:126)x e , e = 1 , . . . , N , include losses at exactly oneedge e , indicated by the dotted red arrows at the terminalnodes, and the flows needed to compensate this loss. (d) Thebasis vector (cid:126)x N +1 represents a lossless cycle flow. only at one particular line, (cid:126)x ( e ) = (2 , . . . , (cid:124) (cid:123)(cid:122) (cid:125) e − , , , . . . , (cid:124) (cid:123)(cid:122) (cid:125) N − e times , , . . . , (cid:124) (cid:123)(cid:122) (cid:125) e − , , , . . . , (cid:124) (cid:123)(cid:122) (cid:125) N − e times ) (cid:62) . (40)In contrast to the tree network we need an additionalbasis vector describing a cycle flow (cid:126)x ( N +1) = (1 , , . . . , (cid:124) (cid:123)(cid:122) (cid:125) N times , . . . , (cid:62) . (41)This set of basis vectors is illustrated in Fig. 6. Allsolution candidates of the dynamic and the flow-loss con-ditions can thus be written as (cid:126)x = (cid:126)x ( s ) + f(cid:126)x ( N +1) + N (cid:88) e =1 α e (cid:126)x ( e ) , (42)where f ∈ R is a parameter giving the cycle flow strength.In terms of the flows and losses this yields F e = F ( s ) e + f + 2 N (cid:88) n = e +1 α n (cid:124) (cid:123)(cid:122) (cid:125) =: F e + α e ,L e = α e . (43)As before, we can now calculate the parameters α e iteratively from e = N to e = 1 using Eq. (39) α ± e = g e b e ( g e + b e ) (cid:20) b e − g e b e F e − σ e (cid:112) b e − F e − g e F e (cid:21) . (44) a b FIG. 7. Simple cycle network with three nodes (blue circles)and three edges (black arrows). (a) Arrows indicate the orien-tation of edges which in turn determines the direction of flows.We consider a network with power injections at the nodes P , P and P and power flows on the edges denoted F , F and F . (b) Example studied in section VI B. The node j = 1 ischosen as a slack node and the (indicated by symbol S ) andthe two other nodes are assumed to be consumer nodes with P , = − P . Arrows again represent the edge orientations andthe values give the specific solution F ( s )1 = P , F ( s )2 = 0 and F ( s )3 = − P . However, we now have to take into account that the quan-tities F e also depend on the parameter f – the cycleflow strength. Hence, each potential solutions is nowcharacterized by the continuous parameter f in additionto the signs σ , . . . , σ N ∈ {− , +1 } . Whether a solu-tion exists and respects the line limits can be determinedfrom Lemma 2, in particular from condition (38). Westress that this condition must be satisfied for all edges e ∈ { , . . . , N } simultaneously, keeping in mind that thequantities F e depend on the values F e +1 , . . . , F N andthe cycle flow strength f . Hence, condition (38) must bechecked iteratively for all e = N, N − , . . . , in depen-dence of the value of f .In a cyclic network we further have to satisfy the geo-metric condition (27), which fixes the remaining contin-uous degree of freedom f . For a single cycle, the windingnumber is given by (cid:36) (cid:126)σ = 12 π M (cid:88) e =1 ∆ σ e e , The phase differences ∆ σ e e and hence the winding numberare determined by the line flows and losses via Eq. (28)and depend on the respective solution branch indicatedby the signs (cid:126)σ . Recall that the geometric condition statesthat the winding number (cid:36) can be an arbitrary integer.Hence there can be multiple solutions for f for a givenset of signs σ , . . . , σ N if the cycle is large enough. Thisroute to multistability was analyzed in detail for losslessnetworks in .1 f a + 4 2 0 2 4 f b + ++ + 4 2 0 2 4 f c + + ++ ++ +++++ + FIG. 8. The possible values of line losses α , , as a function of the cycle flow strength f for the simple three node cycle networkshown in Fig. 7 and parameters P = − , g = 1 and b = 4 . Black dotted line indicates values where α e = g , thus determiningthe sign of angular differences according to Eq. (28). (a) Branching of α into two different solutions referred to as α +3 (darkred, bottom) and α − (dark green, top) for the different signs of the square root as predicted by Eq. (44). (b-c) The solutionsfound for α can be used to subsequently calculate the solutions for α and then α . The signs indicated here in the legend areordered as ( σ ) , ( σ , σ ) and ( σ , σ , σ ) for panels (a), (b) and (c), respectively. B. Example
We analyze here a three-node cycle shown in Fig. 7where node is the slack node. The node-edge incidencematrix I and its modulus E will then be I = +1 0 − − − ⇒ E = +1 0 +1+1 +1 00 +1 +1 . The dynamic condition (24) thus reads (cid:40) P = F − F + L + L P = F − F + L + L . A particular solution of these equations is given by (cid:126)x ( s ) = (cid:20) (cid:126)F ( s ) (cid:126)L ( s ) (cid:21) = ( − P , , + P , , , (cid:62) . and the kernel is spanned by the basis vectors (cid:126)x = (1 , , , , , (cid:62) ,(cid:126)x = (2 , , , , , (cid:62) ,(cid:126)x = (2 , , , , , (cid:62) ,(cid:126)x = (1 , , , , , (cid:62) , which are illustrated in Fig. 6. Hence, the general solu-tion can be written as (cid:126)x = F F F L L L = F ( s )1 + f + 2 α + 2 α + α F ( s )2 + f + 2 α + α F ( s )3 + f + α α α α . (45) The coefficients α , , are calculated as a function of f iteratively starting from N = 3 via Eq. (44) with F = F ( s )3 + f , F = F ( s )2 + f + 2 α ± and F = F ( s )1 + f + 2 α ± +2 α ± . The results are shown in Fig. 8 (a-c) for all differentpossible realizations of the sign vector ( σ , σ , σ ) : for α we have 2 choices, then for α we have = 4 choices (twochoices for each of α and α ) and finally we have = 8 choices for α . For the sake of simplicity, we have chosen P = P = 1 in this example. Notably, all branches of thesolutions must form closed curves when plotted via theparameter f . This is due to the fact that a real solutionof the Eq. (44) can only vanish when the discriminantgoes to zero, i.e., when it collides with another branch ofthe solution.The remaining parameter f is determined by the ge-ometric condition (27). To evaluate this condition andto finally determine all steady states we plot the windingnumber (cid:36) (cid:126)σ ( f ) = 12 π M (cid:88) e =1 ∆ e as a function of f in Fig. 9. The phase differences aregiven by (cf. Eq. (28)) ∆ σ e e = (cid:26) arcsin( F e /b e ) if L e ≤ g e π − arcsin( F e /b e ) if L e > g e . They depend on the solution branch, i.e., on the valuesof the σ e and so does the winding number. For the givencyclic network we find solution branches, which have tobe considered when evaluating the geometric condition,see Fig. 9. Inspecting the winding number (cid:36) (cid:126)σ ( f ) for eachbranch, we find 2 steady states of which one is stable and2 f FIG. 9. The winding number (cid:36) as a function of the cycle flowstrength f for different solution branches in the three nodenetwork depicted in Fig. 7. Solutions require that (cid:36) ∈ Z ,cf. Eq. (27). Colour code as in Fig. 8,c for all panels. one is unstable. Again, the stable fixed point is given bythe (+ , + , +) -branch which has the lowest Ohmic losses.However, we can find two dynamically stable branchesfor higher losses as in the case for the tree network. Forexample, fixing line susceptances and conductances b = g = 3 and power injections P = P = − , we find againtwo dynamically stable branches corresponding to lowlosses (cid:126)σ = (+ , + , +) and high losses (cid:126)σ = (+ , + , − ) . VII. SUMMARY AND DISCUSSION
In this article, we studied solutions to the real powerload-flow equations in AC transmission grids of generaltopology with a special focus on the impact of Ohmiclosses. Extending our previous work , we constructed ananalytical method for computing all load flow solutions,both stable and unstable ones. We demonstrated howto explicitly compute all steady states in two elementarytest topologies: a -node tree and a -node ring.We find that analogous to the lossless case, different so-lutions exist corresponding to different winding numbers(19) along each basis cycle, as well as a choice betweentwo solution branches in each edge. The two branchescorrespond to a state with low losses and phase differ-ences on the respective edge ( + branch) and high lossesand phase difference ( − branch).We show that Ohmic losses have two conflicting effectson the existence and number of steady states: On the onehand, high losses must be compensated by higher flows.Hence, solutions may vanish due to Ohmic losses unlessthe line capacities are also increased. On the other hand,Ohmic losses can stabilize certain solution branches andthus foster multistability. In particular, we demonstratethat two grid topologies that have been proven to exhibit no multistability in the lossless case – trees and -noderings – are multistable in the lossy case for certain pa-rameter values. ACKNOWLEDGMENTS
We thank Tom Brown and Johannes Schiffer for valu-able discussions. We gratefully acknowledge supportfrom the German Federal Ministry of Education and Re-search (grant no. 03EK3055B) and the Helmholtz Asso-ciation (via the joint initiative “Energy System 2050 – AContribution of the Research Field Energy” and the grantno. VH-NG-1025). D.M. acknowledges funding from theMax Planck Society.
Appendix A: Proof of lemma 1
Proof.
The result can be proven by making use of Ger-shgorin’s circle theorem . Recall that the equilibrium islinearly stable if all the eigenvalues µ j of the Laplacian Λ have a positive real part, Re( µ j ) > , ∀ j ∈ { , ..., N − } , except for the eigenvalue µ N = 0 corresponding to aglobal phase shift. According to Gershgorin’s theorem,each eigenvalue µ j is located in a disk in the complexplane with radius R j = (cid:80) (cid:96) (cid:54) = j | Λ j,(cid:96) | centred at Λ j,j . Ifthe condition w jk > is satisfied, we have that | Λ j,(cid:96) | = − Λ j,(cid:96) . Therefore, applying Gershgorin’s theorem resultsin the following inequality | µ j − Λ j,j | ≤ (cid:88) (cid:96) (cid:54) = j | Λ j,(cid:96) | = Λ j,j . This inequality thus predicts that all eigenvalues µ j havereal part greater than or equal to zero Re( µ j ) ≥ . Now itremains to show that the eigenvalue µ N to the eigenvector (1 , ..., (cid:62) is the only zero eigenvalue. Assume that (cid:126)v ∈ R N is an eigenvector with eigenvalue µ = 0 . Assume thatthis vector has its minimum entry at position i , such that v i = min( v j ) , j ∈ { , ...N } and hence v i − v j ≤ , ∀ j .Then we arrive at Λ (cid:126)v ) i = (cid:88) j (cid:54) = i Λ ij ( v i − v j ) . Since the off-diagonal elements Λ ij are all negative bythe assumption of the lemma, it follows that the entriesof the vector at neighbouring nodes equal its minimumvalue v i = v j . We can now apply the same reasoningfor next-nearest neighbours and proceed in the same waythrough the whole network to show that v i = v j , ∀ j ∈ { , ...N } , which proofs that (cid:126)v = (1 , ..., (cid:62) is the only eigenvectorwith vanishing eigenvalue µ = 0 .3 ( , 1)( g , 0) (2 g , 1)(0, 1) r ( ) l ( ) FIG. 10. Illustration of the two function l ( α e ) and r ( α e ) usedin the proof of lemma 2 for arbitrary parameter values g = 0 . , b = 1 . and F = 1 . Appendix B: Proof of lemma 2
Proof.
We first note that if condition (38) is satisfied, thediscriminant in Eq. (39) is non-negative, such that allsolutions are real. The two solutions coalesce if the dis-criminant vanishes, i.e., if b e = F e + 2 g e F e . Conversely,if the condition (38) is not satisfied, the discriminant inEq. (39) is negative, such that no real solution exists.It remains to be shown that if a solution exists, thenit is positive and respects the line limits. To this end, werewrite the flow-loss condition (26) as (cid:18) α e g e − (cid:19) (cid:124) (cid:123)(cid:122) (cid:125) =: l ( α e ) = 1 − (cid:18) α e + F e b e (cid:19) (cid:124) (cid:123)(cid:122) (cid:125) =: r ( α e ) (B1)The two parabola l ( α e ) and r ( α e ) are illustrated inFig. 10. The left-hand side l ( α e ) is non-negative every-where with l ( α e ) ∈ [0 , if α e ∈ [0 , g e ] l ( α e ) > if α e / ∈ [0 , g e ] . The right-hand side is smaller or equal to one with r ( α e ) ∈ [0 , if α e ∈ [ − b e − F e , + b e − F e ] r ( α e ) < if α e / ∈ [ − b e − F e , + b e − F e ] . Hence, we find the necessary condition for the crossingof the two parabola as l ( α e ) = r ( α e ) ∈ [0 , ,L e = α e ∈ [0 , g e ] ,F e = F e + α e ∈ [ − b e , + b e ] . That is, if a solution α e exists, it is guaranteed to bepositive and satisfy the line limits. REFERENCES C. D. Brummitt, P. D. Hines, I. Dobson, C. Moore, and R. M.D’Souza, “Transdisciplinary electric power grid science,” Pro-ceedings of the National Academy of Sciences , 12159–12159(2013). M. Timme, L. Kocarev, and D. Witthaut, “Focus on networks,energy and the economy,” New journal of physics , 110201(2015). G. Filatrella, A. H. Nielsen, and N. F. Pedersen, “Analysis of apower grid using a kuramoto-like model,” The European PhysicalJournal B , 485–491 (2008). D. Witthaut and M. Timme, “Braess’s paradox in oscillator net-works, desynchronization and power outage,” New journal ofphysics , 083036 (2012). F. Dörfler and F. Bullo, “Synchronization and transient stabilityin power networks and nonuniform kuramoto oscillators,” SIAMJournal on Control and Optimization , 1616–1642 (2012). D. Witthaut and M. Timme, “Nonlocal failures in complex sup-ply networks by single link additions,” The European PhysicalJournal B , 377 (2013). F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in com-plex oscillator networks and smart grids,” Proceedings of the Na-tional Academy of Sciences , 2005–2010 (2013). D. Manik, D. Witthaut, B. Schäfer, M. Matthiae, A. Sorge,M. Rohden, E. Katifori, and M. Timme, “Supply networks:Instabilities without overload,” The European Physical JournalSpecial Topics , 2527–2547 (2014). J. W. Simpson-Porco, “A theory of solvability for lossless powerflow equations—part i: Fixed-point power flow,” IEEE Transac-tions on Control of Network Systems , 1361–1372 (2017). J. W. Simpson-Porco, “A theory of solvability for lossless powerflow equations—part ii: Conditions for radial networks,” IEEETransactions on Control of Network Systems , 1373–1385(2017). S. Jafarpour, E. Y. Huang, K. D. Smith, and F. Bullo, “Multi-stable synchronous power flows: From geometry to analysis andcomputation,” arXiv preprint arXiv:1901.11189 (2019). P. M. Anderson and A. A. Fouad,
Power system control andstability (John Wiley & Sons, 2008). A. R. Bergen and D. J. Hill, “A structure preserving model forpower system stability analysis,” IEEE Transactions on PowerApparatus and Systems , 25–35 (1981). Y. Kuramoto, “International symposium on mathematical prob-lems in theoretical physics,” Lecture notes in Physics , 420(1975). S. H. Strogatz, “From kuramoto to crawford: exploring the onsetof synchronization in populations of coupled oscillators,” PhysicaD: Nonlinear Phenomena , 1–20 (2000). J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, andR. Spigler, “The kuramoto model: A simple paradigm for syn-chronization phenomena,” Reviews of modern physics , 137(2005). A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou,“Synchronization in complex networks,” Physics reports , 93–153 (2008). R. Taylor, “There is no non-zero stable fixed point for dense net-works in the homogeneous kuramoto model,” Journal of PhysicsA: Mathematical and Theoretical , 055102 (2012). D. Manik, M. Timme, and D. Witthaut, “Cycle flows and mul-tistability in oscillatory networks,” Chaos: An InterdisciplinaryJournal of Nonlinear Science , 083123 (2017). H.-D. Chiang and M. E. Baran, “On the existence and unique-ness of load flow solution for radial distribution power networks,”IEEE Transactions on Circuits and Systems , 410–416 (1990). A. J. Korsak, “On the question of uniqueness of stable load-flowsolutions,” IEEE Transactions on Power Apparatus and Systems, 1093–1100 (1972). D. A. Wiley, S. H. Strogatz, and M. Girvan, “The size of the sync basin,” Chaos: An Interdisciplinary Journal of NonlinearScience , 015103 (2006). J. Ochab and P. Gora, “Synchronization of coupled oscillators in alocal one-dimensional kuramoto model,” Acta Physica Polonica.Series B, Proceedings Supplement , 453–462 (2010). R. Delabays, T. Coletta, and P. Jacquod, “Multistability ofphase-locking and topological winding numbers in locally coupledkuramoto models on single-loop networks,” Journal of Mathemat-ical Physics , 032701 (2016). R. Delabays, T. Coletta, and P. Jacquod, “Multistability ofphase-locking in equal-frequency kuramoto models on planargraphs,” Journal of Mathematical Physics , 032703 (2017). D. Mehta, N. S. Daleo, F. Dörfler, and J. D. Hauenstein, “Al-gebraic geometrization of the kuramoto model: Equilibria andstability analysis,” Chaos: An Interdisciplinary Journal of Non-linear Science , 053103 (2015). T. Coletta, R. Delabays, I. Adagideli, and P. Jacquod, “Topolog-ically protected loop flows in high voltage ac power grids,” NewJournal of Physics , 103042 (2016). S. Park, R. Zhang, J. Lavaei, and R. Baldick, “Monotonicitybetween phase angles and power flow and its implications forthe uniqueness of solutions,” in
Proceedings of the 52nd HawaiiInternational Conference on System Sciences (2019). A. Zachariah, Z. Charles, N. Boston, and B. Lesieutre, “Dis-tributions of the number of solutions to the network power flowequations,” in (IEEE, 2018) pp. 1–5. K. N. Miu and H.-D. Chiang, “Existence, uniqueness, and mono-tonic properties of the feasible power flow solution for radialthree-phase distribution networks,” IEEE Transactions on Cir-cuits and Systems I: Fundamental Theory and Applications ,1502–1514 (2000). J. Lavaei, D. Tse, and B. Zhang, “Geometry of power flows intree networks,” in (IEEE, 2012) pp. 1–8. S. Bolognani and S. Zampieri, “On the existence and linear ap-proximation of the power flow solution in power distributionnetworks,” IEEE Transactions on Power Systems , 163–172(2015). W. A. Bukhsh, A. Grothey, K. I. McKinnon, and P. A. Trod-den, “Local solutions of the optimal power flow problem,” IEEETransactions on Power Systems , 4780–4788 (2013). B. Cui and X. A. Sun, “Solvability of power flow equationsthrough existence and uniqueness of complex fixed point,” arXivpreprint arXiv:1904.08855 (2019). J. W. Simpson-Porco, F. Dörfler, and F. Bullo, “Voltage col-lapse in complex power grids,” Nature communications , 10790(2016). A. J. Wood, B. F. Wollenberg, and G. B. Sheblé,
Power gener-ation, operation, and control (John Wiley & Sons, 2013). T. Nishikawa and A. E. Motter, “Comparative analysis of existingmodels for power-grid synchronization,” New Journal of Physics , 015012 (2015). H. Daido, “Quasientrainment and slow relaxation in a populationof oscillators with random and frustrated interactions,” Physicalreview letters , 1073 (1992). D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley,“Solvable model for chimera states of coupled oscillators,” Phys-ical review letters , 084103 (2008). D. Witthaut and M. Timme, “Kuramoto dynamics in hamilto-nian systems,” Physical Review E , 032917 (2014). S. H. Strogatz,
Nonlinear Dynamics and Chaos with Student So-lutions Manual: With Applications to Physics, Biology, Chem-istry, and Engineering (CRC Press, 2018). W. Chen, J. Liu, Y. Chen, S. Z. Khong, D. Wang, T. Başar,L. Qiu, and K. H. Johansson, “Characterizing the positivesemidefiniteness of signed laplacians via effective resistances,” in (IEEE, 2016) pp. 985–990. J. Schiffer, D. Goldin, J. Raisch, and T. Sezi, “Synchronization of droop-controlled microgrids with distributed rotational andelectronic generation,” in (2013) pp. 2334–2339. M. E. J. Newman,
Networks: an introduction (Oxford UniversityPress, 2010). C. Godsil and G. F. Royle,
Algebraic graph theory , Vol. 207(Springer Science & Business Media, 2013). H. Ronellenfitsch, M. Timme, and D. Witthaut, “A dual methodfor computing power transfer distribution factors,” IEEE Trans-actions on Power Systems , 1007–1015 (2016). H. Ronellenfitsch, D. Manik, J. Hörsch, T. Brown, and D. Wit-thaut, “Dual theory of transmission line outages,” IEEE Trans-actions on Power Systems , 4060–4068 (2017). J. Hörsch, H. Ronellenfitsch, D. Witthaut, and T. Brown, “Lin-ear optimal power flow using cycle flows,” Electric Power SystemsResearch , 126–135 (2018). R. Diestel,
Graph Theory (Springer, New York, 2010). Tong Wu, Z. Alaywan, and A. D. Papalexopoulos, “Locationalmarginal price calculations using the distributed-slack power-flowformulation,” IEEE Transactions on Power Systems , 1188–1190 (2005).51