arXiv:1412.8054v2 [math.OC] 15 Nov 2016

Power Flow as an Algebraic System
Jakub Marecek, Timothy McCoy, and Martin Mevissen

∗

January 24, 2018

Abstract
Steady states of alternating-current (AC) circuits have been studied
in considerable detail. In 1982, Baillieul and Byrnes derived an upper
bound on the number of steady states in a loss-less AC circuit [IEEE
TCAS, 29(11): 724–737] and conjectured that this bound holds for AC
circuits in general. We prove this is indeed the case, among other results,
by studying a certain multi-homogeneous structure in an algebraisation.

1

Introduction

For more than 60 years [37], steady states of alternating-current (AC) circuits
have been studied in considerable detail. The key problem, sometimes known
as the power flow or load flow problem, considers complex voltages Vk at all
buses k as variables, except for one reference bus (k = 0), where the power
supplied. When one denotes complex admittance matrix Y , complex current
Ik , and complex power Sk at bus k, the steady-state equations are based on:
X
X
∗
∗
Sk = Vk Ik∗ = Vk
Yk,l
Vl∗ =
Yk,l
Vk Vl∗
(1)
l∈N

l∈N

where asterisk denotes complex conjugate. This captures the complex, nonconvex non-linear nature [16, 38, 11, 15, 18] of any problem in the AC model.
In order to obtain an algebraic system from (1), one needs to reformulate the
complex conjugate. In order to do so, one may replace all Vk∗ with independent
variables Uk , and filter for “real” solutions where Uk = Vk∗ = ℜVk − ℑVk ı once
the complex solutions are obtained. Thereby, we obtain a particular structure,
which allows us to prove a variety of results.
In particular, the main contributions of our paper are:
• a reformulation of the steady-state equations to a multi-homogeneous algebraic system
∗ J. Marecek and M. Mevissen are with IBM Research – Ireland, IBM Technology Campus
Damastown, Dublin D15, Ireland. e-mail: jakub.marecek@ie.ibm.com. T. McCoy is with
Google.

1

• analytical results on the number and structure of feasible solutions considering losses, resolving a conjecture of Baillieul and Byrnes [2], which
has been open for over three decades
• empirical results for some well-known instances, including the numbers of
roots, conditions for non-uniqueness of optima, and tree-width.
Our analytical results rely on the work of Morgan and Sommese [32, 36] on multihomogeneous structures. Our empirical results rely on Bertini [6], a leading
implementation of homotopy-continuation methods. As we explain in Section
8, ours is not the first algebraisation of the system (1), cf. [37, 4, 3, 2], and
there is a long history [34, 14, 26, 23, 31, 29, 9, 39, 41] of the use of homotopycontinuation methods.

2

The Problem

In order to make the paper self-contained, we restate of the steady-state equations. Consider a circuit represented by an undirected graph (N, E), where
vertices n ∈ N are called buses and edges {l, m} ∈ E ⊆ N × N are called
branches, and an admittance matrix Y = G + Bı ∈ C|N |×|N |, where the real
part of an element is called conductance G = (glm ) and the imaginary part
susceptance B = (blm ). Each bus k ∈ N is associated with complex voltage
Vk = ℜVk + ℑVk ı, complex current Ik = ℜIk + ℑIk ı, and power Sk = Pk + Qk ı
demanded or generated. Let 0 ∈ N correspond to a reference bus with phase
ℑV0 = 0 and magnitude |V0 | fixed; powers at all other buses are fixed too. (In
a variety of extensions, there are other buses, denoted generators, where voltage magnitude, but not the phase and not the power is fixed.) Each branch
(l, m) ∈ E is associated with the complex power Slm = Plm + Qlm ı. The
key constraint linking the buses is Kirchhoff’s current law, which stipulates the
sum of the currents injected and withdrawn at each bus is 0. Considering the

2

relationship I = Y V , the steady state equations hence are:
Pkg = Pkd + ℜVk

n
X

(ℜyik ℜVi − ℑyik ℑVi )

i=1

+ ℑVk

n
X

(ℑyik ℜVi + ℜyik ℑVi )

(2)

i=1

Qgk = Qdk + ℜVk

n
X

(−ℑyik ℜVi − ℜyik ℑVi )

i=1

+ ℑVk

n
X

(ℜyik ℜVi − ℑyik ℑVi )

(3)

i=1

Plm = blm (ℜVl ℑVm − ℜVm ℑVl )

