Preprint No: ADP-15-34/T936

1

Toward Topologically Based Upper Bounds on the
Number of Power Flow Solutions

arXiv:1509.09227v1 [math.OC] 30 Sep 2015

Daniel K. Molzahn∗ , Dhagash Mehta† , and Matthew Niemerg§

Abstract—The power flow equations, which relate power injections and voltage phasors, are at the heart of many electric power
system computations. While Newton-based methods typically
find the “high-voltage” solution to the power flow equations,
which is of primary interest, there are potentially many “lowvoltage” solutions that are useful for certain analyses. This paper
addresses the number of solutions to the power flow equations.
There exist upper bounds on the number of power flow solutions;
however, there is only limited work regarding bounds that are
functions of network topology. This paper empirically explores
the relationship between the network topology, as characterized
by the maximal cliques, and the number of power flow solutions. To facilitate this analysis, we use a numerical polynomial
homotopy continuation approach that is guaranteed to find all
complex solutions to the power flow equations. The number of
solutions obtained from this approach upper bounds the number
of real solutions. Testing with many small networks informs the
development of upper bounds that are functions of the network
topology. Initial results include empirically derived expressions
for the maximum number of solutions for certain classes of
network topologies.

I. I NTRODUCTION
The power flow equations model the steady-state relationship between the voltages and power injections in a power
system. Solutions to these nonlinear equations correspond to
the equilibria of the underlying differential-algebraic equations
describing power system dynamic behavior. Multiple power
flow solutions may exist [1]. Typically, lightly loaded systems
have many solutions which disappear in bifurcations with
increasing loading. After reaching the system’s maximum
loadability limit, the power flow equations are infeasible. This
paper presents a empirical analysis leading towards development of topologically based upper bounds on the number of
power flow solutions.
For typical operating conditions, there is a single “highvoltage” solution corresponding to the stable equilibrium of
reasonable dynamical models. The voltage magnitudes associated with this solution are usually close to nominal values,
giving rise to the name “high-voltage” solution.
There exist mature techniques which often succeed in calculating the high-voltage solution. Typical power flow solvers
employ Newton-Raphson and Gauss-Seidel algorithms [2].
∗

Argonne

National

Laboratory,

Energy

Systems

Division.

dmolzahn@anl.gov
† University of Notre Dame, Dept. of Applied and Computational
Mathematics and Statistics; University of Adelaide, School of Physical
Sciences, Dept. of Physics, Centre for the Subatomic Structure of Matter. Support from NSF-ECCS award ID 1509036 and an Australian Research Council DECRA fellowship no. DE140100867. dmehta@nd.edu,
Preprint No: ADP-15-34/T936

§

Fields

Institute

for

Research

research@matthewniemerg.com

in

Mathematical

Sciences:

For a sufficiently close initialization, Newton-based techniques yield quadratic convergence to a power flow solution.
However, the capabilities of Newton-based techniques are
very dependent on the initialization. Reasonable initializations
(e.g., the solution to a related problem or a “flat start” of
1∠0◦ per unit voltages) often result in convergence to the
high-voltage solution. Obtaining an appropriate initialization is
more difficult for parameters outside typical operating ranges,
as may occur due to high penetrations of renewable generation and during contingencies. Illustrating this initialization
challenge, [3] demonstrates that the basins of attraction for
Newton-based methods are fractal in nature.
This challenge motivates the development of alternative
solution techniques. Many alternative algorithms are Newtonbased with modifications that improve robustness to the choice
of initialization. For instance, [4] describes an “optimal multiplier” which modifies the step size to ensure non-divergence
of the iterations.
Many recent approaches do not rely on Newton-based
iteration. For instance, the Holomorphic Embedding Load
Flow method [5], which uses analytic continuation theory from
complex analysis, is claimed to be capable of reliably finding
a stable solution for any feasible set of power flow equations.
A monotone operator approach in [6] identifies regions with at
most one power flow solution. Any solution contained within
these regions can be quickly calculated. Convex relaxation
techniques can calculate power flow solutions [7] and certify
infeasibility [8]–[10]. Progress has also been made using
active/reactive power “decoupling” approximations [11], [12].
Many of these approaches focus primarily on calculation of
the high-voltage power flow solution. However, other power
flow solutions generally exist, often at lower voltages. Locating
multiple power flow solutions is important for many types
of power system stability assessments [13]–[16]. There may
also exist multiple stable power flow solutions, particularly in
the presence of power flow reversal conditions on distribution
systems [17]. The literature includes attempts to calculate
multiple power flow solutions using a Newton-based algorithm
with a range of carefully chosen initializations [18], [19]. A
semidefinite relaxation is used in [20] to identify multiple
power flow solutions. However, these approaches may not find
all solutions.
A numerical continuation approach that claims to find all
power flow solutions was presented in [21]. Since it scales with
the number of actual rather than potential power flow solutions, this approach is computationally tractable for realistic
power systems. However, the robustness proof indicating that
this approach finds all solutions to all power flow problems
is flawed [22], and [23] presents a counterexample. This

