1

Robust Steady State Analysis of the Power Grid
Amritanshu Pandey1, Graduate Student Member, IEEE, Marko Jereminov1, Graduate Student Member,
IEEE, Martin R. Wagner1, Graduate Student Member, IEEE, David M. Bromberg1, Graduate Student
Member, IEEE, Gabriela Hug2, Member, IEEE and Larry Pileggi1, Fellow, IEEE

Abstract—A robust methodology for obtaining the steady-state
solution of the power grid is essential for reliable operation as well
as planning of the future transmission and distribution grid. At
present, disparate methods exist for steady-state analysis of the
transmission (power flow) and distribution power grid (threephase power flow). All existing alternating current (AC) power
flow and three-phase power flow analyses formulate a non-linear
problem that generally lacks the ability to ensure convergence to
the correct physical solution from an arbitrary set of initial
guesses. In this paper, we demonstrate that our equivalent circuit
formulation approach can model both the positive sequence model
of the transmission grid and three-phase model of the distribution
grid without loss of generality. Moreover, when combined with
novel circuit simulation techniques, it can robustly solve for the
steady-state solution for both these network models independent
of the choice of initial conditions. Examples for positive sequence
transmission and three-phase distribution systems, including 75k+
nodes Eastern Interconnection transmission test cases and 8k+
nodes taxonomy distribution test cases, are solved from arbitrary
initial guesses to demonstrate the robustness of our approach.
Index Terms— circuit simulation methods, continuation
methods, convergence problems, equivalent circuit approach,
power flow, robust convergence, steady-state analysis, three-phase
power flow, Tx stepping method

I. INTRODUCTION

A

n interconnected electric grid is a network of
synchronized power providers and consumers that are
connected via transmission and distribution lines and
operated by one of multiple entities. Robust operation of this
electric grid is of utmost importance for maintaining a country’s
economy and the well-being of its citizens.
To operate the grid robustly under all conditions, as well as
adequately plan for the future, it is a must that one can analyze
the grid robustly off-line and in real-time. At present, numerous
analysis methods exist for robust operation and planning of the
grid that can be broadly categorized into one of the following
categories: i) steady-state analysis in frequency domain (power
flow, three phase power flow, and harmonic analyses), ii)
transient and steady-state analysis in time domain, iii) analysis
for optimal dispatch of resources, and iv) other market dispatch
based analyses. Amongst these analyses, fundamental
frequency based steady-state analysis (power flow and threephase power flow) is essential for the day-to-day operation as
well as future planning of the grid. Furthermore, the steady-

state analysis also forms the backbone for other analyses such
as the optimal power flow problem and the state estimation.
At present, the steady-state simulation domain in power
industry is broadly divided spatially into high-voltage
transmission systems and sub-station level voltage distribution
systems. Disparate methods exist for analyzing these two
(transmission and distribution) systems. The steady-state
solution for the high voltage transmission system is obtained
via positive sequence AC power flow analysis (also referred to
as power flow analysis in this paper), whereas the steady-state
operating point for the distribution system is obtained via threephase AC power flow analysis. The industry standard for
solving the positive sequence AC power flow problem is the
‘PQV’ formulation [1], wherein nonlinear power mismatch
equations are solved for bus voltage magnitudes and angles that
further define the steady-state operating point of the system. On
the other hand, either the backward-forward sweep method [2]
or the current injection method (CIM) [3] are primarily used for
obtaining the steady-state solution of the three-phase power
flow problem for the distribution grid.
All of the existing methods for solving the positive sequence
power flow problem as well as the three-phase power flow
problem suffer from lack of robustness [5]-[6] in terms of
convergence for the nonlinear problem formulation. The ‘PQV’
based formulation for the positive sequence power flow
problem is known to diverge or converge to non-physical
solutions for ill-conditioned [2] and large scale (>50k buses)
systems [15]. Similarly, the backward-forward sweep method
that was proposed to solve radial and weakly meshed
distribution systems with high R/X ratio [2] has difficulties
converging for heavily meshed systems with more than a single
source of generation [7]. On the other hand, the CIM method
based on Dommel’s work in 1970 [4] has challenges with
incorporating multiple PV buses in the system [8]-[9]. Broadly
speaking, of the many known drawbacks of the existing power
flow and three-phase power flow analyses that are contributing
towards a lack of robustness, the two that are the most
detrimental are convergence to non-physical or unacceptable
solution [10] and divergence [5].
The factors that are the most fundamental toward making
these problems challenging are the use of non-physical
representations for modeling the power grid components, and
in the case of the ‘PQV’ formulation, the use of inherently nonlinear power mismatch equations to formulate the problem. The
non-physical models, such as the approximate model of a

1
Authors are with the Electrical and Computer Engineering Department, Carnegie Mellon University, Pittsburgh, PA 15213 USA, (e-mail: {amritanp, mjeremin,
mwagner1, dbromber, pileggi}@andrew.cmu.edu).
2
Author is with the Power System Laboratory, ETH Zurich, (e-mail: hug@eeh.ee.ethz.ch).

Submitted manuscript for IEEE Transactions on Power Systems (under review)

2