(4)

2

2

+ glm (ℜVl + ℑVm − ℑVl ℑVm − ℜVl ℜVm )
Qlm = blm (ℜVl ℑVm − ℑVl ℑVm − ℜVl 2 − ℑVl 2 )
+ glm (ℜVl ℑVm − ℜVm ℑVl − ℜVm ℑVl )
−

b̄lm
(ℜVl 2 + ℑVl 2 )
2

(5)

Additionally, one can optimise a variety of objectives over the steady states.
In one commonly used objective function, one approximates the costs of real
power P0 generated at the reference bus 0 by a quadratic function f0 :
cost := f0 (P0 ).

(6)

(In a variety of extensions, in which there are other buses where the power
is not fixed, there would be a quadratic function for each such bus and the
quadratic function of power would be summed across all of these buses.) In
the Lp -norm loss objective, one computes a norm of the vector D obtained by
summing apparent powers S(u, v) + S(v, u)∀(u, v) ∈ E for:
||D||p =

 X

|S(u, v) + S(v, u)|

(u,v)∈E

p

1/p

.

(7)

The usual ||D||1 is denoted loss below. We consider these objectives only in
Section VII, while our results in Sections IV–VI apply independently of the use
of any objective function whatsoever.

3

Definitions from Algebraic Geometry

In order to state our results, we need some definitions from algebraic geometry.
While we refer the reader to [2, 27] for the basics, we present the concepts
introduced in the past three decades, not yet widely covered by textbooks. For
a more comprehensive treatment, please see [32, 36, 35].
3

Let n ≥ 0 be an integer and f (z) be a system of n polynomial equations in
z ∈ C n with support (A1 , . . . , An ):

P
α1 α2
αn

α∈A1 f1α z1 z2 · · · zn
 f1 (z) =
..
(8)
. P


α1 α2
αn
fn (z) =
α∈An fnα z1 z2 · · · zn ,

where coefficients fiα are non-zero complex numbers. It is well known [12] that
the polynomials define n projective hypersurfaces in a projective space CPn .
Bézout theorem states that either hypersurfaces intersect in an infinite set
with some component of positive dimension, or the number of intersection
points, counted with multiplicity, is equal to the product d1 · · · dn , where di
is the degree of polynomial i. We call the product d1 · · · dn the usual Bézout
number.
The usual Bézout number can be improved by considering:
Definition 1 (Structure). Any partition of the index set {1, . . . , n} into k sets
I1 , . . . , Ik defines a structure. There, Zj = {zi : i ∈ Ij } is known as the group
of variables for each set Ij . The associated degree dij of a polynomial fi with
respect to group Zj is
X
def
αl .
(9)
dij = max
α∈Ai

l∈Ij

We say that fi has multi-degree (di1 , . . . , din ).
Whenever for some j, for all i, the same dij is attained for all α ∈ Ai , we
call the system (11) homogeneous in the group of variables Zj . The projective
space associated to the group of variables Zj in a structure has dimension

|Ij | − 1 if (11) is homogeneous in Zj , and
def
aj =
(10)
|Ij |
otherwise.
Pk
Definition 2 (Multi-homogeneous Bézout Number). Assuming n = j=1 aj ,
the multi-homogeneous Bézout number Béz(A1 , . . . , An ; I1 , . . . , Ik ) is defined as
Qk
a
the coefficient of the term j=1 ζj j , where aj is the associated dimension (10),
Q n Pk
within the polynomial i=1 j=1 dij ζj , in variables ζj , j = 1 . . . k with coefficients dij are the associated degrees (9); that is (d11 ζ1 + d12 ζ2 + . . . + d1k ζk )
(d21 ζ1 + d22 ζ2 + . . . + d2k ζk ) · · · (d2n ζ1 + d2n ζ2 + . . . + dnk ζk ).
Consider the example of Wampler

 p1 (z) =
p2 (z) =

p3 (z) =

[42] in x ∈ C 3 :
x21 + x2 + 1,
x1 x3 + x2 + 2,
x2 x3 + x3 + 3,

(11)

with the usual Bézout number of 8. Considering the partition {x1 , x2 }, {x3 },
where d11 = 2, d12 = 0, d21 = d22 = d31 = d32 = 1. the monomial ζ12 ζ21 is
4

