ExaGridPF: A Parallel Power Flow Solver for
Transmission and Unbalanced Distribution Systems
Bin Wang, John Bachan, Cy Chan
Lawrence Berkeley National Laboratory
{wangbin, jdbachan, cychan}@lbl.gov
Abstract—This paper investigates parallelization strategies for
solving power flow problems in both transmission and
unbalanced, three-phase distribution systems by developing a
scalable power flow solver, ExaGridPF, which is compatible with
existing high-performance computing platforms. NewtonRaphson (NR) and Newton-Krylov (NK) algorithms have been
implemented to verify the performance improvement over both
standard IEEE test cases and synthesized grid topologies. For
three-phase, unbalanced system, we adapt the current injection
method (CIM) to model the power flow and utilize SuperLU to
parallelize the computing load across multiple threads. The
experimental results indicate that more than 5 times speedup ratio
can be achieved for synthesized large-scale transmission
topologies, and significant efficiency improvements are observed
over existing methods for the distribution networks.
Index Terms—Power Flow, parallel computing, Newton-Raphson,
Newton-Krylov, High performance computing (HPC).

I.

INTRODUCTION

Solving power flow is the backbone for many applications
in power system analysis, such as planning, operation, transient
analysis, and contingency analysis. The objective of solving
power flow is to determine the voltage profile, the current flow,
and power loss, etc., across the power grid topology. Previous
research has focused on modeling approaches [1], [2] in
different systems and convergence properties of different
algorithms, using multi-core processors or graphics processing
units (GPU) [3]. As smart grid technologies advance, the power
system is becoming more observable, enabling modeling and
planning over larger networks with more details. The increased
penetration of distributed energy resources (DERs) with their
associated intermittency and uncertainty [4], such as electric
vehicles (EVs) [5]–[9], solar PV generation [4], battery energy
storage system (BESS) [4], [10], etc, will require operators to
solve the power flow problem over large-scale grid topologies
much more frequently to minimize grid constraint violations.
Essentially, the power flow problem solves the nonlinear
power system equations due to the nonlinear properties of the
ZIP load model [11], which is widely utilized by many power
flow applications. Iterative algorithms have been proposed to
solve the nonlinear equations, such as Gauss-Seidel [12],
Newton-Raphson [12], Forward-Backward-Sweep [11], Z-Bus
This paper is sponsored by the LDRD funding by Lawrence Berkeley
National Laboratory

method [2], etc. According to [12], Newton’s method is
numerically superior to the Gauss-Seidel method and less prone
to divergence for ill-conditioned problems, because of its
quadratic convergence properties. For transmission systems,
single line models are often utilized because most transmission
lines are balanced, so they can be modeled using a simplified
single-phase, balanced model. Voltage magnitude and angle are
utilized to model the single-phase system. On the other hand,
three-phase, unbalanced networks are very common in the
distribution system, where the single-phase, balanced model
assumption is no longer valid for the power flow calculation.
Open-source modeling tools have been developed to solve
the power flow problems on different platforms, such as
GridLab-D [13], OpenDSS [14], MAT-POWER [15], etc.
Those applications are mostly desktop-based with limited
computing resources. Newton-Raphson power flow algorithms
have been implemented in Gridlab-D based on current injection
method (CIM) using SuperLU [16] for solving sparse linear
equations. In addition, three-phase Forward-Backward-Sweep
(FBS) has been implemented in Gridlab-D with fast
convergence performance. OpenDSS also works on desktop
environments developed in the Delphi programming language,
while MAT-POWER is developed in the MATLAB
environment. However, none of these platforms can fully
exploit the capability of the existing HPC platforms, such as
Cori and Edison [17] at LBNL.
This paper investigates some parallelization strategies of
solving power flow problems in both transmission and
unbalanced distribution systems. This paper serves to 1)
introduce the open-source implementation of ExaGridPF, which
is a large-scale power flow solver designed for use on HPC
platforms, 2) explore the parallelization strategies to speed/scale
up the iterative power flow algorithms using multiple threads on
multi-core computers, 3) present the experimental results and
performance improvement over standard IEEE test cases and
synthesized larger ones. This paper is organized as follows:
Section II summarizes the formulations the power flow
algorithms for transmission and distribution systems,
respectively; Section III discusses the parallelization strategies
of iterative algorithms using existing computing architectures;
Section IV presents the experimental results and discussions;
And Section V concludes the paper with future work.