2

Preprint No: ADP-15-34/T936

counterexample also invalidates the robustness of the related
approach in [24] for calculating all power flow solutions that
have a certain stability property. Recent work [25] presents a
related method that improves the robustness of [21].
Although not yet scalable to large systems, there are approaches which are guaranteed to find all power flow solutions.
For instance, interval analysis [26] and Gröbner bases [27]
methods are applicable to systems with up to five buses.
Other methods guaranteed to find all solutions to systems of
polynomial equations are the eigenvalue technique in [28] and
the “moment/sum-of-squares” relaxations in [29].
The most computationally tractable methods that reliably
find all solutions are based on numerical polynomial homotopy
continuation (NPHC). These methods use continuation to trace
all complex solutions from a selected polynomial system for
which all solutions can be easily characterized to the solutions
of the specified target system. In this context, the power flow
equations are transformed by splitting real and imaginary parts
of the voltage phasors to obtain polynomials in real variables.
Thus, only the real solutions to these polynomial equations
are physically meaningful. Nevertheless, ensuring recovery of
all real solutions requires a number of continuation traces
that depends on an upper bound for the number of complex
solutions. Existing techniques are tractable for systems with
up to 14 buses [30], [31] (and the equivalent of 18 buses for
the related Kuramoto model [32]).
The computational complexity of NPHC is unknown and is
related to Smale’s 17th open problem of the 21st century [33].
As such, the complexity of NPHC is considered in terms of the
maximal number of complex solutions of a polynomial system
in a fixed parametrized family of coefficients. Obtaining tighter
upper bounds on the number of possible power flow solutions
for a specified network thus has the potential to improve the
tractability of NPHC. Upper bounds on the number of power
flow solutions are also theoretically interesting in their own
right.
After introducing the power flow problem, this paper reviews existing bounds and clarifies inaccurate statements in the
literature on the maximum number of power flow solutions.
To improve upon existing bounds, this paper then uses NPHC
to perform an empirical analysis regarding the number of
complex power flow solutions. By testing a variety of network
topologies, the results inform conjectured network-structure
dependent upper bounds.

II. OVERVIEW OF THE P OWER F LOW E QUATIONS
The nonlinear power flow equations model the steady-state
relationship between voltages and power injections in a power
system. Denote the voltage phasors for a n-bus system as V ∈
Cn . We use rectangular coordinates Vi = Vdi + jVqi , where
Vd , Vq ∈ Rn and j is the imaginary unit, to obtain a system
of polynomial equations in real variables Vd and Vq .
Let Y = G + jB represent the network admittance matrix. Using “active/reactive” representation of complex power
injections Pi + jQi , power balance at bus i yields

Pi =Vdi

n
X
k=1

(Gik Vdk − Bik Vqk ) + Vqi

n
X

(Bik Vdk + Gik Vqk )

k=1

(1a)
Qi =Vdi

n
X
k=1

(−Bik Vdk − Gik Vqk )+Vqi

n
X

(Gik Vdk − Bik Vqk )

k=1

(1b)

The squared voltage magnitude at bus i is
|Vi |2 = Vdi2 + Vqi2

(1c)