to be looked up in the polynomial 2ζ1 (ζ1 + ζ2 )2 . The corresponding multihomogeneous Bézout number is hence 4 and this is the minimum across all
possible structures.
In general, the multi-homogeneous Bézout number Béz(A1 , . . . , An ; I1 , . . . , Ik )
is an upper bound on the number of isolated roots of (11) in CPa1 × · · · × CPak ,
and thereby an upper bounds the number of isolated finite complex roots of (11).
There are a variety of additional methods for computing the multi-homogeneous
Bézout number, e.g., [42]. In the particular case where A = A1 = · · · = An , we
denote
 Y

k
n
def
a
(12)
dj j ,
Béz(A1 , . . . , An ; I1 , . . . , Ik ) =
a1 a2 · · · ak
j=1

where dj = dij (equal for each i) and the multinomial coefficient


n!
n
def
=
a1 a2 · · · ak
a1 ! a2 ! · · · ak !
Qk
Pk
is the coefficient of j=1 ζjak in (ζ1 + · · · + ζk )n with n = j=1 aj , as above.
In summary, the multi-homogeneous Bézout number provides a sharper
bound on the numberQof isolated solutions of a system of equations than the
n
usual Bézout number i=1 di = d1 · · · dn . In the famous example of the eigenvalue problem [43], it is known that the Bézout number is 2n , whereas there
exists a structure with multi-homogeneous Bézout number of n. We hence
study the multi-homogeneous structure within the steady state equations of
alternating-current circuits.

4

The Multi-Homogeneous Structure

Notice that in order to obtain an algebraic system from the steady-state equations (2–5), one needs to reformulate the complex conjugate. In order to do so,
one may replace all vn∗ with independent variables un , and later filter for solutions where un = vn∗ once the complex solutions are obtained. We denote such
solution “real”. Let G be the set of slack generators for which |vn | is specified,
and assume 0 ∈ G corresponds to a reference node with phase 0. Notice that
the use of variables vn and un produces a multi-homogeneous structure with
variable groups {vn } and {un }:
vn

X

Yn,k uk + un

X

Yn,k uk − un

k

vn

k

X

∗
Yn,k
vk = 2pn

n∈N \G

X

∗
Yn,k
vk = 2qn

n∈N \G

k

k

vn un = |vn |2
v0 = |v0 |, u0 = |v0 |

n ∈ G − {0}
(13)
5

For example for the two-bus network, we obtain:
∗
∗
v1 (Y1,0 u0 + Y1,1 u1 ) + u1 (Y1,0
v0 + Y1,1
v1 ) = 2p1
∗
∗
v1 (Y1,0 u0 + Y1,1 u1 ) − u1 (Y1,0
v0 + Y1,1
v1 ) = 2q1

v0 = |v0 |, u0 = |v0 |

(14)

Using the algebraic system, one can formulate a number of structural results
concerning power flows.

5

An Analysis for s = 1

For the particular multi-homogeneous structure, which is the partition of the
variables into several groups in (13), we can bound the number of isolated
solutions:
Theorem 1. With exceptions on a parameter set of measure zero, the alternatingcurrent power flow (13) has a finite number of complex solutions, which is
bounded above by:


2n − 2
(15)
n−1
Proof. Each equation in the system (13) is linear in the V variables and also
in the V ∗ variables, giving rise to a natural multi-homogeneous structure of
mult-idegree (1, 1). Since the slack bus voltage is fixed at a reference value, the
system has 2n − 2 such equations in (n − 1, n − 1) variables. By the multihomogeneous form of Bézout’s Theorem (see e.g. Theorem 8.4.7 in [36]), the
total number of solutions in multi-projective space CPn−1 × CPn−1 is precisely
the stated bound, counting multiplicity. Some subset of these lie on the affine
patch C n−1 × C n−1 ⊂ CPn−1 × CPn−1 , giving the result.
Notice that this applies also to some well-known instances of alternatingcurrent optimal power flows (ACOPF). For example, the instances of Lesieutre
et al. [22] and Bukhsh et al. [8] have only a single “slack” bus, whose active
and reactive powers are not fixed, and hence the result applies. Notice that the
exception of measure-zero set is necessary; cf. Example 4.1 of [2].
As we will illustrate in the next section, this bound is tight in some cases.
Deciding whether the bound on the number of roots obtained using a particular
multi-homogeneous structure is tight for a particular instance is, nevertheless,
hard. This can be seen from:
Theorem 2 (Theorem 1 of Malajovich and Meer [27]). There does not exist
a polynomial time algorithm to approximate the minimal multi-homogeneous
Bézout number for a polynomial system (11) up to any fixed factor, unless P =
NP.