II.

&D

𝛥𝑃, = 𝑃,

POWER FLOW ALGORITHMS

Solving power flow problem is required by many power
system analysis and applications to determine voltage
magnitude and angles (in each phase) across the power grid
topology, as well as the current flow in each transmission line,
etc., given the active (𝑃 ) and reactive (𝑄 ) power drawn or
injected into phase(s) at each bus. Due to the nonlinear
properties of load modeling, i.e. ZIP load, iterative methods are
needed to solve the nonlinear problem. Numerical
approximation methods have been proposed to accelerate the
computation and guarantee the convergence. Since most
transmission systems are balanced, single-phase balanced
models are suitable for the analysis, while the current injection
method has been developed for unbalanced distribution systems
with more details. The node-voltage representation is commonly
used, which holds true for both the single-phase and three-phase
algorithms.
𝐼$%& = 𝑌$%& ⋅ 𝑉$%&

(1)

/

𝑌,, =

(2)

𝑦,.
.01

𝑌,. = 𝑌., = −𝑦,.

(3)

where 𝐼$%& and 𝑉$%& are the nodal current and voltage vectors,
respectively. 𝑌$%& denotes the bus admittance matrix. The
diagonal elements in 𝑌,, are given by equation (2), while the offdiagonal elements in 𝑌$%& are given by (3), where 𝑦,. is the line
admittance from bus 𝑖 to bus 𝑗.
A. Single-phase power flow modeling
The polar form is adopted in the single-phase modeling of
power flow algorithms, where the current for a specific bus 𝑖 can
be expressed by equation (4), where 𝜃,. is the angle different
between bus 𝑖 and bus 𝑗. 𝛿, and 𝛿. are the phase angles of bus 𝑖
and bus 𝑗, respectively. The complex form of power at bus 𝑖 is
expressed in (5); Thus, the real and imaginary parts can be
expressed separately by (6) and (7):
/

𝐼, =

𝑌,. 𝑉. ∠𝜃,. + 𝛿.

(4)

.09
/

𝑃, − 𝑗 ∗ 𝑄, = 𝑉, ∗ 𝐼, = 𝑉, ∠ − 𝛿,

𝑌,. 𝑉. ∠𝜃,. + 𝛿.

(5)

/

𝑃, =

𝑉, 𝑉. 𝑌,. cos (𝜃,. − 𝛿, + 𝛿. )

(6)

.09
/

𝑄, = −

𝑉, 𝑉. 𝑌,. sin (𝜃,. − 𝛿, + 𝛿. )

(7)

.09

Based on the above equations, the Jacobian matrix, which
represents the linear relationship between small changes of 𝑃,
and 𝑄, , i.e. Δ𝑃, and Δ𝑄, , and the voltage angle and magnitude
deviations, i.e. Δ𝛿, and Δ 𝑉, , will be constructed as follows:

(8)

&D

𝛥𝑄, = 𝑄, − 𝑄,EFG
𝐽
Δ𝑃
= 9
𝐽K
Δ𝑄

(9)

𝐽J Δ𝛿
𝐽L Δ 𝑉

(10)

Each element of the Jacobian can be obtained by taking the
first-order derivatives of 𝑃 and 𝑄 with respect to the
corresponding voltage angle and magnitude deviations. In the
runtime, both sides of the equation (10) and corresponding
elements in the Jacobian for buses with power consumption or
injection. The iterative process is illustrated as follows:
Single-phase Newton-Raphson Iteration
1
2
3
4
5
6
7