generator that is represented via positive sequence or threephase PV bus, can result in convergence to a non-physical
solution or divergence. Similarly, the inherent non-linearities in
the ‘PQV’ formulation almost always cause divergence for
large (>50k) and ill-conditioned test cases [15] when solved
from an arbitrary set of initial conditions. This lack of a physicsbased formulation, along with the methods that can constrain
the non-physics based models in their physical space, is what
renders the existing power flow and three-phase power flow
problem and solution approaches to be “non-robust.”
To develop a robust solver for the steady-state solution of
the power grid, it is imperative that the solver can efficiently
and effectively navigate through these challenges while
converging to a solution that is both meaningful and correct.
Intuitionally and physically, both the transmission and
distribution electric grids correspond to an electric circuit.
Therefore, our approach toward solving the power flow and
three-phase power flow problem is to treat them as such and
perform both analyses using circuit simulation fundamentals.
Toward this objective we propose a two-pronged approach.
First, the use of equivalent circuit formulation with true state
variables of currents and voltages [11]-[13] to model both the
transmission and distribution power grid (Sect. III.). Secondly,
the use of circuit simulation methods [14]-[15] to ensure robust
convergence to correct physical solutions (Sect. IV.).
II. BACKGROUND
A power grid in its simplest form can be represented by a set
of 𝒩 buses, where the sets of generators 𝒢 and load demands
ℒ are subsets of 𝒩, which are further connected by a set of line
elements, 𝒯 and set of transformers 𝑥𝑓𝑚𝑟𝑠 . Furthermore, there
is a set of slack buses (one for each island in the system)
represented by ϩ. In addition to these, the power grid may
contain other elements, such as shunts, flexible alternating
current transmission system (FACTS), etc. The objective of
steady-state analysis of the power grid is to model the
fundamental frequency component of the power grid and solve
for the complex voltages at its buses. The high voltage
transmission network of the grid generally operates under
balanced conditions, and therefore, the steady-state solution of
the transmission network is obtained via positive sequence
power flow model and analysis. In contrast, the distribution
network of the power grid can operate under unbalanced
conditions, and therefore we must apply three-phase power
flow model and analysis to find the steady-state solution of the
distribution grid. In the following sub-sections, we discuss the
current formulations used for steady-state analysis of
transmission and distribution networks and highlight their
limitations.
A. ‘PQV’ based Formulation for Positive Sequence Power
Flow Problem
The ‘PQV’ based power flow formulation is the industry
standard for solving for the steady-state solution of the high
voltage transmission network. In this formulation, a set of
2(𝒩 − ϩ) − 𝒢 power mismatch equations are solved for
unknown complex voltage magnitudes and angles of the system
using the Newton Raphson (NR) method. The set of power
mismatch equations are defined as follows:

𝒩

𝑃 − 𝑃 = |𝑉 |

|𝑉 | 𝐺 𝑐𝑜𝑠  + 𝐵 𝑠𝑖𝑛 

(1)

|𝑉 | 𝐺 𝑠𝑖𝑛  − 𝐵 𝑐𝑜𝑠 

(2)

𝒩

𝑄 − 𝑄 = |𝑉 |

where, 𝑃 + 𝑗𝑄 and 𝑃 + 𝑗𝑄 are the complex generation and
complex load at the node 𝑖 and 𝐺 + 𝑗𝐵 is the complex
admittance between the nodes 𝑖 and 𝑙.
In order to solve for unknown complex voltages 𝑉 ∠ in
the system, the real and reactive power mismatch equations
given by (1)-(2) are solved for the set of (𝒩 − ϩ − 𝒢 ) buses in
the system, whereas only real mismatch equations (1) are
solved for the set of buses with generators 𝒢 connected to it.
Importantly, this ‘PQV’ formulation is inherently non-linear,
since the set of network constraints result in non-linear power
mismatch equations independent of physics of the models used.
For example, in the ‘PQV’ formulation, a linear network
consisting of linear models for slack bus, transmission lines and
loads would correspond to a non-linear set of power mismatch
equations, a feature that could result in convergence difficulties
for systems even trivial in size.
B. Current Injection Method for Three-Phase Power Flow
Problem
Until recently, the backward forward sweep method was the
most commonly used method for the steady-state analysis of the
radial and weakly meshed distribution systems [2]. The method
was preferred over the ‘PQV’ formulation due to the radial
nature of the distribution grid and high R/X ratios of the
distribution lines, both of which are known to cause
convergence difficulties for the NR method [2]. However, the
backward forward sweep method itself is prone to convergence
difficulties for systems that are highly meshed or have multiple
sources [7].
The current injection method (CIM) for the three-phase
power flow problem [3] was proposed to address challenges
associated with the ‘PQV’ formulation and the backwardforward sweep method. In the CIM formulation, the non-linear
current mismatch equations for the system buses are solved via
NR method for each individual phase with complex rectangular
real and imaginary voltages (𝑉 + 𝑗𝑉 ) as the unknown
variables. The current mismatch equations for the three-phase
power flow problem are defined as follows [3]:
𝛥𝐼

=

𝑃

𝑉 + 𝑄
𝑉

𝑉

+ 𝑉

(3)

𝒩

−
𝛥𝐼 =

𝑃

𝐺 𝑉 −𝐵 𝑉

𝑉 − 𝑄
𝑉

𝑉

+ 𝑉

(4)

𝒩

−

𝐺 𝑉 −𝐵 𝑉

where 𝛥𝐼 + 𝑗𝛥𝐼 is the net current mismatch in phase 𝛺 at
node 𝑖 and 𝑃
+𝑗 𝑄
is the specified complex power
injection at node 𝑖. The set Ω includes phases a, b and c.
Although, the CIM method is known to improve the
convergence properties for heavily and weakly meshed three-

Submitted manuscript for IEEE Transactions on Power Systems (under review)

3

phase radial distribution systems with high R/X ratio, the
method is known to diverge for test-cases with high penetration
of PV buses [8]. Traditionally, the number of PV buses in the
distribution system were limited to a small number; however,
with the advent of large scale installation of distributed
generation (DGs) and voltage control devices in the distribution
system this is no longer true. Therefore, it is essential that a
robust three-phase power flow formulation can robustly handle
high penetration of PV buses and voltage control devices in the
system.