6

We can, however, show there exists a certain structure among these solutions:
Corollary 1. If there exists a feasible solution of the alternating-current power
flow, then the solution has even multiplicity greater or equal to 2 or another
solution exists.
Proof. The finite number of solutions to the power flows problem of Theorem 1 is
even. Observe that (U, Û ) is a solution of the system (13) if and only if (Û ∗ , U ∗ )
is a solution. This implies that the non-“real” solutions, that is solutions for
which U ∗ 6= Û , necessarily come in pairs. It follows that “real” feasible solutions
are also even in number, counting multiplicity. The result follows.
Note that a solution having multiplicity greater than 1 is a special case that
is highly unlikely in a real system. Moreover, it is easily detected, since the
Jacobian at a solution is nonsingular if and only if the solution has multiplicity
1.

6

Alternating-Current Optimal Power Flow

One may also make the following observations about the alternating-current
optimal power flows, i.e., the problem of optimising an objective over the steady
state:
Remark 1. For the alternating-current power flows, where powers are fixed at
all but the reference bus, whenever there exists a real feasible solution, except
for a parameter set of measure zero, one can enumerate all feasible solutions in
finite time.
Indeed: By Theorem 1, we know there exist a finite number of isolated
solutions to the system (13). By the homotopy-continuation method of Sommese
et al. [36, 6], we can enumerate the roots with probability 1, which allows us
to pick the global optimum, trivially. Notice that Bertini, the implementation
of the method of Sommese [6], makes it possible to check that all roots are
obtained. Notice that the addition of inequalities can be accommodated by
filtering the real roots.
Nevertheless, this method is not practical, as there may be too many isolated
solutions to enumerate. Generically, this is indeed the case, whenever there
are two or more generators with variable output, i.e., buses, whose active and
reactive power is not fixed:
Corollary 2. In the alternating-current optimal power flow problem, i.e., with
s > 1, where powers are variable outside of the reference bus and there are no additional inequalities, the complex solution set is empty or positive-dimensional,
except for a parameter set of measure zero. When the complex solution set
is positive-dimensional, if a smooth real feasible solution exists, then there are
infinitely many real feasible solutions.

7

Proof. For each slack bus after the first, the system has two variables but only
one equation. With s + 1 slack buses, the rank of the Jacobian and hence the
dimension of the complex solution set will be at least s by Lemma 13.4.1 in [36].
Furthermore, if a real feasible solution U exists and the solution set is smooth
at this point, then the local dimension of the complex and real solution sets are
equal at U . Therefore, since the complex solution set is positive-dimensional,
so is the real set at U , and so infinitely many real feasible solutions exist.
Although there are a variety of methods for studying positive-dimensional
systems, including the enumeration of a point within each connected component
[5, 33] and studying the critical points of the restriction to the variety of the
distance function to such points [1], we suggest the method of moments [20, 13]
may be more suitable for studying the feasible set of alternating-current optimal
power flows. It has been shown recently [13] that it allows for very small errors
on systems in dimension of over five thousand.
As is often the case in engineering applications, one may also be interested
in the distance of a point to the set of feasible solutions. Again, by considering
our algebraisation, one can bound the probability this distance is small using
the theorem of Lotz [25] on the zero-set V of multivariate polynomials. This
could be seen as the converse of the results of [24].

8

Method

No. |N | of buses
3
4
5
6

7

8

9

10

11

12

13

14

9

Bézout’s upper bound
A BKK-based upper bound
Theorem 1

16
8
6

64
40
20

256
192
70

1024
864
252

4096
3712
924

16384
15488
3432

65536
63488
12870

262144
257536
48620

1048576
1038336
184756

4194304
4171776
705432

16777216
16728064
2704156

67108864
67002368
10400600

Generic lower bound

6

20

70

252

924

3432

12870

48620

184756

705432

2704156

10400600

Table 1: The maximum number of steady states in a circuit with a fixed number of buses.

Treewidth tw(P )

No. |X| of solutions

minx∈X (cost(x))

No. of min. wrt. cost

avgx∈X (cost(x))

maxx∈X (cost(x))

minx∈X (loss(x))

No. of min. wrt. loss

avgx∈X (loss(x))