Each bus is classified as PQ, PV, or slack. PQ buses, which
typically correspond to loads and are denoted by the set PQ,
treat Pi and Qi as specified quantities and enforce the active
power (1a) and reactive power (1b) equations. PV buses, which
typically correspond to generators and are denoted by the set
2
PV, enforce (1a) and (1c) with specified Pi and |Vi | . The
associated reactive power Qi may be computed as an “output
quantity” via (1b).1 Finally, a single slack bus is selected with
specified Vi (typically chosen such that the reference angle
is 0◦ ; i.e., Vqi = 0). The set S denotes the slack bus. The
active power Pi and reactive power Qi at the slack bus are
determined from (1a) and (1b).
III. R EVIEW OF B OUNDS ON THE N UMBER OF P OWER
F LOW S OLUTIONS
Zero active and reactive power injections at a bus correspond to either zero current injection or zero voltage magnitude. Systems composed solely of PQ buses with the exception
of a single slack bus can thus achieve at least 2(n−1) real
solutions when all power injections at PQ buses are zero.
This result has led some publications (e.g., [21], [34], [35])
to claim that there exist at most 2(n−1) real power flow
solutions. However, there are systems with more than 2(n−1)
real solutions: [36]–[38] present three-bus systems with six
solutions, which is larger than 2(n−1) = 4. Thus, while the
approaches in [21], [34], [35] may find many power flow
solutions, they do not find all solutions for all systems.
Bézout’s theorem with 2n − 2 degree-two polynomials
bounds the number of power flow solutions at 2(2n−2) [39].
The number of complex solutions for sparse systems of
polynomial equations with generic coefficients is given by
the BKK bound [40], [41]. For this particular power flow
formulation for systems with PQ buses, the BKK bound is the
same as the total degree bound from Bézout’s theorem (i.e.,
2(2n−2) ). However, each unique network topology restricts
the power flow equations to a parametrized family in the
coefficient space, resulting in a deficiency in the number of
solutions and the potential for tighter bounds. By exploiting
the structure of the power flow equations for systems without
PQ buses, [39] proves a tighter bound of


2n − 2
κn =
(2)
n−1

where ·· is the binomial coefficient function. The number
of complex solutions for general systems (i.e., a mix of PV
and PQ buses) is also bounded by κn [42].
1 This

paper does not consider reactive-power limited generators.

Preprint No: ADP-15-34/T936

Maximum Number of Complex Solutions vs. Average Clique Size
κ

7

900

Maximum Number of Complex Solutions

An upper bound of 2(n−1) + m, where m is an unspecified
network-dependent parameter, is conjectured in [38], but no
guidance is provided as to the value of m.
Under some technical assumptions, [43] shows that the maximum number of certain (“type-1”) unstable powerflow solu
2 |PV|
tions associated with a given stable equilibrium is 2
,
|PV|
where |PV| is the number of PV buses. However, [43] does
not bound the total number of solutions.
The aforementioned bounds are not functions of the network
topology. The sole topology-dependent bound is given in [44]:
for a system comprised of certain subnetworks which share a
single bus, the number of complex solutions is the product
of the number of solutions for each subnetwork. In order to
significantly improve on the bound κn for complete networks,
we speculate that developing topologically based bounds is
possible for a broader class of networks. The remainder of this
paper presents an empirical investigation intended to inform
such topologically based bounds.

3

n=3
n=4
n=5
n=6
n=7

800
700
600
500
400
300

κ6

200
100
0
2

3

4
5
Average Clique Size

6

κ5
κ4
κ3
7

Fig. 1. Maximum number of complex solutions vs. average clique size.
Dotted lines correspond to the upper bound κn for complete networks.

IV. OVERVIEW OF THE N UMERICAL P OLYNOMIAL
H OMOTOPY A LGORITHM
With the NPHC method [41], [45], [46], one first creates a
start system that is in the same parametrized family as the
target system. Then, one easily computes all the solutions
of the start system and deforms these solutions to all the
solutions of the target system. This deformation of solutions is
done by creating a one-parameter embedding between the start
system and the target system known as a homotopy. Under
suitable conditions, each isolated solution of the target system
is connected by a homotopy path to at least one solution
of the start system. By tracking the paths over the complex
numbers, one can be guaranteed to find all of the isolated,
complex solutions of a polynomial system using NPHC. Once
all the complex solutions are found, we then consider only
the solutions with a small enough imaginary part and deem
these to be real. For all of our computations, we use Bertini
v1.5 [47].
V. E MPIRICAL R ESULTS AND D ISCUSSION
To discover potential topology-dependent bounds, we used
NPHC to calculate all complex solutions for many small test
cases. In this section, we describe the method employed for
generating test cases, present the results, and discuss some
patterns discovered thus far in the data.
A total of 50,000 systems (10,000 for each size from
three to seven buses) were randomly constructed. The systems
were created by sampling from a uniform distribution for the
2
number of lines (n + 1 to n 2+n ), with a topology developed
from a random spanning tree [48] augmented with additional
lines whose terminal buses were randomly selected. See the
appendix for further test case construction details.
We characterize the power system network topology using
its maximal cliques. A clique is a subgraph where all nodes
are linked to one another (i.e., a complete subgraph), and a
maximal clique is clique that is not a subgraph of another
clique. The set of maximally sized cliques of a given graph can