In case of three-phase power flow, three such additional
variables and constraints are introduced.
As an example, the equivalent circuit for the positive
sequence model for a PV bus used in power flow is shown in
Fig. 1 for the (k+1)th iteration of NR. It is constructed by
linearizing the set of equations (5)-(6) for the positive sequence
parameters and then representing the resulting equations using
fundamental circuit elements. To construct the PV bus
equivalent circuit for three-phase power flow problem, three
such circuits are first constructed and then are connected in
grounded-wye configuration.

III. EQUIVALENT CIRCUIT FORMULATION

Real Circuit

We apply the equivalent circuit approach in [11]-[15] for
steady-state analysis of the transmission and distribution power
grid to tackle the challenges exhibited by the existing
formulations. This approach for generalized modeling of the
power system in steady-state (i.e. power flow and three-phase
power flow) represents both the transmission and distribution
power grid elements in terms of equivalent circuit elements
without loss of generality. It was shown that each of the power
system components can be directly mapped to an equivalent
circuit model based on the underlying relationship between
current and voltage state variables. Importantly, this
formulation can represent any physics based model or
measurement based semi-empirical models as a sub-circuit, as
shown in [18], [21] and [22], and these models can be combined
hierarchically with other circuit abstractions to build larger
aggregated models. In the following section, we discuss
equivalent circuit representations of generic power system
components for both the positive sequence power flow problem
and the three-phase power flow problem. Note that throughout
the paper, the symbol superscript Ω in the mathematical
expressions of the equivalent circuit models represents a phase
from the set Ω of three phases a, b and c for the three-phase
problem and represents the positive sequence (p) component for
the power flow problem.
A. PV Bus or the Generator Model
In the equivalent circuit approach, the generator (PV) bus
model is modeled via a complex current source [14]. To enable
the application of NR, this complex current source is split into
real and imaginary current sources (𝐼 and 𝐼Ω𝐼𝐺 , respectively).
This is necessary due to the non-analyticity of complex
conjugate functions [11]. The resulting equations for the PV
model in the power flow and three-phase power problem are:
𝐼

=

𝐼

=

𝑃 𝑉

+𝑄 𝑉

(𝑉 ) + (𝑉 )
𝑃 𝑉 −𝑄 𝑉
(𝑉 ) + (𝑉 )

(5)

Imag. Circuit

𝐼

𝐼

+

+

𝑉

𝑉
_

_

Figure 1: Equivalent Circuit Model for PV generator model.

B. Voltage Regulation of a Bus
Numerous power grid elements such as generators, FACTS
devices, transformers, shunts etc., can control a voltage
magnitude at a given node in the system. Moreover, they can
control the voltage magnitude (𝑉 ) at either their own node 𝒪
or a remote node 𝒲 in the system. In equivalent circuit
formulation, we represent the control of the voltage magnitude
by a control circuit (Fig. 2), which is governed by the following
expression:
𝐹𝒲 ≡ (𝑉 ) − (𝑉

𝒲

) − (𝑉 𝒲 ) = 0

(7)

The circuit in Fig. 2 is derived from the linearized version
of (7). For the power flow problem, it is stamped for each node
𝒲 in the system whose voltage is being controlled such that
there exists at least one single path between the node 𝒲 and the
equipment’s node 𝒪 that is controlling it. Similarly, for threephase power flow three of these circuits are stamped for each
node 𝒲. The additional unknown variables for these additional
constraints are dependent on the power system device that is
controlling the voltage magnitude. For example, the additional
unknown variable for a generator is its reactive power 𝑄 ,
whereas in the case of transformers, it is the transformer tap
𝑡𝑟 , and for FACTS devices it can be the firing angle 𝜑 . The
previous section already described how the additional unknown
variable 𝑄 for PV buses is integrated into the respective
equivalent circuits for generators.

(6)

Additional constraints that allow the generators to control the
voltage magnitude either at its own node or any other remote
node in the system are modeled by a control circuits, as shown
in the following subsection. In the case of power flow problem,
a single control circuit is needed whereas for the three-phase
power flow problem, three such control circuits are needed for
each PV bus in the system. The reactive power 𝑄 of the
generator acts as the additional unknown variable for the
additional constraint that is introduced due to voltage control.

𝑉

Figure 2: Voltage magnitude constraint control equivalent circuit.

C. ZIP Load Model
In this section, we derive the positive sequence and threephase model for the ZIP load. The ZIP load models the

Submitted manuscript for IEEE Transactions on Power Systems (under review)

4

aggregated load in the system as a mix of constant impedance
𝑍 + 𝑗𝑍 , constant current 𝐼 + 𝑗𝐼 , and constant
power 𝑆 + 𝑗𝑆 load behavior, which can be mathematically
represented as follows:
𝑃

=𝑍

𝑉

+𝐼

𝑉

+𝑆

𝑄

=𝑍

𝑉

+ 𝐼

𝑉

+𝑆

(8)
(9)

In the equivalent circuit approach, the equations for the ZIP
load model can be re-written as:
𝐼

𝑆 𝑉 + 𝑆 𝑉

= 𝑌 𝑉 −𝑌 𝑉 +

𝑉
+
𝐼

𝐼

+ 𝐼

𝑆 𝑉 − 𝑆 𝑉
𝑉

𝐼

+ 𝐼

(10)

. 𝑐𝑜𝑠( − 𝐼 )

= 𝑌 𝑉 +𝑌 𝑉 +
+

+ 𝑉

+ 𝑉

(11)

. 𝑠𝑖𝑛( − 𝐼 )

where:
𝐼

𝐼

= 𝑡𝑎𝑛

 = 𝑡𝑎𝑛
𝑌𝛺𝑃 + 𝑗𝑌𝛺𝑄 =