maxx∈X (loss(x))

Bukhsh et al. [8]
Klos and Wojcicka [18]
Lavaei and Low [21]
Bukhsh et al. [8]
McCoy et al.
McCoy et al.
Bukhsh et al. [8]
Grainger and Stevenson [45]
Lesieutre et al. [22]
McCoy et al.
McCoy et al.
McCoy et al.
Bukhsh et al. [8]
McCoy et al.
McCoy et al.
McCoy et al.
Bukhsh et al. [8]
Chow [45]
McCoy et al.
McCoy et al.
McCoy et al.
McCoy et al.
McCoy et al.
McCoy et al.

No. |E| of branches

10

case2w
case3KW
case3LL
case3w
case4
case4ac
case4cyc
case4gs
case5w
case6ac
case6ac2
case6b
case6cyc
case6cyc2
case6cyc3
case7
case8cyc
case9
case9g
caseK4
caseK4sym
caseK6
caseK6b
caseK6sym

Source

No. |N | of buses

Instance

2
3
3
3
4
4
4
4
5
6
6
6
6
6
6
7
8
9
9
4
4
6
6
6

1
3
3
2
3
4
4
4
6
6
6
6
6
6
6
7
8
9
9
6
6
15
15
15

1
2
2
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
3
3
5
5
5

2
6
2
2
4
6
4
10
2
23
22
30
30
12
12
2
60
16
2
8
6
36
40
48

8.42
-0.0
1502.07
5.88
1502.51
2005.25
3001.61
2316.37
2003.57
4005.44
4005.32
3005.7
3005.7
1506.87
4502.42
1667.31
6506.11
3504.44
2604.41
2008.62
2009.28
2018.2
6002.79
6003.01

1
1
1
1
1
2
1
1
1
1
1
5
4
1
1
1
1
1
1
1
2
1
1
1

9.04
1250.0
1502.19
5.94
1505.19
3337.5
3003.56
3681.07
2004.6
5574.23
5554.03
3606.3
3606.3
3131.55
4505.95
1668.36
8557.72
5692.02
4103.36
3006.32
3339.29
5075.47
6017.01
6017.75

9.66
1500.0
1502.31
6.01
1507.88
4005.51
3005.52
4595.8
2005.63
6013.26
6012.82
4508.28
4508.28
4508.24
4508.01
1669.41
9511.04
6505.02
5602.31
4005.51
4005.91
6030.0
6019.41
6019.86

0.71
-0.0
0.22
0.55
0.11
0.01
0.01
0.37
0.47
0.02
0.02
0.01
0.01
0.03
0.02
0.08
0.01
0.03
0.08
0.01
0.01
0.02
0.01
0.01

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

1.02
-0.0
0.34
0.58
2.79
2.44
1.96
27.24
1.5
6.26
5.97
3.9
3.9
4.15
3.55
0.26
4.52
1.37
0.26
4.6
3.96
17.16
14.23
14.75

1.33
0.0
0.46
0.61
5.48
3.78
3.92
39.72
2.53
10.51
10.21
5.88
5.88
6.56
5.61
0.45
7.84
1.87
0.45
7.03
7.28
27.25
16.62
16.86

Table 2: Properties of the instances tested.

variable_group V0, V1; variable_group U0, U1;
I0 = V0*Yv0_0 + V1*Yv0_1; I1 = V0*Yv1_0 + V1*Yv1_1;
J0 = U0*Yu0_0 + U1*Yu0_1; J1 = U0*Yu1_0 + U1*Yu1_1;
fv0 = V0 - 1.0; fv1 = I1*U1 + J1*V1 + 7.0;
fu0 = U0 - 1.0; fu1 = -I1*U1 + J1*V1 - 7.0*I;
Figure 1: A Bertini encoding of ACPF on the two-bus instance of Bukhsh et al.
2012, where the impedance of a single branch is 0.04 + 0.2i.

7

Computational Illustrations