be computed using the Bron-Kerbosch algorithm [49]. Despite
being NP-Complete [50], calculating the maximal cliques of
sparse graphs representing power system networks is tractable
(e.g., identifying the maximal cliques of the 2383-bus Polish
system [51] only required 35 seconds on a laptop).
Fig. 1 shows the maximum number of complex solutions
versus the average clique size (i.e., for a system P
with maximal
m
cliques C1 , . . . , Cm , the average clique size is i=1 |Ci | /m,
where |Ci | is the size of clique i) for each of the 144
unique combinations of clique structures (i.e., unique sets
{|C1 | , . . . , |Cm |}) among the 50,000 test cases. Larger average clique size (roughly speaking, a more tightly connected
network) is correlated with more complex solutions, although
the dependence is not monotonic.
The dotted lines in Fig. 1 show the values of κn . The results
for the complete networks show that the bounds κn on the
maximum number of complex solutions are exact for systems
with up to seven buses. For systems with more than three
buses, it is an open question whether there exist power flow
equations with a number of real solutions equal to κn .2
Fig. 1 illustrates that the number of complex power flow
solutions strongly depends on the network structure. The
fact that power systems typically have sparse topologies with
relatively small maximal cliques suggests that bounds based
on the network structure have the potential to be significantly
tighter than bounds developed for complete networks.
We next describe bounds for two classes of networks. The
first class is block networks. The blocks of a graph are the
maximally sized subgraphs such that the subgraph remains
connected after removal of any single node. A network with
the property that all blocks are cliques is a block network [52].
2 Two-bus systems can have two real solutions and the three-bus systems
in [36]–[38] have six real solutions, so this question is answered affirmatively
for κ2 and κ3 . A complete lossless four-bus network of PV buses with
identical line reactances has at most 14 real solutions [39], which is strictly
less than κ4 = 20. A variety of perturbations to the electrical parameters for
this system still resulted in at most 14 real solutions.

4

Preprint No: ADP-15-34/T936

3
2

3

1

C1

4
1

2

C2

4

5
(a) Five-bus system with maximal cliques C1 = {1, 2}, C2 = {1, 3},
C3 = {1, 4}, and C4 = {1, 5}. Maximum number of solutions is κ(1) =
κ2 κ2 κ2 κ2 = 2 · 2 · 2 · 2 = 16.

1

(a) Four-bus system with maximal cliques C1 = {1, 2, 3} and C2 = {2, 3, 4}.
3 κ3
The (conjectured) maximum number of solutions is κ(2) = κ(2−1)
= 6·6
=
2
2
18.

4

5

3

C3
3

2

1

(b) Four-bus system with maximal cliques C1 = {1, 2, 3} and C2 = {3, 4}.
Maximum number of solutions is κ(1) = κ3 κ2 = 2 · 6 = 12.

1
2

5

C1

4

C2
C4

7

4
3

2

6
6

(b) Six-bus system with maximal cliques C1 = {1, 2, 3}, C2 = {2, 3, 4},
C3 = {3, 4, 5}, and C4 = {2, 4, 6}. The (conjectured) maximum number of
3 κ3 κ3
= 162.
= 6·6·6·6
solutions is κ(2) = κ3 κ(4−1)
23
2

(c) Seven-bus system with maximal cliques C1 = {1, 2, 3}, C2 = {3, 4, 5},
and C3 = {5, 6, 7}. Maximum number of solutions is κ(1) = κ3 κ3 κ3 =
6 · 6 · 6 = 216.

1

4

5

2

3

6

Fig. 2. Topologies for example test cases with a block network structure
(i.e., the upper bound κ(1) in (3) is applicable).

C1

First shown in [44], the number of solutions for this class of
networks, denoted κ(1) , is given by the product of the number
of solutions for each clique:
κ(1) =

m
Y

κ|Ck | .

(3)

i=1

Fig. 2 shows example block network topologies with corresponding κ(1) values.
The second class of networks is defined based on the
corresponding clique graph (i.e., the graph where the maximal
cliques are the nodes and the edges represent the buses shared
by the maximal cliques). Specifically, the definition of the
second class of networks is that the clique graph for those
maximal cliques which share exactly two buses (i.e., a single
edge) is a tree.3 The NPHC results empirically suggest (but do
not prove) that the number of solutions for this second class
of networks, denoted κ(2) , is at most
κ(2) =