(12)

𝐼
𝑉
𝑉

(13)

1
𝑍𝛺𝑃

(14)

+ 𝑗𝑍𝛺𝑄

circuit theoretic approach in [21]-[24] and aims to create a
linear load model that can capture the true measure and
sensitivity of the aggregated load in the system. The complex
governing equation of the generalized load current for the BIG
load model is given by:
𝐼
+

+𝑗 𝐼

𝑉

_

+

+

_

Figure 3: Real circuit for a) wye connected ZIP Load Model (on left) b) delta
(D) connected ZIP load model (on right).

It is important to note that both the ZIP and PQ load models
result in non-linear network constraints for both the ‘PQV’ and
CIM formulations. In the ‘PQV’ formulation the non-linearities
in the network constraints are due to the use of power mismatch
equations whereas in the CIM, the non-linearities are due to PQ
and ZIP model equations. These added non-linearities are one
of the primary causes of divergence and convergence to nonphysical solutions with these formulations. To address this
problem, we have proposed an accurate and yet linear BIG load
model [21]-[23].
D. BIG Linear Load Model
The BIG aggregated load model was proposed based on the

+𝑗 𝛼

((𝐺

) + 𝑗(𝐵

(15)

) )

where 𝛼
+𝑗 𝛼
represents the base value for the
modeled aggregated load and the corresponding complex
admittance ((𝐺 ) + 𝑗(𝐵 ) ) captures the voltage
sensitivities. For instance, a negative conductance in
conjunction with complex current 𝛼
+𝑗 𝛼
mimics
the inverse current/voltage sensitivity relationship, similar to
constant power (PQ) load behavior and positive conductance in
conjunction with complex current source will represent the
positively correlated current/voltage sensitivity relationship,
similar to the impedance load behavior. Both the positive and
negative impedances capture the change in load with voltage
with respect to the portion of the load that is modeled by the
current source. Fig. 4 shows the positive sequence (p) BIG load
model. Similar to the ZIP load model, the three-phase BIG load
model can be constructed by connecting the equivalent circuits
of individual phases in wye or delta formation.
Real Circuit

𝐼

For the load model given in (10) through (14), the constant
impedance part of the load is linear, whereas the constant
current and constant power part of the aggregated load is
nonlinear. Once, (10)-(11) are linearized, they are used to
construct the equivalent circuit models for both the power flow
and three-phase power flow problem. The constructed threephase model of the ZIP load model can either be connected in
wye or delta formation. As an example, ZIP load model
connected in wye and delta formation is shown in Fig. 3.

= 𝛼

+𝑗 𝑉

Imag. Circuit

𝐼

+

+

𝑉

𝑉
_

𝐺

−𝐵

α

_

𝐺

𝐵

α

Figure 4: Equivalent circuit of a BIG load model

IV. CIRCUIT SIMULATION METHODS
Decades of research in circuit simulation have demonstrated
that circuit simulation methods can be applied for determining
the DC state of highly non-linear circuits using NR. These
techniques have been shown to make NR robust and practical
for large-scale circuit problems [19] consisting of billions of
nodes. Most notable is the ability to guarantee convergence to
the correct physical solution (i.e. global convergence) and the
capability of finding multiple operating points [25]. We have
previously proposed analogous techniques for ensuring
convergence to the correct physical solution for the power flow
problem [14]-[15]. In this section, we provide a short overview
of these techniques for the power flow problem and further
extend these techniques for their use in the three-phase power
flow problem.
A. General Methods
1) Variable Limiting
The solution space of the system node voltages in a power
flow and three-phase power flow problem are well defined.
While solving these problems, a large NR step may step out of
this solution space and result in either divergence or
convergence to a non-physical solution. It is, therefore,
important to limit the NR step before an invalid step out of the
solution space is made. In [14], we proposed the variable
limiting method to achieve the postulated goal. In this

Submitted manuscript for IEEE Transactions on Power Systems (under review)

5

technique, the state variables that are most sensitive to initial
guesses are damped when the NR algorithm takes a large step
out of the pre-defined solution space. Note, however, that not
all the system variables are damped for the variable limiting
technique, as is done for traditional damped NR. Circuit
simulation research has shown that damping the most sensitive
variables provides superior convergence compared to damped
NR in general [19].
To apply variable limiting in our prototype simulator for the
power flow and three-phase power flow problem, the
mathematical expressions for the PV nodes in the system are
modified as follows:
𝐼

= 𝜍

𝜕𝐼
𝜕𝑉

𝑉

−𝑉

+𝐼

∆

𝜕𝐼
+𝜍
𝜕𝑉

𝑉

𝜕𝐼
+
𝜕𝑄

−𝑉

(16)
𝑄

−𝑄

∆

where, 0 ≤ ς ≤ 1. The magnitude of ς is dynamically varied
through heuristics such that convergence to the correct physical
solution is achieved in the most efficient manner. The heuristics
depend on the largest delta voltage (∆𝑉𝛺𝑅𝐺 , ∆𝑉𝛺𝐼𝐺 ) step during
subsequent NR iterations. If during subsequent NR iterations, a
large step (∆𝑉𝛺𝑅𝐺 , ∆𝑉𝛺𝐼𝐺 ) is encountered, then the factor ς is
decreased. The factor ς is scaled back up if consecutive NR
steps result in monotonically decreasing absolute values for the
largest error.
2) Voltage Limiting
An equally simple, yet effective, technique is to limit the
absolute value of the delta step that the real and imaginary
voltage vectors can make during each NR iteration. This is
analogous to the voltage limiting technique used for diodes in
circuit simulation, wherein the maximum allowable voltage
step during NR is limited to twice the thermal voltage of the
diode [20]. Similarly, for the power flow and three-phase power
flow analyses, a hard limit is enforced on the normalized real
and imaginary voltages in the system. The mathematical
implementation of voltage limiting in our formulation is as
follows:
(𝑉 )