In order to illustrate Theorem 1, we first present the maximum number of
steady states in a circuit with a fixed number of buses in Table , and compare
it to the values of our upper bound, Bézout-based upper bound, and the BKKbased upper bound [10, 31]. Notice that the maximum number of steady states
in a circuit with a fixed number of buses is achieved when (N, E) is a clique.
Notice further that the generic lower bound, obtained as the number of solutions
found by tracing the paths, matches the upper bound of Theorem 1 throughout
Table 7.
To illustrate Proposition 1, we have enumerated the steady states using
Bertini, a versatile package for homotopy-continuation methods by Sommese
et al. [36]. See Figure 1 for an example of Bertini input corresponding to
the example above (14), with constants Yu i j representing Yi,j and Yv i j
∗
representing Yi,j
. The results are summarised in Table .
To illustrate Theorem 1 further, we present the values of our upper bound on
a collection of instances widely known in the power systems community. The
instances are mostly available from the Test Case Archive of Optimal Power
Flow (OPF) Problems with Local Optima of Bukhsh1 , while some have appeared
in well-known papers, e.g. [18], and some are available in recent distributions of
Matpower [45], a well-known benchmark. In particular, we present the numbers
of distinct roots of the instances. In all cases, where Theorem 1 applies, the
number of solutions found by tracing the paths matches the upper bound of
Theorem 1, certifying the completeness. In other cases, one could rely on Bertini
certificates of completeness of the search.
Empirically, we observe there exists a unique global optimum in all these instances tested with respect to the L1 -loss objective. For the generation cost objective, however, there are a number of instances (case4ac, caseK4sym, case6b,
case6cyc), where the global optima are not unique. The case of caseK4sym is a
particularly good illustration, where the symmetry between two generators and
two demand nodes in a complete graph results in multiple global optima.
In order to provide material for further study of structural properties, we
present tree-width of the instances in Table 7 in column tw(N ). Notice that
for the well-known small instances, tree-width is 1 or 2, e.g., 1 for the instance
1 http://www.maths.ed.ac.uk/OptEnergy/LocalOpt/ ,

11

accessed November 30th, 2014.

in Figure 1, and 2 for the instance of Lesieutre et al. As the instances grow,
however, this need not be the case: Kloks [17] shows treewidth is not bounded
even in sparse random graphs, with high probability. In complete graphs, such
as caseK4sym with tree-width 3 above, tree-width grows linearly in the number
of vertices.

8

Related Work

There is a long history of study of the number and structure of solutions of
power flows [37, 4, 3, 2]. [37] considered the Bézout bound. [10, 31] considered
a bound based on the work of Bernstein [7] and Kushnirenko [19]. [2, 4] derived
the same expression as in Theorem 1 using intersection theory, but in the lossless
AC model. They highlight that the number of solutions in an alternatingcurrent model with losses is an important open problem. We note that Theorem
1 subsumes Theorem 4.1 of [2] as a special case. Finally, we note that [18]
present a lower bound without a proof and recent papers [39, 44] bound certain
distinguished solutions, but not all solutions; cf. [28].
There is also a long history of applications of homotopy-continuation methods in power systems [34, 14, 26, 23], although often, e.g. in [26], the set-up
of the homotopy restricted the methods to a heuristic, which could not enumerate all the solutions of the power flow [30]. Recently, these have attracted
much interest [31, 29, 9] following the work of Trias [39, 40, 41]. See [28] for an
overview.

9

Conclusions

We hope that the structural results provided will aid the development of faster
solvers for the related non-linear problems [15]. Arguably, one could:
• By using Theorem 1 in the construction of start systems for homotopycontinuation methods [43], allow for larger zero-dimensional systems to be
studied.
• Extend Corollary 2 to finding at least one point in each connected component [33, 1].
• Extend the homotopy-continuation methods to consider inequalities within
the tracing, rather than only in the filtering phase, which could improve
their computational performance considerably.
• Develop methods for the optimal power flow problem, whose complexity
would be superpolynomial only in the tree-width and the number of buses.
The latter two may be some of the most important challenges within the analysis
of circuits and systems.

12

Acknowledgements Parts of this work have been done while Tim was visiting IBM Research. Jakub would like to thank Isaac Newton Institute for
Mathematical Sciences at the University of Cambridge for their generous support for his visits. Dhagash Mehta has kindly provided a variety of suggestions
as the related work.