1
2m−1

m
Y

κ|Ck | .

C2

C3

7

(c) Seven-bus system with maximal cliques C1 = {1, 2, 3, 4}, C2 =
{3, 4, 5, 6}, and C3 = {5, 6, 7}. The (conjectured) maximum number of
κ4 κ3
= 600.
solutions is κ(2) = κ4(3−1)
= 20·20·6
22
2

Fig. 3. Topologies (left) and clique graphs (right) for example test cases
where the maximal cliques which share two buses are arranged in a tree (i.e.,
the conjectured upper bound κ(2) in (4) is applicable).

Note that other factors, such as electrical parameters and
the number and locations of PV buses, may reduce the actual
number of complex solutions for a specific problem below the
expression in (4).
To compute the (conjectured) maximum number of solutions
for systems where some subnetworks have a block topology
and others are arranged such that the maximal cliques have

(4)

3

C1

4

C4

C5

C3

i=1

The left hand sides in Fig. 3 show example topologies of
systems for which κ(2) is applicable. The right hand sides
of Fig. 3 show the associated clique graphs, where the nodes
represent the maximal cliques and each edge represents a bus
shared by the corresponding maximal cliques.
3 As illustrated in Figs. 3 – 5, “tree” in this context refers to the lack of
loops in connections between maximal cliques which share an edge; i.e.,
these maximal cliques, but not necessarily the buses, are connected in a tree.

1

2

7

5

6

C2

Fig. 4. Seven-bus system with maximal cliques C1 = {1, 2}, C2 = {3, 4},
C3 = {5, 6}, C4 = {2, 3, 7}, and C5 = {3, 5, 7}. The (conjectured)
3 κ3
maximum number of solutions is κ(1) κ(2) = κ2 κ2 κ2 κ(2−1)
= 2·2·2· 6·6
=
2
2
144. The loop in the clique graph is associated with the maximal cliques
sharing a single bus; the restriction of the clique graph to a tree only applies
to those cliques sharing two buses.

Preprint No: ADP-15-34/T936

5

two shared buses (e.g., Fig. 4), multiply the values from (3)
and (4) for the appropriate subnetworks.
The expressions κ(1) and κ(2) suggest an extension to
systems with maximal cliques that share
Qm more than two buses:
for an integer α ≥ 3, κ(α) = β1α i=1 κ|Ck | , where βα is
some integer function of α (e.g., βα = αm−1 ). However,
this extension does not hold. See, for instance, the five-bus
system in Fig. 5a, which has two size-four cliques that share
α = 3 buses. Since κ4 = 20, we would expect that the
number of solutions for the five-bus system would have the
form κβ4 κ3 4 = 400
β3 . The empirical results show a maximum
of 68 solutions for this topology, and 400
β3 = 68 implies that
β3 = 100/17 = 5.8823 . . ., which is not an integer.
Figs. 5b and 5c show two other example topologies for
which we have not yet been able to determine a pattern in
the maximum number of complex solutions. Both of these
topologies have a loop in the clique graph among maximal
cliques that share two buses. Thus, κ(2) is not applicable.
VI. C ONCLUSION
The power flow equations are intrinsic to many power
system analyses. Yet, despite their importance, we lack some
basic information such as the number of solutions to these
equations for general systems. One practical implication is

1

3

2

C1
5

C2

4

(a) Five-bus system with maximal cliques C1 = {1, 2, 4, 5} and C2 =
{2, 3, 4, 5}. Maximum of 68 solutions found among all test cases with
isomorphic topologies.

3

C1

1

C2

5

4

C3

2

(b) Five-bus system with maximal cliques C1 = {1, 2, 3}, C2 = {2, 3, 4},
and C3 = {2, 3, 5}. Maximum of 52 solutions found among all test cases
with isomorphic topologies.

1

C1

3
2

C2

4

5

C3
C4

(c) Five-bus system with maximal cliques C1 = {1, 2, 3}, C2 = {1, 3, 4},
C3 = {3, 4, 5}, and C4 = {2, 3, 5}. Maximum of 66 solutions found among
all test cases with isomorphic topologies.
Fig. 5. Example test cases without known expressions for upper bounds on
the number of solutions.