Iteration 𝑘 = 0;
Do
𝑘 = 𝑘 + 1;
Calculate power mismatch by equation (8), (9)
Update Jacobian and compute Δ𝛿,P and Δ 𝑉,P by (10);
𝛿,PQ9 = 𝛿,P + Δ𝛿,P , 𝑉,PQ9 = 𝑉,P + Δ 𝑉,P ;
While max Δ𝛿,P > 𝑒𝑟𝑟Z or max Δ 𝑉,P > 𝑒𝑟𝑟[
,∈V

,∈V

B. Three-phase power flow (Current Injection Method)
Unlike the single-phase representation of the transmission
lines, the three-phase system is more complicated with the
addition of phase information. The mutual admittance of each
branch is represented by a 3×3 matrix, taking into account the
impact of phases over each other. The details of line models can
be found in [1], [11]. For three-phase systems, we use the current
injection method to formulate the power flow problem based on
the admittance matrix construction methods (CIM) [18]
discussed in [1], [11]. The admittance matrix will contain entries
for each phase. In this paper, we use the indexing strategies in
[1], where each available phase will be indexed by the bus
number 𝑖 and the phase order 𝜙 in Φ, . For instance, the index of
phase 𝑎 on bus 𝑛 will be 3 ⋅ (𝑖 − 1) + 1. Similarly, for phase
𝑏 of bus 𝑖, the index will be 3 ⋅ (𝑖 − 1) + 2. Thus, the individual
element with bus 𝑖 and bus 𝑗 in the admittance matrix is
represented as:
FF
𝐺,.FF + 𝑗𝐵,.

F$
𝐺,.F$ + 𝑗𝐵,.

FE
𝐺,.FE + 𝑗𝐵,.

$F
𝑌,. = 𝐺,.$F + 𝑗𝐵,.

$$
𝐺,.$$ + 𝑗𝐵,.

$E
𝐺,.$E + 𝑗𝐵,.

𝐺,.E$

𝐺,.EE

𝐺,.EF

.09

− 𝑃,EFGE,

+

EF
𝑗𝐵,.

+

E$
𝑗𝐵,.

+

(11)

EE
𝑗𝐵,.

The specified active and reactive power for each bus, i.e.
&,&D
&,&D
&
𝑃,
and 𝑄, , can be obtained by subtracting the load part 𝑃f,
&
from the power output of generator 𝑃G, in (12) and (13). Note
&
that 𝑃f,
is 0 if there’s no generator attached to such bus. The ZIP
load model is given by (14) and (15).
&,&D

&
= 𝑃f,
− 𝑃G,&

&,&D

&
= 𝑄f,
− 𝑄G,&

𝑃,

𝑄,

&
&
𝑃G,& = 𝑃g,,
+ 𝑃h,, 𝑉, + 𝑃i,,
𝑉,

(12)
(13)
J

(14)

&
&
𝑄G,& = 𝑄g,,
+ 𝑄h,, 𝑉, + 𝑄i,,
𝑉,

𝛥𝐼,&
&,&D

𝐼,

=

=

&,&D
𝐼,

&,&D
𝑃,

−

J

(15)

𝐼,&,EFGE

+𝑗∗
𝑉,&

𝐼,&,EFGE =

(16)

&,&D ∗
𝑄,

=

𝑌,.&j 𝑉.j

(18)

F,,
$,,
E,,
F,,
$,,
E,, r
𝛥𝐼hof
, 𝛥𝐼hof
, 𝛥𝐼hof
, 𝛥𝐼pq
, 𝛥𝐼pq
, 𝛥𝐼pq

F,,
$,,
E,,
F,,
$,,
E,,
𝛥𝑉,E = 𝛥𝑉pq
, 𝛥𝑉pq
, 𝛥𝑉pq
, 𝛥𝑉hof
, 𝛥𝑉hof
, 𝛥𝑉hof

(19)
r

(20)

The current mismatches are calculated by equation (16),
&,&D
&,&D
where 𝑃,
and 𝑄,
are the specified active and reactive
powers for a given phase 𝑠 at bus 𝑖. Accordingly, the specified
current is given by (17). Due to the Kirchhoff's circuit laws
(KCL), the current injected into bus 𝑖 should be equal to the
summation of currents that flow out from this bus. Therefore,
the alternative way of computing the injected current is to
calculate the summation of the current flowing out into other
buses (phases) from phase 𝑠 , by equation (18). Due to the
complexity of having the relation between the power (active and
reactive) and voltage profile (magnitude and angles) in threephase systems, current-voltage relation in equation (1)
represented by the admittance matrix 𝑌, is augmented to build
the Jacobian instead. However, in the current injection method
[18], real and imaginary parts of current and voltage on each bus
are concatenated so that both LHS and RHS of the Jacobian are
real vectors. Therefore, the concatenated changes of current and
voltage vectors are given in equation (19) and (20).
Consequently, each block in Jacobian, which corresponds to
one bus, is converted from a 3×3 matrix to 6×6. The updated
block is represented by (21), where 𝐵,. and 𝐺,. are 3×3 matrice
of real and imaginary parts of 𝑌,. . Therefore, the augmented
Jacobian matrix with doubled size is given by (22).
𝐽,. =
𝐽99
Δ𝐼9
𝐽J9
Δ𝐼J
=
⋮
⋮
𝐽V9
Δ𝐼V

𝐵,.
𝐺,.
𝐽9J
𝐽JJ
⋮
𝐽VJ

𝐺,.
−𝐵,.
⋯ 𝐽9V
Δ𝑉9
⋯ 𝐽JV
Δ𝑉J
⋅
∗
⋮
𝐽PP
⋮
Δ𝑉V
⋯ 𝐽VV

∗
𝐽PP
=

(17)

.∈mn j∈kl

𝛥𝐼,E

have not yet added power injection from distributed generations,
i.e. PV buses.

(21)

(22)

However, due to the specified active and reactive power
consumption/injection on specific buses, the Jacobian is
obtained by updating the diagonal blocks, equation (23) with
adjusting elements for each bus, i.e. 𝑎, , 𝑏, , 𝑐, , 𝑑, , which are 3×3
matrice, depending on the ZIP load model. The adjusting
elements represent the first-order linear relation between small
changes of augmented current and voltage elements, given the
active 𝑃 and reactive power 𝑄 consumed at specific buses. The
approaches to generate the adjusting elements are explicitly
discussed in the appendix section of [18]. Since the primary
purpose of this paper is to explore the computational aspects, we

𝐵PP − 𝑎P
𝐺PP − 𝑐P

𝐺PP − 𝑏P
−𝐵PP − 𝑑P

(23)

During the runtime, the current mismatch is evaluated,
which is similar to the approach for the single-phase system, and
then the voltage mismatch Δ𝑉,E will be subsequently computed
by solving the linear equation, Δ𝐼 E = 𝐽 ⋅ Δ𝑉 E . The iterative
solving process stops when the maximum voltage deviation falls
below a minimal threshold. The iterative process is as follows:
Three-phase Newton-Raphson (CIM) Iterations
1
2
3
4
5
6
7
8

Iteration 𝑘 = 0;
Do
𝑘 = 𝑘 + 1;
Obtain Δ𝐼,E by (12) – (20)
Update 𝐽 and compute Δ𝑉,E by solving equation (22);
𝑉,E = 𝑉,E + Δ𝑉,E ;
Convert 𝑉,E into complex form 𝑉, ;
While max Δ𝑉,E > 𝑒𝑟𝑟[
,∈V

III.

PARALLELIZATION STRATEGIES OF EXAGRIDPF

We investigate parallelization of the power flow solver using
two techniques. Both solvers utilize the same outer Newton
iteration, while using different strategies to solve the inner
Jacobian system.
A. Newton-Raphson with SuperLU Direct Linear Solver
We leverage the SuperLU [16] sparse direct solver library to
investigate parallelizability of the power flow solver over
multiple threads. There are two stages of the solver:
factorization of the Jacobian into lower and upper triangular
matrices, and solving the two resulting triangular systems. Care
must be taken to permute the matrix to reduce the amount of fillin that results during factorization, which can dramatically
increase the number of non-zeros in the sparse matrix. We
utilize the approximate minimum degree reordering from the
COLAMD [19] package to limit the fill-in.
Since this is a direct method, there is no iteration necessary.
However, the total amount of computation may still be high due
to the inherent cost of the matrix factorization. The main benefit
of using the slower direct method is that it is more robust for illconditioned matrices, allowing the power flow to be solved in
cases where the Newton-Krylov method fails to converge.
B. Newton-Krylov with Iterative Linear Solver
For the Newton-Krylov solver, we utilize an iterative
method to solve the Jacobian system during each Newton
iteration. Specifically, we tested a bi-conjugate gradient
stabilized Krylov method [20], where each computational
routine required is parallelized over multiple threads. Each
computational thread is assigned a block of contiguous values in
the vector being updated. For example, one of the main

computations in a Krylov method is multiplication of the sparse
Jacobian with a dense vector. In the sparse matrix vector
multiplication, each thread reads a block of rows in the sparse
matrix to update its assigned block of values in the vector, which
can be calculated independently. Since no matrix factorization
is required, this method can be very efficient so long as the
matrix is well conditioned to ensure convergence of the iteration
in a small number of steps.
One drawback of this method is that it does not guarantee
convergence for ill-conditioned problems. For example, we
observed non-convergent behavior in our replicated three-phase
distribution problem. For these scenarios, the use of the direct
SuperLU solver still allowed us to solve the power flow
problem, though at increased computational cost. As with the
Newton-Raphson method, we permute the Jacobian matrix, but
the reason to do so is different since we are not factoring the
matrix. By permuting the sparse matrix so that the non-zero
elements are close to the diagonal, threads access and reuse
elements of the vector that are close to one another, increasing
the temporal and spatial data locality of access for each thread.
This enhances the data access performance of the machine by
more effectively utilizing the processor caches.
IV.

the solver time with the baseline, and show the results in Figure
2. As indicated by figure, more than 5 times speedup can be
achieved for our largest test case in the parallelized
computational routine over multiple threads.

Figure 1 Krylov solver time versus thread count

RESULTS AND DISCUSSION

We ran multiple experiments for both the transmission and
distribution systems to evaluate the computational performance
of ExaGridPF. We present some speed up metrics for running
the core Newton-Krylov solver on single-phase systems, as well
as a comparison of the execution time of the three-phase
unbalanced distribution system power flow solver utilizing the
Newton-Raphson solver with the Z-Bus method [1]
implemented in the MATLAB environment.
For the single-phase transmission system model, we
constructed large-scale cases by replicating the IEEE-30 test
case with multiple generator buses across the topologies. To
connect replicated blocks, we inserted random links to connect
the replicated transmission network blocks. For instance, the
buses of the second block are indexed from #31 to #60, and a
link to connect bus number #60 to bus #29 may be created as
there is an existing link between bus #29 and bus #30 in the
original topology. Similar approaches to create large-scale
topologies by replicating standard test cases with randomization
can be found in previous research [21]. Therefore, circles
(closed loops) can be found in the replicated topology, making
the solving task challenging for the direct methods, such as
Forward Backward Sweep or Z-Bus.
For each topology size, we evaluate the performance with
different numbers of threads on a single compute node of the
Cori computer at NERSC. The results are illustrated in Figure 1
(note the logarithmic Y-axis). The Krylov solver for the
synthesized grid topology that has 48,000 grid buses takes 1.2
seconds with 1 thread and 210 milliseconds with 32 threads.
Using the single thread performance as the baseline, we
calculate the speedup ratio for each experiment by comparing

Figure 2 Krylov solver parallel speedup
Table 1 Solving time results from ExaGridPF and Z-Bus in Matlab

Z-Bus
Bus
NO.

ExaGridPF

Win
i7/16GB

Mac
i5/8GB

Solving

Other

1.1217

1.4113

0.319107

0.047195

1.8K

3.9944

5.4393

0.389616

0.059578

3.6K

17.6397

38.928

0.799942

0.092241

0.9K

Cori (1 node)

7.2K

×

×

1.12324

0.218415

14.4K

×

×

2.31822

0.450141

28.8K

×

×

4.9977

0.954925

57.6K

×

×

11.8052

2.047031

115K

×

×

26.0988

4.346155

230K

×

×

55.291

9.980889

For the three-phase, unbalanced distribution systems, the
replication is based on IEEE 906 LV test case. We scale up the
grid topology by replicating all buses after the first bus, which
is connected to the substation transformer, whose parameters are

adjusted to prevent a voltage drop constraint violation. We
carried out a number experiments with varying number of
replicated buses to evaluate the solver performance. For
comparison, we evaluated the Z-Bus method in MATLAB on a
Windows and Mac computer. The total solver times are shown
in Figure 3. The ExaGridPF solver time increases roughly
linearly with the number of grid buses, while the solving time of
the benchmark solver (Z-Bus) exponentially increases, which
could be caused by many factors. For the largest three-phase
distribution system test case, we reduce the total solving time to
less than 1 minute, as indicated by the results in Table 1. Time
used for other miscellaneous operations, including Jacobian
initiation, Jacobian update, intermediate parameters update, etc.,
are recorded in the “Other” column in Table 1.

[2]
[3]

[4]

[5]

[6]
[7]

[8]

[9]
[10]
[11]

Figure 3 Newton-Raphson solver time versus Z-Bus method

However, as the problem size (block number) increases, the
Z-Bus solver cannot proceed due to memory shortage. × in
Table 1 indicates the lack of memory to perform the Z-Bus
calculation. In contrast, ExaGridPF scales very well, which
enables the user to perform more frequent power flow
evaluations required by many trending applications in smart grid
power systems analysis.
V.

ACKNOWLEDGEMENT

The research that leads to the results in this paper has been
supported by LDRD funding at Computational Research
Division (CRD), Lawrence Berkeley National Laboratory for
the ExaGrid project.
REFERENCES
[1]

[14]
[15]

[16]

CONCLUSION

The future work will be to further explore the strategies to
parallelize computations, especially for distribution systems
with high fidelity. The results shown in this paper are for sharedmemory computation. We are currently investigating the use of
distributed-memory computation for solving power flow
calculations over multiple nodes.
VI.

[12]
[13]

M. Bazrafshan and N. Gatsis, “Comprehensive Modeling of Three-Phase
Distribution Systems via the Bus Admittance Matrix,” IEEE Trans.
Power Syst., vol. PP, no. 99, pp. 1–1, 2017.

[17]
[18]

[19]
[20]
[21]

M. Bazrafshan and N. Gatsis, “Convergence of the Z-Bus Method for
Three-Phase Distribution Load-Flow with ZIP Loads,” IEEE Trans.
Power Syst., vol. PP, no. 99, pp. 1–1, 2017.
V. Roberge, M. Tarbouchi, and F. Okou, “Parallel Power Flow on
Graphics Processing Units for Concurrent Evaluation of Many
Networks,” IEEE Trans. Smart Grid, vol. 8, no. 4, pp. 1639–1648, Jul.
2017.
H. Nazaripouya, B. Wang, Y. Wang, P. Chu, H. R. Pota, and R. Gadh,
“Univariate time series prediction of solar power using a hybrid waveletARMA-NARX prediction method,” in 2016 IEEE/PES Transmission
and Distribution Conference and Exposition (T D), 2016, pp. 1–5.
B. Wang, Y. Wang, H. Nazaripouya, C. Qiu, C. C. Chu, and R. Gadh,
“Predictive Scheduling Framework for Electric Vehicles With
Uncertainties of User Behaviors,” IEEE Internet Things J., vol. 4, no. 1,
pp. 52–63, Feb. 2017.
B. Wang et al., “Predictive scheduling for Electric Vehicles considering
uncertainty of load and user behaviors,” in 2016 IEEE/PES Transmission
and Distribution Conference and Exposition (T D), 2016, pp. 1–5.
B. Wang, Y. Wang, C. Qiu, C. C. Chu, and R. Gadh, “Event-based
electric vehicle scheduling considering random user behaviors,” in 2015
IEEE International Conference on Smart Grid Communications
(SmartGridComm), 2015, pp. 313–318.
B. Wang, B. Hu, C. Qiu, P. Chu, and R. Gadh, “EV charging algorithm
implementation with user price preference,” in Innovative Smart Grid
Technologies Conference (ISGT), 2015 IEEE Power Energy Society,
2015, pp. 1–5.
Y. Xiong, B. Wang, C. Chu, and R. Gadh, “Distributed Optimal Vehicle
Grid Integration Strategy with User Behavior Prediction,” ArXiv Prepr.
ArXiv170304552, 2017.
Y. Wang, W. Shi, B. Wang, C.-C. Chu, and R. Gadh, “Optimal operation
of stationary and mobile batteries in distribution grids,” Appl. Energy,
vol. 190, pp. 1289–1301, Mar. 2017.
W. H. Kersting, Distribution system modeling and analysis, 3. ed. Boca
Raton, Fla.: CRC Press, 2012.
H. Saadat, Power system analysis. United States: PSA Pub., 2010.
“GridLAB-D
Simulation
Software.”
[Online].
Available:
http://www.gridlabd.org/. [Accessed: 29-Sep-2017].
“OpenDSS,”
SourceForge.
[Online].
Available:
https://sourceforge.net/projects/electricdss/. [Accessed: 29-Sep-2017].
R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas,
“MATPOWER: Steady-State Operations, Planning, and Analysis Tools
for Power Systems Research and Education,” IEEE Trans. Power Syst.,
vol. 26, no. 1, pp. 12–19, Feb. 2011.
J. Demmel, J. Gilbert, and X. Li, “An Asynchronous Parallel Supernodal
Algorithm for Sparse Gaussian Elimination,” SIAM J. Matrix Anal. Appl.,
vol. 20, no. 4, pp. 915–952, Jan. 1999.
“Computational
Systems.”
[Online].
Available:
http://www.nersc.gov/users/computational-systems/. [Accessed: 29-Sep2017].
P. A. N. Garcia, J. L. R. Pereira, S. Carneiro, V. M. da Costa, and N.
Martins, “Three-phase power flow calculations using the current
injection method,” IEEE Trans. Power Syst., vol. 15, no. 2, pp. 508–514,
May 2000.
T. A. Davis, J. R. Gilbert, S. I. Larimore, and E. G. Ng, “A Column
Approximate Minimum Degree Ordering Algorithm,” ACM Trans Math
Softw, vol. 30, no. 3, pp. 353–376, Sep. 2004.
H. van der Vorst, “Bi-CGSTAB: A Fast and Smoothly Converging
Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems,”
SIAM J. Sci. Stat. Comput., vol. 13, no. 2, pp. 631–644, Mar. 1992.
V. Jalili-Marandi, Z. Zhou, and V. Dinavahi, “Large-Scale Transient
Stability Simulation of Electrical Power Systems on Parallel GPUs,”
IEEE Trans. Parallel Distrib. Syst., vol. 23, no. 7, pp. 1255–1266, Jul.
2012.