References
[1] P. Aubry, F. Rouillier, and M. S. E. Din. Real solving for positive dimensional systems. Journal of Symbolic Computation, 34(6):543 – 560, 2002.
[2] J. Baillieul and C. Byrnes. Geometric critical point analysis of lossless power system models. IEEE Transactions on Circuits and Systems,
29(11):724–737, Nov 1982.
[3] J. Baillieul and C. I. Byrnes. Remarks on the number of solutions to the
load flow equations for a power system with electrical losses. In Decision
and Control, 1982 21st IEEE Conference on, pages 919–924, Dec 1982.
[4] J. Baillieul, C. I. Byrnes, and R. B. Washburn. An algebraic-geometric and
topological analysis of the solution to the load-flow equations for a power
system. In Decision and Control including the Symposium on Adaptive
Processes, 1981 20th IEEE Conference on, pages 1312–1320, Dec 1981.
[5] S. Basu, R. Pollack, and M.-F. Roy. A new algorithm to find a point in
every cell defined by a family of polynomials. In Quantifier elimination and
cylindrical algebraic decomposition, pages 341–350. Springer Vienna, 1998.
[6] D. Bates, J. Hauenstein, A. Sommese, and C. Wampler. Numerically Solving Polynomial Systems with Bertini:. Software, Environments, and Tools.
Society for Industrial and Applied Mathematics, 2013.
[7] D. N. Bernstein. The number of roots of a system of equations. Funkcional.
Anal. i Prilozen., 9(3):1–4, 1975.
[8] W. Bukhsh, A. Grothey, K. McKinnon, and P. Trodden. Local solutions
of the optimal power flow problem. IEEE Trans. Power Syst., 28(4):4780–
4788, 2013.
[9] S. Chandra, D. Mehta, and A. Chakrabortty. Equilibria analysis of power
systems using a numerical homotopy method. In Power Energy Society
General Meeting, 2015 IEEE, pages 1–5, July 2015.
[10] T. Chen and D. Mehta. On the Network Topology Dependent Solution
Count of the Algebraic Load Flow Equations. 2015. ArXiv e-prints
1512.04987.

13

[11] H.-D. Chiang, C.-W. Liu, P. P. Varaiya, F. F. Wu, and M. G. Lauby.
Chaos in a simple power system. IEEE Transactions on Power Systems,
8(4):1407–1417, Nov 1993.
[12] W. Fulton. Intersection Theory. Springer New York, Jan. 1998.
[13] B. Ghaddar, J. Marecek, and M. Mevissen. Optimal power flow as a
polynomial optimization problem. IEEE Transactions on Power Systems,
31(1):539–546, Jan 2016.
[14] S. X. Guo and F. M. A. Salam. The real homotopy-based method for
computing solutions of electric power systems. In Circuits and Systems,
1992. ISCAS ’92. Proceedings., 1992 IEEE International Symposium on,
volume 6, pages 2737–2740 vol.6, May 1992.
[15] I. A. Hiskens. Analysis tools for power systems-contending with nonlinearities. Proceedings of the IEEE, 83(11):1573–1587, Nov 1995.
[16] I. A. Hiskens and R. J. Davy. Exploring the power flow solution space
boundary. IEEE Transactions on Power Systems, 16(3):389–395, Aug 2001.
[17] T. Kloks. Only few graphs have bounded treewidth. In T. Kloks, editor,
Treewidth, volume 842 of Lecture Notes in Computer Science, pages 51–60.
Springer Berlin Heidelberg, 1994.
[18] A. Klos and J. Wojcicka. Physical aspects of the nonuniqueness of load flow
solutions. International Journal of Electrical Power and Energy Systems,
13(5):268–276, 1991.
[19] A. Kushnirenko. Newton polytopes and the Bézout theorem. Funct. Anal.
Appl., 10:233–235, 1976.
[20] J. B. Lasserre. Convergent sdp-relaxations in polynomial optimization with
sparsity. SIAM Journal on Optimization, 17(3):822–843, 2006.
[21] J. Lavaei and S. Low. Zero duality gap in optimal power flow problem.
Power Systems, IEEE Transactions on, 27(1):92 –107, feb. 2012.
[22] B. Lesieutre, D. Molzahn, A. Borden, and C. DeMarco. Examining the
limits of the application of semidefinite programming to power flow problems. In Communication, Control, and Computing (Allerton), 2011 49th
Annual Allerton Conference on, pages 1492–1499, 2011.
[23] C.-W. Liu, C.-S. Chang, J. A. Jiang, and G. H. Yeh. Toward a cpflow-based
algorithm to compute all the type-1 load-flow solutions in electric power
systems. IEEE Transactions on Circuits and Systems I: Regular Papers,
52(3):625–630, March 2005.

14