= 𝑚𝑖𝑛 𝑚𝑎𝑥 (𝑉 ) + 𝛿 𝑚𝑖𝑛 ∆(𝑉 ) , ∆𝑉

𝑚𝑖𝑛𝑚𝑎𝑥 =

𝑉
𝑉

,
,
𝑥,

𝑖𝑓 𝑥 > 𝑉
𝑖𝑓 𝑥 < 𝑉
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

(17)
(18)

and 𝛿 = 𝑠𝑖𝑔𝑛 ∆ 𝑉
and 𝐶 ∈ {𝑅, 𝐼} represents the
placeholder for real and imaginary parts.
Analogously, other system variables such as the reactive
power 𝑄 of the PV buses, can be limited by limiting the
calculated currents 𝐼 + ∆ 𝐼
at NR step k + 1 and then
finding the new 𝑄
from inverse function (𝑓 ) of (𝐼Ω𝐶 +
∆(𝐼Ω𝐶 )k ).
B. Homotopy Methods
Limiting methods may fail to ensure convergence for certain
ill-conditioned and large test systems when solved from an
arbitrary set of initial guesses. To ensure convergence for these
network models to the correct physical solutions independent
of the choice of initial conditions, we propose the use of
homotopy methods. Homotopy methods that have been

proposed in the past [5], [16]-[17] have suffered from
convergence to low voltage solutions and divergence.
Furthermore, none of the previously proposed homotopy
methods are known to scale up to test systems that are of the
scale of the European or the US grids and in general are not
extendable to the three-phase power flow problem.
In the homotopy approach, the original problem is replaced
with a set of sub-problems that are sequentially solved. The set
of sub-problems exhibit certain properties, namely, the first
sub-problem has a trivial solution and each incrementally
subsequent problem has a solution very close to the solution of
the prior sub-problem. Mathematically this can be described via
the following expression:
ℋ(𝑥, 𝜆) = (1 − 𝜆)ℱ(𝑥) + 𝜆𝒢(𝑥)

(19)

where 𝜆  [0, 1].
The method begins by replacing the original problem
ℱ(𝑥) = 0 with ℋ(𝑥, 𝜆) = 0. The equation set 𝒢(𝑥) is a
representation of the system that has a trivial solution. The
homotopy factor 𝜆 has the value of 1 for the first sub-problem
and therefore the initial solution is equal to trivial solution
of 𝒢(𝑥). For the final sub-problem that corresponds to the
original problem, the homotopy factor 𝜆 has the value of zero.
To generate sequential sub-problems, the homotopy factor is
dynamically decreased in small steps until it has reached the
value of zero.
In this paper, we discuss two homotopy methods that were
specifically developed for the power systems analysis:
1) Tx Stepping
We proposed the “Tx Stepping” method in [15] specifically
for the power flow problem. In this section, the method is
briefly discussed and is further extended for the three-phase
power flow problem.
a) General Approach
In Tx stepping method, the series elements in the system
(transmission lines, transformers etc.) are first “virtually”
shorted to solve the initial problem that has a trivial solution.
Specifically, a large conductance (≫ 𝐺 ) and a large
susceptance (≫ 𝐵 ) are added in parallel to each transmission
line and transformer model in the system. In case of three-phase
power flow, a large self-impedance (≫ Y𝑖𝑙ΩΩ ) is added in parallel
to each phase of the transmission line and transformer model.
Furthermore, the shunts in the system, are open-circuited by
modifying the original shunt conductance and susceptance
values. Importantly, the solution to this initial problem results
in high system voltages (magnitudes), as they are essentially
driven by the slack bus complex voltages and the PV bus
voltage magnitudes due to the low voltage drops in the lines and
transformers (as expected with virtually shorted systems).
Similarly, the solution for the bus voltage angles will lie within
an -small radius around the slack bus angle. Subsequently, like
other continuation methods, the formulated system problem is
then gradually relaxed to represent the original system by
taking small increment steps of the homotopy factor (𝜆) until
convergence to the solution of the original problem is achieved.
Mathematically, the line and transformer impedances during
homotopy for the power flow is expressed by:
𝑖𝑙 ∈ {𝒯 , 𝑥𝑓𝑚𝑟𝑠} ∶ 𝐺 + 𝑗𝐵 = (𝐺 + 𝑗𝐵 )(1 + 𝜆𝛾) (20)

Submitted manuscript for IEEE Transactions on Power Systems (under review)

6

and for the three-phase problem:
𝑌
𝑌
𝑌

𝑌
𝑌
𝑌
𝑌
𝑌
𝑌
Y (1 + 𝛾𝜆)
=
Y
Y

Y
Y (1 + 𝛾𝜆)
Y

Y
Y
Y (1 + 𝛾𝜆)

(21)