that numerical polynomial homotopy continuation methods
require more computational effort than necessary to find all
power flow solutions. After reviewing upper bounds on the
number of power flow solutions for complete networks, this
paper presented results from a computational experiment using
numerical polynomial homotopy continuation to calculate all
complex power flow solutions for 50,000 small test cases.
These results illustrate the potential for improving existing
bounds by exploiting network topology. The results validated
the only topology-dependent bound in the existing literature,
which is applicable to block networks, and informed the development of conjectured bounds for another class of network
topologies.
Future work includes scaling the computations to larger
systems, uncovering other patterns to develop further empirical
expressions relating the maximum number of power flow
solutions to the network topology, and identifying physical
explanations for these expressions.
A PPENDIX
For each test case, line resistance R, reactance X, and shunt
susceptance b values for Π-model equivalent circuits were
randomly sampled from Gaussian distributions with mean and
standard deviation of µR = 0.03, σR = 0.03; µX = 0.10,
σX = 0.03; and µb = 0.005, σb = 0.001 respectively, all in
per unit using a 100 MVA base, with any negative values
sampled for line resistances instead set to zero. Lines had
an 8% probability of being transformers with tap ratio τ and
phase-shift θ sampled from a Gaussian distribution with mean
and standard deviation values of µτ = 1, στ = 0.02 per unit
and µθ = 0◦ , σθ = 3◦ , respectively. A bus was specified to
be a generator with 30% probability, with the first generator
selected as the slack bus. If no buses were selected to be
generators, a random bus was assigned a generator and chosen
to be the slack bus. The generators’ voltage setpoints were
selected by randomly sampling from a uniform distribution
between 0.90 and 1.10 per unit, and active power injections
were sampled from Gaussian distribution with mean and
standard deviation values of µPg = 200, σPg = 30 MW.
(Generator buses had no associated loads.) Buses had a 70%
probability of being specified as loads (PQ buses). Loads
have a constant active and reactive power component sampled
from a Gaussian distribution with mean and standard deviation
µPd = 50, σPd = 20 MW and µQd = 30, σQd = 20 MVAr,
respectively, as well as a constant impedance component
sampled from a Gaussian distribution with mean and standard
deviation µPz = 0.1, σPz = 0.03 per unit and µQz = 0.1,
σQz = 0.05 per unit for shunt conductance and susceptance,
respectively.
R EFERENCES
[1] A. Klos and A. Kerner, “The Non-Uniqueness of Load Flow Solution,”
in 5th Power Syst. Comput. Conf. (PSCC), 1975.
[2] J. Glover, M. Sarma, and T. Overbye, Power System Analysis and
Design. Cengage Learning, 2011.
[3] J. Thorp and S. Naqavi, “Load Flow Fractals,” in IEEE 28th Ann. Conf.
Decis. Contr. (CDC), Dec. 1989.
[4] S. Iwamoto and Y. Tamura, “A Load Flow Calculation Method for IllConditioned Power Systems,” IEEE Trans. Power App. Syst., vol. PAS100, no. 4, pp. 1736–1743, Apr. 1981.

6