[24] C.-W. Liu and J. S. Thorp. A novel method to compute the closest unstable
equilibrium point for transient stability region estimate in power systems.
IEEE Transactions on Circuits and Systems I: Fundamental Theory and
Applications, 44(7):630–635, Jul 1997.
[25] M. Lotz. On the volume of tubular neighborhoods of real algebraic varieties.
Proc. Amer. Math. Soc., 143:1875–1889, 2015.
[26] W. Ma and J. S. Thorp. An efficient algorithm to locate all the load
flow solutions. IEEE Transactions on Power Systems, 8(3):1077–1083, Aug
1993.
[27] G. Malajovich and K. Meer. Computing minimal multi-homogeneous
Bézout numbers is hard. Theory of Computing Systems, 40(4):553–570,
2007.
[28] D. Mehta, D. K. Molzahn, and K. Turitsyn. Recent advances in computational methods for the power flow equations. In 2016 American Control
Conference (ACC), pages 1753–1765, July 2016.
[29] D. Mehta, H. D. Nguyen, and K. Turitsyn. Numerical polynomial homotopy
continuation method to locate all the power flow solutions. IET Generation,
Transmission Distribution, 10(12):2972–2980, 2016.
[30] D. K. Molzahn, B. C. Lesieutre, and H. Chen. Counterexample to a
continuation-based algorithm for finding all power flow solutions. IEEE
Transactions on Power Systems, 28(1):564–565, Feb 2013.
[31] D. K. Molzahn, D. Mehta, and M. Niemerg. Toward topologically based
upper bounds on the number of power flow solutions. In 2016 American
Control Conference (ACC), pages 5927–5932, July 2016.
[32] A. Morgan and A. Sommese. A homotopy for solving general polynomial
systems that respects m-homogeneous structures. Appl. Math. Comput.,
24(2):101–113, 1987.
[33] F. Rouillier, M.-F. Roy, and M. S. E. Din. Finding at least one point in each
connected component of a real algebraic set defined by a single equation.
Journal of Complexity, 16(4):716–750, 2000.
[34] F. M. A. Salam, L. Ni, S. Guo, and X. Sun. Parallel processing for the
load flow of power systems: the approach and applications. In Decision
and Control, 1989., Proceedings of the 28th IEEE Conference on, pages
2173–2178 vol.3, Dec 1989.
[35] I. R. Shafarevich. Basic algebraic geometry. Springer-Verlag, springer study
edition edition, 1977. Translated from the Russian by K. A. Hirsch; Revised printing of Grundlehren der mathematischen Wissenschaften, Vol.
213, 1974.

15

[36] A. J. Sommese and C. W. Wampler. The numerical solution of systems of
polynomials arising in Engineering and Science. World Scientific, 2005.
[37] C. J. Tavora and O. J. M. Smith. Equilibrium analysis of power systems.
IEEE Transactions on Power Apparatus and Systems, PAS-91(3):1131–
1137, May 1972.
[38] J. S. Thorp and S. A. Naqavi. Load flow fractals. In Decision and Control,
1989., Proceedings of the 28th IEEE Conference on, pages 1822–1827 vol.2,
Dec 1989.
[39] A. Trias. The holomorphic embedding load flow method. In Power and
Energy Society General Meeting, 2012 IEEE, pages 1–8, July 2012.
[40] A. Trias. System and method for monitoring and managing electrical power
transmission and distribution networks, US Patent 7519506 and 7979239,
2009 and 2010, number = US 7519506 and 7979239, type = Patent, version
= , location = US,.
[41] A. Trias and J. L. Marı́n. The holomorphic embedding loadflow method for
dc power systems and nonlinear dc circuits. IEEE Transactions on Circuits
and Systems I: Regular Papers, 63(2):322–333, Feb 2016.
[42] C. W. Wampler. Bezout number calculations for multi-homogeneous polynomial systems. Applied Mathematics and Computation, 51(2):143 – 157,
1992.
[43] C. W. Wampler. An efficient start system for multi-homogeneous polynomial continuation. Numerische Mathematik, 66(4):517–524, 1993/94.
[44] T. Wang and H. D. Chiang. On the number of system separations in
electric power systems. IEEE Transactions on Circuits and Systems I:
Regular Papers, 63(5):661–670, May 2016.
[45] R. Zimmerman, C. Murillo-Sánchez, and R. Thomas. Matpower: Steadystate operations, planning, and analysis tools for power systems research
and education. Power Systems, IEEE Transactions on, 26(1):12–19, 2011.

16