where, 𝐺 , 𝐵 , and 𝑌ΩΩ 𝑖𝑙 are the original system impedances
and 𝐺 , 𝐵 , and 𝑌ΩΩ are the system impedances used while
iterating from the trivial problem to the original problem. The
parameter 𝛾 is used as a scaling factor for the conductances and
susceptances. If the homotopy factor (𝜆) takes the value one,
the system has a trivial solution and if its takes the value zero,
the original system is represented.
Along with ensuring convergence for a problem, Tx stepping
avoids the undesirable low voltage solutions for the positive
sequence power flow and three-phase power flow problem
since the initial problem results in a solution with high system
voltages, and each subsequent step of the homotopy approach
continues and deviates ever so slightly from this initial solution,
thereby guaranteeing convergence to the high voltage solution
for the original problem.
b) Handling of Transformer Phase Shifters and Taps
To “virtually short” a power system, we must also account
for transformer taps 𝑡𝑟 and phase shifting angles 𝛩 . In a
“virtually” shorted condition, all the nodes in the system must
have complex voltages that are near the slack bus or PV bus
complex voltages, which can be intuitively defined by a small
epsilon norm ball around these voltages. Therefore, to achieve
the following form, we must modify the transformer taps and
phase shifter angles such that at 𝜆 = 1, their turns ratios and
phase shift angles correspond to a magnitude of 1 pu and 0° ,
respectively. Subsequently, the homotopy factor 𝜆 is varied
such that the original problem is solved with original
transformer tap and phase shifter settings. This can be
mathematically expressed as follows:
𝑖 ∈ 𝑥𝑓𝑚𝑟𝑠 ∶ 𝑡𝑟 = 𝑡𝑟 + 𝜆 1 − 𝑡𝑟
(22)
𝑖 ∈ 𝑥𝑓𝑚𝑟𝑠 ∶ 𝛩 = 𝛩 − 𝜆𝛩
(23)
c) Handling of Voltage Control for Remote Buses
To achieve a trivial solution during the first step of Tx
stepping it is essential that we also handle remote voltage
control appropriately. Remote voltage control refers to a device
on node 𝒪 in the system controlling the voltage of another node
𝒲 in the system. This behavior is highly non-linear and if not
handled correctly can result in divergence or converge to low
voltage solution. Existing commercial tools for power flow and
three-phase power flow analyses have difficulties dealing with
this problem and suffer from lack of robust convergence when
modeling remote voltage control in general. With Tx stepping
we can handle this problem efficiently and effectively. We first
incorporate a “virtually short path” between the controlling
node (𝒪) and the controlled node (𝒲) at 𝜆 = 1, such that the
device at the controlling node can easily supply the current
needed for node 𝒲 to control its voltage. Then following the
homotopy progression, we gradually relax the system such that
additional line connecting the controlling node (𝒪) and
controlled node (𝒲) is open at 𝜆 = 0.

d) Implementation of Tx Stepping in Equivalent Circuit
Formulation
Unlike traditional implementations of homotopy methods, in
our equivalent circuit formulation we do not directly modify the
non-linear set of mathematical equations, but instead embed a
homotopy factor in each of the equivalent circuit models for the
power grid components. In doing so we allow for incorporation
of any power system equipment into the Tx stepping approach
within our equivalent circuit formulation framework, without
loss of generality. Furthermore, we ensure, that the physics of
the system is preserved while modifying it for the homotopy
method. Fig. 5 demonstrates how the homotopy factor is
embedded into the equivalent circuit of the transformer.

Figure 5: Homotopy factor embedded in transformer equivalent circuit.

2) Dynamic Power Stepping
Another homotopy technique that can ensure robust
convergence for systems that have a low percentage of constant
voltage nodes in the system is the dynamic power stepping
method. Existing distribution systems tend to belong to these
types of systems and therefore, dynamic power stepping can be
applied to robustly obtain the steady-solution of the distribution
grid by solving the three-phase power flow problem. This
method has been previously described for the positive-sequence
power flow problem in [14] and is analogous to the source
stepping and gmin stepping approaches in standard circuit
simulation solvers.
In the power stepping method, the system loads and
generation are scaled back by a factor of 𝛽 until the
convergence is achieved. If these loads and generations are
scaled down all the way to zero, then the constraints for the PQ
buses in the system result in linear network constraints.
Similarly, current source non-linearities of the PV buses that
are due to the constant real power are also eliminated.
Therefore, by applying the power stepping factor, the nonlinearities in the system are greatly eased and convergence is
easily achieved. Upon convergence, the factor is gradually
scaled back up to unity to solve the original problem. In this
method, as in all continuation methods, the solution from the
prior step is used as the initial condition for the next step. The
mathematical representation of dynamic power stepping for the
three-phase power flow and positive sequence power flow
problem is as follows:
𝐺 ∈ 𝑃𝑉: 𝑃 = 𝛽𝑃
𝐿 ∈ 𝑃𝑄: 𝑃 = 𝛽𝑃 𝑎𝑛𝑑 𝑄 = 𝛽𝑄

(24)
(25)

where, PQ are all load nodes and PV are all generator nodes.
V. RESULTS
Example cases were simulated in our prototype solver
SUGAR (Simulation with Unified Grid Analyses and
Renewables) to demonstrate that the equivalent circuit
approach along with circuit simulation techniques facilitates a
robust framework for positive sequence power flow and the

Submitted manuscript for IEEE Transactions on Power Systems (under review)

7

three-phase power flow analyses. The example cases for
positive sequence power flow analyses include known illconditioned test cases and large network models that represent
different operating conditions of the eastern interconnection of
the US grid. For the three-phase power flow analysis, example
cases include a set of standard distribution taxonomy cases [26]
and meshed transmission grid test cases that were modified
from positive sequence to three-phase network models. The
results that follow demonstrate that the proposed framework
along with the use of circuit simulation methods can ensure
convergence to a correct physical solution for all the power
flow and three-phase power flow cases, independent of the
choice of the initial guess and thus overcomes the challenges
faced by existing formulations.
A. Positive Sequence Power Flow Results
1) Ill-Conditioned and Large Test cases
A convergence sweep was run on the ill-conditioned 13659
bus PEGASE test case using the SUGAR solver and a standard
commercial tool and their results were compared. Fig. 6 shows
that SUGAR was able to robustly converge to a correct physical
solution independent of the choice of the initial conditions,
whereas the standard tool was highly sensitive to the choice of
the initial guess and could converge to the correct physical
solution only from a few samples for the initial guess.