[5] A. Trias, “The Holomorphic Embedding Load Flow Method,” in IEEE
PES General Meeting, Jul. 2012.
[6] K. Dvijotham, S. Low, and M. Chertkov, “Solving the Power Flow
Equations: A Monotone Operator Approach,” arXiv:1506.08472, 2015.
[7] R. Madani, J. Lavaei, and R. Baldick, “Convexification of Power Flow
Problem over Arbitrary Networks,” IEEE 54th Ann. Conf. Decis. Contr.
(CDC), Dec. 2015.
[8] D. Molzahn, B. Lesieutre, and C. DeMarco, “A Sufficient Condition
for Power Flow Insolvability With Applications to Voltage Stability
Margins,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2592–2601, Aug.
2013.
[9] D. Molzahn, V. Dawar, B. Lesieutre, and C. DeMarco, “Sufficient
Conditions for Power Flow Insolvability Considering Reactive Power
Limited Generators with Applications to Voltage Stability Margins,” in
IREP Symp., Aug. 25-30 2013.
[10] D. Molzahn, I. Hiskens, and B. Lesieutre, “Calculation of Voltage
Stability Margins and Certification of Power Flow Insolvability using
Second-Order Cone Programming,” in 47th Hawaii Int. Conf. Syst. Sci.
(HICSS), Jan. 5-8 2016.
[11] M. Ilic, “Network Theoretic Conditions for Existence and Uniqueness
of Steady State Solutions to Electric Power Circuits,” in IEEE Int. Symp.
Circuits Syst. (ISCAS), May 1992.
[12] F. Dörfler and F. Bullo, “Novel Insights into Lossless AC and DC Power
Flow,” in IEEE PES General Meeting, Jul. 2013.
[13] V. Venikov, V. Stroev, V. Idelchick, and V. Tarasov, “Estimation of Electrical Power System Steady-State Stability in Load Flow Calculations,”
IEEE Trans. Power App. Syst., vol. 94, no. 3, pp. 1034–1041, May 1975.
[14] 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., vol. PAS-102, no. 5, pp. 1115–1125, May
1983.
[15] M. Ribbens-Pavella and F. Evans, “Direct Methods for Studying Dynamics of Large-Scale Electric Power Systems-A Survey,” Automatica,
vol. 21, no. 1, pp. 1–21, 1985.
[16] H.-D. Chiang, F. Wu, and P. Varaiya, “Foundations of Direct Methods
for Power System Transient Stability Analysis,” IEEE Trans. Circuits
Syst., vol. 34, no. 2, pp. 160–173, Feb. 1987.
[17] H. Nguyen and K. Turitsyn, “Appearance of Multiple Stable Load Flow
Solutions under Power Flow Reversal Conditions,” in IEEE PES General
Meeting, Jul. 2014.
[18] T. Overbye and R. Klump, “Effective Calculation of Power System LowVoltage Solutions,” IEEE Trans. Power Syst., vol. 11, no. 1, pp. 75–82,
Feb. 1996.
[19] R. Klump and T. Overbye, “A New Method for Finding Low-Voltage
Power Flow Solutions,” in IEEE PES Summer Meeting, 2000.
[20] B. Lesieutre, D. Molzahn, A. Borden, and C. DeMarco, “Examining the
Limits of the Application of Semidefinite Programming to Power Flow
Problems,” in 49th Annu. Allerton Conf. Commun., Control, Comput.,
Sept. 2011, pp. 1492–1499.
[21] W. Ma and S. Thorp, “An Efficient Algorithm to Locate All the Load
Flow Solutions,” IEEE Trans. Power Syst., vol. 8, no. 3, p. 1077, 1993.
[22] H. Chen, “Cascaded Stalling of Induction Motors in Fault-Induced
Delayed Voltage Recovery (FIDVR),” MS Thesis, Univ. Wisconsin–
Madison, ECE Depart., 2011. [Online]. Available: http://minds.
wisconsin.edu/handle/1793/53749
[23] D. Molzahn, B. Lesieutre, and H. Chen, “Counterexample to a
Continuation-Based Algorithm for Finding All Power Flow Solutions,”
IEEE Trans. Power Syst., vol. 28, no. 1, pp. 564–565, 2013.
[24] 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 Trans. Circuits Syst. I, Reg.
Papers, vol. 52, no. 3, pp. 625–630, Mar. 2005.
[25] B. Lesieutre and D. Wu, “An Efficient Method to Locate All the Load
Flow Solutions – Revisited,” in 53rd Annu. Allerton Conf. Commun.,
Control, Comput., Sept. 29 - Oct. 2 2015.
[26] H. Mori and A. Yuihara, “Calculation of Multiple Power Flow Solutions
with the Krawczyk Method,” in IEEE Int. Symp. Circuits Syst. (ISCAS),
30 May - 2 Jun. 1999.
[27] J. Ning, W. Gao, G. Radman, and J. Liu, “The Application of the
Groebner Basis Technique in Power Flow Study,” in North American
Power Symp. (NAPS), Oct. 2009.
[28] P. Dreesen and B. De Moor, “Polynomial Optimization Problems are
Eigenvalue Problems,” in Model-Based Control: Bridging Rigorous
Theory and Advanced Technology, P. Hof, C. Scherer, and P. Heuberger,
Eds. Springer US, 2009, pp. 49–68.
[29] J.-B. Lasserre, Moments, Positive Polynomials and Their Applications.
Imperial College Press, 2010, vol. 1.

Preprint No: ADP-15-34/T936