Figure 6: Convergence sweep comparison for 13659 node PEGASE testcase
between SUGAR and Standard tool with ▽ showing divergence and *
showing convergence

A similar convergence sweep was performed for larger test
cases (> 75k+ nodes) that represent different loading and
operating scenarios for eastern interconnection of the US grid.
Simulations were performed on three different test cases for 15
different initial conditions each. Results are shown in Table 1.
The set of initial conditions for all buses were identical and
were uniformly sampled from:
𝑉
∈ [−40, 40] , 𝑉
∈ [0.9, 1.1].
(26)
TABLE 1: CONVERGENCE PERFORMANCE FOR LARGE EASTERN
INTERCONNECTION TEST CASES
SUGAR
Standard Tool
Case
# Nodes
#
#
#
#
Name
Converge Diverge Converge Diverge
2017SLL
80778
0
15
15
0
Eastern1
76228
0
15
15
0
2021SUM
81904
0
15
15
0

For the larger eastern interconnection test cases, the runtime
per iteration is less than 0.4 seconds and is comparable to other
simulation tools out in the market. The total computation time
in general is dependent on the choice of initial conditions. A
sufficiently close initial condition may result in convergence

within 7 iterations whereas a totally random set of initial
guesses may take up to 100 iterations with Tx stepping method.
2) Contingency Analysis
In the next set of results, we performed a set of contingency
analyses with our tool SUGAR and a standard commercial tool
for two test cases that represent different network configuration
of the eastern interconnection of the US grid. The initial guess
for solving the contingency cases was chosen to be the
operating point prior to the contingency. The set of
contingencies in the experiment includes loss of generation
(ℒ ) and loss of branches (ℒ ). The results are summarized in
Table 2, which highlight the need for continuation methods to
solve such problems robustly.
TABLE 2: COMPARISON OF CONTINGENCIES OF LARGE TEST CASES
Case
# Nodes Contingency*
Standard Tool
SUGAR
2ℒ
Diverge
Converge
Case 1
76228
Diverge
Converge
2ℒ + 2ℒ
2ℒ
Diverge
Converge
Case 2
78201
2ℒ + 2ℒ
Diverge
Converge
*The number in front of ℒ and ℒ represents the equipment outage count. (For
e.g. 2ℒ represents that two generators were lost during this contingency).

B. Three-Phase Power Flow Results
1) Taxonomical Test Cases
Table 3 documents the results obtained from SUGAR threephase solver for standard taxonomical cases along with an
additional large meshed test case. The results validate the
accuracy of the SUGAR solver by comparing the results
obtained against a standard distribution power flow tool
GridLAB-D. Slight differences (less than 1e-3) in results for
taxonomical cases between SUGAR and GridLAB-D are due to
default values used for unspecified parameters (e.g. neutral
conductor resistance) in GridLAB-D. In order to demonstrate
that SUGAR can also solve a large meshed three-phase test
systems, an additional case was run. For this test case, we
extended the positive sequence model of the 9241-node
PEGASE test case to a three-phase model and documented the
simulation results for the case with SUGAR. This example did
not converge in GridLAB-D.
TABLE 3: SUGAR THREE-PHASE RESULTS FOR TAXONOMICAL CASES
Deviation from GridLAB-D
Iter.
Cases
#Nodes
Max. ΔVmag
Max. ΔVang
Count
[pu]
[°]
R1-12.47-1
2455
5
8.73E-04
9.94E-03
R2-12.47-3
2311
5
6.56E-04
1.32E-02
R3-12.47-3
7096
5
1.94E-03
3.89E-02
R4-12.47-1
2157
5
6.81E-04
9.61E-03
R5-12.47-5
2216
5
5.44E-05
4.20E-03
case9241pegase*
12528
5
NA#
NA#
* 9241 bus PEGASE transmission test case was extended to three-phase model
#The following case did not run in GridLAB-D

2) Ill-Conditioned Test Cases
To solve certain hard to solve ill-conditioned three-phase
test cases, we made use of homotopy methods. To demonstrate
one such example, we extended the standard 145 bus
transmission system model into a balanced three-phase network
model. Fig. 7 plots the convergence results for this test case
with and without the use of dynamic power stepping. Without
the use of dynamic power stepping, the test system did not
converge within maximum number of allowable iterations;
however, with the use of dynamic power stepping, the system
robustly converged to correct physical solution.

Submitted manuscript for IEEE Transactions on Power Systems (under review)

8

[9]
[10]

[11]

[12]

Figure 7: Convergence of 145 bus test case for three-phase power flow with
(bottom) and without (top) power stepping. For the power stepping case, the
green dotted line represents the change in continuation factor λ

VI. CONCLUSIONS
In this paper, we have demonstrated that an equivalent
circuit approach with the use of novel circuit simulation
methods can robustly solve for the steady-state solution of the
transmission and distribution grid without loss of generality.
This proposed formulation and the analogous circuit simulation
methods can be generically applied to both the positive
sequence power flow problem and the three-phase power flow
problem. Importantly, our approach toward steady-state
analyses of transmission and distribution grid ensures robust
convergence to correct physical solutions, and in doing so
enables robust contingency analyses, statistical analyses, and
security constrained optimal power flow analyses.

[13]

[14]
[15]

[16]
[17]

[18]

VII. REFERENCES
[1]
[2]

[3]

[4]

[5]

[6]
[7]

[8]

W. F. Tinney and C. E. Hart, “Power flow solutions by
Newton’s method,” IEEE Transactions on PAS, Vol. 86, No.
11, pp. 1449-1460, Nov. 1967.
D. Shiromohammadi, H. W. Hong, A. Semlyen, G. X. Luo, “A
compensation-based power flow method for weakly meshed
distribution and transmission networks”, IEEE Transactions
on Power Systems, Vol. 3, No. 2, May 1988.
P. A. N. Garcia, J. L. R. Pereria, S. Carneiro Jr., M. P. Vinagre,
F. V. Gomes, “Improvements in the Representation of PV
Buses on Three-Phase Distribution Power Flow”, IEEE
Transactions on Power Systems, Vol. 19, No. 2, April 2004.
H.W. Dommel, W.F. Tinney, and W.L. Powell, “Further
developments in Newton's method for power system
applications,” IEEE Winter Power Meeting, Conference Paper
No. 70 CP 161-PWR New York, January 1970.
W. Murray, T. T. De Rubira, and A. Wigington, “Improving
the robustness of Newton-based power flow methods to cope
with poor initial conditions”, North American Power
Symposium (NAPS), 2013.
A. G. Exposito, E. R. Ramos, “Reliable Load Flow Technique
for Radial Distribution Networks”, IEEE Transactions on
Power Systems, Vol. 14, No. 3, August 1999.
J. Liu, M. M. A. Salama, R. R. Mansour, “An efficient power
flow algorithm for distribution systems with polynomial load”,
International Journal of Electrical Engineering Education,
39/4.
Whei-Min Lin, Tung-Sheng Zhan, Ming-Tong Tsay, “Multiple
–Frequency Three-Phase Load Flow for Harmonic Analysis,”
IEEE Transactions on Power Systems, Vol. 19, No. 2, May
2004.

[19]
[20]
[21]

[22]

[23]

[24]

[25]

[26]

Mary B. Cain et al.,”History of Optimal Power Flow and
Formulations,” Federal Energy Regulatory Commission,
December 2012.
L. DeMarco, T. J. Overbye, “Low voltage power flow solutions
and their role in exit time-based security measures for voltage
collapse”, 27th Conference of Decision and Control, Austin,
Texas, Dec. 1988.
D. Bromberg, M. Jereminov, L. Xin, G. Hug, L. Pileggi, “An
Equivalent Circuit Formulation of the Power Flow Problem
with Current and Voltage State Variables”, PowerTech
Eindhoven, June 2015.
M. Jereminov, D. M. Bromberg, A. Pandey, L. Xin, G. Hug, L.
Pileggi, “An equivalent circuit formulation for three-phase
power flow analysis of distribution systems” T&D Conference,
2016 IEEE PES.
M. Jereminov, D. M. Bromberg, L. Xin, G. Hug, L. Pileggi,
“Improving Robustness and Modeling Generality for Power
Flow Analysis,” T&D Conference and Exposition, 2016 IEEE
PES.
A. Pandey, M. Jereminov, G. Hug, L. Pileggi, “Improving
Power Flow Robustness via Circuit Simulation Methods,”
IEEE PES General Meeting, Chicago, 2017.
A. Pandey, M. Jereminov, G. Hug, L. Pileggi, “Robust
Convergence of Power Flow using Tx Stepping Method with
Equivalent Circuit Formulation,” Power Systems Computation
Conference, 2018. (Accepted to PSCC 2017)
S. Cvijic, P. Feldmann, M. Ilic, “Applications of homotopy for
solving AC power flow and AC optimal power flow”, IEEE
PES General Meeting, San Diego, July 2012.
Mehta, H. D. Nguyen, K. Turitsyn, “Numerical Polynomial
Homotopy Continuation Method to Locate All the Power Flow
Solutions”, IET Generation, Transmission and Distribution,
Vol. 10, No. 12, August 16.
A. Pandey, M. Jereminov, X. Li, G. Hug, L. Pileggi, “Unified
Power System Analyses and Models using Equivalent Circuit
Formulation,” IEEE PES Innovative Smart Grid Technologies,
Minneapolis, USA, 2016.
L. Pillage (Pileggi), R. Rohrer, C. Visweswariah, Electronic
Circuit & System Simulation Methods, McGraw-Hill, Inc.,
New York, NY, USA, 1995.
L. Nagel, SPICE2: A Computer Program to Simulate
Semiconductor Circuits, EECS Department, University of
California Berkeley, 1975.
M. Jereminov, A. Pandey, H. A. Song, B. Hooi, C. Faloutsos,
L. Pileggi “Linear load model for robust power system
analysis”, IEEE PES Innovative Smart Grid Technologies,
Torino Italy, September 2017.
A. Pandey, M. Jereminov, X. Li, G. Hug, L. Pileggi,
“Aggregated Load and Generation Equivalent Circuit Models
with Semi-Empirical Data Fitting,” IEEE Green Energy and
Systems Conference (IGESC 2016) November 2016.
H. A. Song, B. Hooi, M. Jereminov, A. Pandey, L. Pileggi, C.
Faloutsos,” PowerCast: Mining and forecasting power grid
sequences” Joint European Conference on Machine Learning
and Knowledge Discovery in Databases, Springer 2017.
B. Hooi, H. A. Song, A. Pandey, M. Jereminov, L. Pileggi, and
C. Faloutsos. StreamCast: Fast and Online Mining of Power
Grid Time Sequences. Submitted to Proceedings of the 2018
SIAM International Conference on Data Mining. Society for
Industrial and Applied Mathematics.
L. B. Goldgeisser and M. M. Green, “Using Continuation
Methods to Improve Convergence of Circuits with High
Impedance Nodes,” in IEEE International Symposium on
Circuits and Systems, 2000, Geneva Switzerland.
Modern Grid Initiative: Distribution Taxonomy Final Report,
Pacific Northwest National Laboratory, November 15th, 2008.

Submitted manuscript for IEEE Transactions on Power Systems (under review)