[30] D. Mehta, H. Nguyen, and K. Turitsyn, “Numerical Polynomial Homotopy Continuation Method to Locate All The Power Flow Solutions,”
arXiv:1408.2732, 2014.
[31] S. Chandra, D. Mehta, and A. Chakrabortty, “Exploring the Impact of
Wind Penetration on Power System Equilibrium Using a Numerical
Continuation Approach,” arXiv:1409.7844, 2014.
[32] D. Mehta, N. S. Daleo, F. Dörfler, and J. D. Hauenstein, “Algebraic Geometrization of the Kuramoto Model: Equilibria and Stability Analysis,”
Chaos, vol. 25, no. 5, 2015.
[33] S. Smale, “Mathematical Problems for the Next Century,” Math. Intell.,
vol. 20, no. 1, pp. 7–15, 1998.
[34] R. Guedes, L. Alberto, and N. Bretas, “Power System Low-Voltage
Solutions Using an Auxiliary Gradient System for Voltage Collapse
Purposes,” IEEE Trans. Power Syst., vol. 20, no. 3, pp. 1528–1537,
Aug. 2005.
[35] Y. Feng, “Solving for the Low-Voltage/Large-Angle Power-Flow Solutions by Using the Holomorphic Embedding Method,” dissertation,
Arizona State Univ., Aug. 2015, http://repository.asu.edu/attachments/
158002/content/Feng asu 0010E 15211.pdf.
[36] C. Tavora and O. J. Smith, “Equilibrium Analysis of Power Systems,”
IEEE Trans. Power App. Syst., vol. PAS-91, no. 3, pp. 1131–1137, 1972.
[37] J. Baillieul, “The Critical Point Analysis of Electric Power Systems,” in
IEEE 23rd Ann. Conf. Decis. Contr. (CDC), vol. 23, 1984.
[38] A. Klos and J. Wojcicka, “Physical Aspects of the Nonuniqueness of
Load Flow Solutions,” Int. J. Elec. Power Energy Syst., vol. 13, no. 5,
pp. 268–276, 1991.
[39] J. Baillieul and C. Byrnes, “Geometric Critical Point Analysis of
Lossless Power System Models,” IEEE Trans. Circuits Syst., vol. 29,
no. 11, 1982.
[40] A. J. Sommese, J. Verschelde, and C. W. Wampler, “Advances in
Polynomial Continuation for Solving Problems in Kinematics,” J. Mech.
Design, vol. 126, no. 2, pp. 262–268, 2004.
[41] D. Bates, J. Hauenstein, A. Sommese, and C. Wampler, Numerically
Solving Polynomial Systems with Bertini. SIAM, 2013, vol. 25.
[42] S. Guo and F. Salam, “The Number of (Equilibrium) Steady-State
Solutions of Models of Power Systems,” IEEE Trans. Circuits Syst. I,
Fundam. Theory Appl., vol. 41, no. 9, pp. 584–600, Sept. 1994.
[43] H.-D. Chiang and T. Wang, “On the Number of System Separations in
Power System,” in IEEE Int. Symp. Circuits Syst. (ISCAS), May 24-27
2015.
[44] S. Guo and F. Salam, “Determining the Solutions of the Load Flow of
Power Systems: Theoretical Results and Computer Implementation,” in
IEEE 29th Ann. Conf. Decis. Contr. (CDC), Dec. 1990, pp. 1561–1566.
[45] T. Y. Li, “Solving Polynomial Systems by the Homotopy Continuation
Method,” Handbook of numerical analysis, vol. 11, pp. 209–304, 2003.
[46] A. J. Sommese and C. W. Wampler, The Numerical Solution of Systems
of Polynomials Arising in Engineering and Science. World Scientific,
2005, vol. 99.
[47] D. Bates, J. Hauenstein, A. Sommese, and C. Wampler, “Bertini: Software for Numerical Algebraic Geometry,” Available at bertini.nd.edu,
doi: dx.doi.org/10.7274/R0H41PB5.
[48] A. Broder, “Generating Random Spanning Trees,” in 30th Annu. Symp.
Found. Comput. Sci., Oct. 1989, pp. 442–447.
[49] C. Bron and J. Kerbosch, “Algorithm 457: Finding All Cliques of an
Undirected Graph,” Commun. ACM, vol. 16, no. 9, pp. 575–577, Sept.
1973.
[50] R. Karp, Reducibility Among Combinatorial Problems. Springer, 1972.
[51] R. Zimmerman, C. Murillo-Sánchez, and R. Thomas, “MATPOWER:
Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education,” IEEE Trans. Power Syst., vol. 99, pp.
1–8, 2011.
[52] F. Harary, “A Characterization of Block-Graphs,” Canad. Math. Bull.,
vol. 6, no. 1, pp. 1–6, 1963.

