© 2012 IEEE. Personal use of this material is permitted.

Permission from IEEE must be obtained for all

other uses, in any current or future media, including reprinting/republishing this material for advertising or
promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or
reuse of any copyrighted component of this work in other works.
IEEE PES GENERAL MEETING, JULY 2012 1

The Holomorphic Embedding Load Flow Method

Antonio Trias

Abstract—The Holomorphic Embedding Load Flow is a novel well-known the load flow equations have multiple solutions,
general-purpose method for solving the steady state equations of and only one of them corresponds to the real operative state
power systems. Based on the techniques of Complex Analysis, of the electrical system. Unless one provides a starting point
it has been granted two US Patents. Experience has proven it
is performant and competitive with respect established iterative sufficiently close to the correct solution, iterative schemes may
methods, but its main practical features are that it is non-iterative not just fail to converge, but converge to a spurious solution.
and deterministic, yielding the correct solution when it exists Although these problems are well known, there is a certain
and, conversely, unequivocally signaling voltage collapse when it lack of awareness among practitioners. The reasons are clear:
does not. This paper reviews the embedded load flow method under normal operating conditions, convergence problems do
and highlights the technological breakthroughs that it enables:
reliable real-time applications based on unsupervised exploratory not appear often, and even the flat profile provides a good
load flows, such as Contingency Analysis, OPF, Limit-Violations initial seed. Some heuristics have been devised to come up
solvers, and Restoration plan builders. We also report on the with better starting seeds [9], [10], which helps to minimize
experience with the method in the implementation of several non-convergence cases, and efforts have been made to under-
real-time EMS products now operating at large utilities. stand and characterize the regions of convergence [11]. In any
Index Terms—Load flow analysis, power system modeling, case, non-convergence is mostly a problem in real-time, where
power system simulation, power engineering computing, energy there is no time to manually tune the seed until convergence
management, decision support systems, power system restoration.
is reached, or, more critically, to review the solution in detail
to check for possible spurious convergence at a few buses.
In the author’s opinion, it is not sufficiently recognized how
I. I NTRODUCTION these reliability problems have hindered the development of

M OST known load flow methods, particularly those used

in software for large real-world networks, are based on
general-purpose numerical iterative techniques (for a recent
real-time applications for the operators of power transmission
and distribution networks. Intelligent real-time applications
have a tremendous potential in power systems because, in
review, see [1]). The earliest methods were based on Gauss- contrast to other infrastructures under human control, there
Seidel (GS) iteration [2], which has slow convergence rates is an underlying physical model that is for the most part
but very small memory requirements. GS may still be used completely descriptive and accurate. However, real-time usage
when the other methods fail to converge starting from the requires that the physical model be solved in a fully 100%
flat profile, but methods based on Newton-Raphson (NR) [3] deterministic and reliable way. The Holomorphic Embedding
are judged to be better in general, because of their quadratic load flow method is posed to fill in this important gap.
convergence properties. However, even after applying sparse The rest of the paper is organized as follows. Section II
linear algebra techniques for efficiency [4], the full NR method describes the problems of convergence of iterative methods
is computationally too expensive due to the reevaluation of the and discusses the contexts in which such problems become
Jacobian matrix and the corresponding factorizations that need relevant. Section III exposes the Holomorphic Embedding
to be done at each iteration. Several techniques exploit the Load Flow Method (HELM), from the concepts of complex
weak coupling between active power and voltage magnitude embedding to the practical mechanics of the method. Sec-
on the one hand, and reactive power and phase angles on the tion IV exemplifies the method step by step, on the two-
other, in order to derive more efficient schemes [5], [6] where bus model. Finally, Section V briefly reports on the real-
the Jacobian matrix only needs to be factorized once. Of all world experience in the development of intelligent real-time
the various decoupled methods based on NR, the so-called applications that benefit from the full reliability of the method.
Fast Decoupled Load Flow (FDLF) formulation of Stott and
Alsac [7] has become the most successful and it is almost a
II. T HE PROBLEM OF CONVERGENCE IN ITERATIVE
de-facto standard in the industry, either in its original form or
METHODS
in one of its variants [8].
To a greater or lesser degree, all these iterative methods Numerical iterative methods are widely used in numerical
suffer from the same convergence pitfalls: on the one hand, analysis for the solution of nonlinear equations. But contrary
there is no guarantee that the iteration will always converge, to linear systems, where convergence is ensured and unique,
as this depends on the choice of the initial point; on the nonlinear systems exhibit a wide array of potential problems.
other, since the system has multiple solutions, it is not always Convergence is only ensured when the initial seed is an
possible to control to which solution it will converge. As it is estimate “sufficiently close” to the actual solution (so that the
conditions for the Banach fixed-point theorem are satisfied),
A. Trias is with Aplicaciones en Informatica Avanzada SA, Sant Cugat del
Vallès, Spain. e-mail: triast@aia.es but there are no practical criteria for quantifying how close
Manuscript received November 30, 2011. Reviewed February 8, 2012. this needs to be in general. Then systems having multiple

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 2

problem. These use the numerical techniques of homotopic

continuation, or “path-following” [19], which allow to faith-
fully track a known solution around the bifurcation point. The
problem is that, in real-time EMS applications, one commonly
needs to calculate load flow solutions for hypothetical what-if
scenarios that are not smoothly connected to any previously
known one. This happens for instance in many contingency
studies, and most notably in any application that explores
the state space by simulating SCADA actions. Therefore
using a continuation load flow for real-time applications is
Fig. 1. Fractal basins of attraction for the two-bus load flow problem, under fundamentally unfeasible, let alone slow.
the FDLF method. Initial seeds that lead to the correct solution are shown in
green; to the spurious solution in red; and non-convergence in black.
III. T HE H OLOMORPHIC E MBEDDED L OAD F LOW
M ETHOD
solutions pose additional problems, since the iterations may As evidenced in the previous section there is a need for a
converge to a non-desired solution. When different solutions direct, fully reliable load flow method. One theoretical possi-
get too close to each other, there are no general methods to bility is solving the equations exactly in closed form, using
deterministically select the desired one, unless one can exploit polynomial elimination techniques (resultants and Groebner
particular symmetries in the system. Basis) with the help of computer algebra packages [20],
[21]. This has very strong limitations, since the memory and
Moreover, for systems in complex variables (as it is the
computational costs make it impossible to get past 5 or 6
case of the load flow problem) the boundaries between the
buses. A more interesting approach from the practical point
basins of attraction under the iterative scheme are typically
of view is the so-called Series Load Flow [22], [23], based on
fractal. As two or more solutions get closer, their respective
an earlier idea by Sauer [24]. This method uses the Taylor
basins become more and more intertwined, so that, eventually,
expansion of the voltage variables as functions of all the
the neighborhood of any given solution becomes peppered
specified parameters of the problem (power injections, voltage
by points attracted to different ones. In the case of the
magnitudes), calculated on a point in which the solution is
Newton-Raphson method, this phenomenon has been studied
known. Summation of the Taylor series then allows to extend
extensively and has been given the name of Newton fractals
the solution to other scenarios, in a process that is limited by
[12]. It is a straightforward exercise to verify that Newton
the radius of convergence of the series.
fractals are indeed present in the load flow problem, using
Although it was developed completely independently, the
either the full NR or the FDLF method, as several authors
Holomorphic Embedding method presented in this paper is
have shown [13], [14]. Figure 1 shows an example using the
somewhat related to these ideas, but with one key difference:
two-bus model, where the results can be contrasted against
the Series Load Flow uses real variables and therefore cannot
the two exact solutions. The figure, which represents the
guarantee convergence of the Taylor series in general, for
complex plane for the voltage variable V , is obtained simply
arbitrary ranges. By contrast the Holomorphic Embedding
by coloring each point according to where the iteration leads
Load Flow method described in this paper is based on Com-
to. Typically, one would assign different colors for at least
plex Analysis. This seemingly minor technicality makes an
these three possible outcomes: the correct load flow solution,
enormous difference, and not just in terms of gaining new
the spurious low-voltage solution, and no-convergence. One
theoretical insights. It is only by working in the complex field
could also use colors for convergence to non-solutions (more
and using the wonderful properties of holomorphic functions,
rare, but it also happens), and for other numerical artifacts
that the method achieves its desired properties. In short, the
such as the convergence of phase angles on high multiples of
Holomorphic Embedding Load Flow provides a procedure to
2π. The author provides more examples on the web site [15].
compute, with mathematically proven guarantees of success,
For the load flow problem the consequence is that, for
the right solution to the desired accuracy (within the con-
arbitrary electrical scenarios, it is hopelessly difficult to devise
straints of the computer arithmetic accuracy); and otherwise
any a priori mechanism to select a good seed for the iterative
signals unambiguously if the system has no solution (voltage
method. Some efforts have been devoted to improve the
collapse).
seed [9], [10], but they cannot guarantee success. The problem
becomes progressively more evident as the network load
increases: when one bus approaches its maximum loadability A. Holomorphic Embedding
point (voltage collapse), a spurious solution gets closer to In the following, the main steps of the Holomorphic Em-
the real one, in what may be understood as a saddle-node bedding method are reviewed. The method is based on the
bifurcation scenario [16], [17]. Under these stressed conditions concept of embedding as general problem-solving program:
it becomes increasingly difficult for iterative methods to ensure the original problem is first submerged (embedded) in a larger
convergence to the correct solution, or any convergence at all. problem, in which the solution is easy to find under some
There is a class of methods, commonly referred to as the reference conditions. Then one expects to be able to use
Continuation Load Flow [18], which try to circumvent this such reference solution away from those conditions, the goal

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 3

being to compute the solution to the original problem. In this proposes to study the following system of algebraic equations:
case the method proposes embedding the original algebraic X sSi∗
equations in a holomorphic functional extension of them, Yik Vk (s) = sIiload +
V i (s)
which allows us to exploit the nice properties endowed by k
complex analyticity.
X sSi
Yik∗ V k (s) = sIi∗load + (3)
For the sake of clarity, the method will be exposed in Vi (s)
k
the case where all buses are of type PQ. The treatment of
where Vi (s), V i (s) are now independent complex functions
PV nodes and other types of controls is touched upon in
representing the two degrees of freedom of the voltages. Note
Section V-A. Consider then the following general form for
that when
a load flow problem:
V i (s) = Vi∗ (s∗ ) (4)
X Si∗ these equations are just complex conjugates of each other, as
Yik Vk = Iiload + (1)
Vi∗ expected. However, it should be emphasized that the converse
k
is not true: the condition (4) is not implied by the algebraic
where Yik is the generalized admittance containing branch ad- system (3), and in fact there may exist solutions which do not
mittances, bus shunt admittances, and any constant-impedance satisfy the condition and therefore are not physical solutions
injections. Symmetry in Yik is not required, so phase shifting to the original load flow (1). Therefore the embedding method
transformers are allowed. The right hand side is left with consists in solving the algebraic system (3) and requiring
the constant-injection and constant-power components of a the additional condition (4), which will be referred to as the
general ZIP load model. The proposed embedding consists in reflection condition.
introducing a complex parameter s into (1) so that the voltages From here on, the following terminology is adopted:
become functions of this new complex variable. As it will be • solutions to the algebraic embedded system which do not
shown, it is essential to the method to work in the complex satisfy (4) will be referred to as ghost solutions. They are
field. The embedding can be done in various ways, but the simply not a solution of the load flow equations.
method explicitly proposes the following form: • solutions that do satisfy the reflection condition will be
referred to as physical solutions. As it will be seen
X sSi∗
Yik Vk (s) = sIiload + (2) below, these can be either “normal” (corresponding to
Vi∗ (s∗ ) correct operating conditions, of which there may be
k
only one), or “anomalous” (corresponding to unstable
This particular embedding satisfies the first requirement of operating conditions, of which there may be several,
the method: at s = 0 all injection terms vanish and the in general). Since anomalous solutions originate in the
system is trivially solvable by linear algebra. It represents the physical state corresponding to low voltage magnitude
system under no load or generation, just swing-bus sources (low load impedance), these will be referred to as black
propagating voltage everywhere (and in the absence of shunts, solutions, and the normal solution will be referred to as
the solution is exactly |Vi | = 1, θi = 0 everywhere). As it will the white solution.
be shown below, this reference point used by the method can
Up to this point it remains to be shown that the embedding
be unambiguously defined.
in (3) does in fact define Vi (s) and V i (s) as holomorphic
Secondly, it is required for the embedding to be holomor- functions. It turns out that, since (3) are algebraic, elimination
phic, that is, it should define the voltages Vi to be holomorphic techniques based on the theory of resultants and Gröbner
functions in the embedding parameter s. This is done in Basis [25] guarantee that all variables Vi can be successively
order to benefit from all the power of complex analysis, in eliminated in terms of the remaining ones, until a polynomial
particular the process of analytical continuation that will allow equation in V1 is obtained:
to obtain the objective state from the reference one. It should
N
be strongly remarked that the denominator on the right hand X
side of (2) has the form V ∗ (s∗ ), and not V ∗ (s). It can be P(V1 ) = pn (s)V1n = 0 (5)
n=0
shown that, because complex conjugation does not leave the
Cauchy-Riemann equations invariant, this is the only choice and the V 1 , V2 , V 2 , V3 , etc., are expressed explicitly as poly-
that allows V (s) to have a chance of being holomorphic. It is nomials in all the previous ones in a triangular manner, which
then useful to define V (s) ≡ V ∗ (s∗ ), so that the embedded allows the obtention of all other Vi and V i from each solution
system becomes: V1 , by simple progressive back-substitution. The degree N of
this polynomial is in general rather large (of order exponential
X sSi∗ in the number of variables Vi in the original system), but
Yik Vk (s) = sIiload +
k
V i (s) always finite, and the coefficients pn (s) are polynomial in s.
This is precisely the definition of an algebraic curve. There-
However this change of notation does not hide the fact that fore all Vi and V i are proved to be holomorphic functions
is cumbersome to proceed with the analysis of the embedded everywhere except on a finite number of points, known as
equations by requiring V i (s) = Vi∗ (s∗ ) throughout the treat- the exceptional set of the algebraic curve. These exceptional
ment. In order to escape this difficulty, the embedding method points are those values of s on which the polynomial equation

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 4

for V1 exhibits a null derivative ∂P/∂V1 , and will play an formal power series appears on both sides:
important role on the discussion about analytic continuation X ∞
X ∞
X
further below. Yik ck [n]sn = sIiload + sSi∗ d∗i [n]sn (6)
k n=0 n=0

It is now straightforward to realize that (6) provides a way

B. Power series to progressively obtain the coefficients order by order, up to
Obtaining and solving (5) for large networks is unfeasible in any desired level. The procedure is bootstrapped by making
practice. Instead, one should work with the power series of the s = 0 in (6), which yields the linear system:
holomorphic functions, as is commonly done when calculating X
Yik ck [0] = 0
algebraic curves. The power series provides the so-called germ
k
of the analytical function, as it allows to calculate the function
well away from the convergence radius of the series, by the from which the zero-th order coefficients are obtained. Taking
well-known mechanism of analytic continuation. the derivative of (6) with respect to s and again evaluating at
As seen above, the algebraic problem has multiple solutions. s = 0: X
Indeed at s = 0, multiple solutions as can be found for (3) if Yik ck [1] = Iiload + Si∗ d∗i [0] (7)
k
one allows Vi (0) = 0 and V i (0) = 0 at one or more buses i.
However, those solutions are either ghost (unphysical) or black But the coefficients dk [n] can be obtained from knowledge
(non-vanishing injection at s = 0, meaning a short-circuit at of the coefficients ck [m] (up to m = n) and dk [m] (up to
the bus). On the other hand, by requiring V i (0) 6= 0 and m = n−1), since they are related by the convolution formulas:
V i (0) 6= 0 for all buses i, all injection terms in (3) vanish at X∞
! ∞
X
!
s = 0 and the resulting linear system has a unique solution, −1 n n
1 = V (s)V (s) = c[n]s d[n]s
dubbed the white solution. This is the operational solution n=0 n=0
of a power system under a no-load, no-generation scenario, 1
X 2
X
where only the swing bus is providing a voltage source. The = c0 d0 + s c[1 − m]d[m] + s2 c[2 − m]d[m]
HELM method proposes that the operational solution to the m=0 m=0
n
load flow problem is the maximal analytical continuation of n
X
+··· + s c[n − m]d[m] + . . .
this reference solution at s = 0.
m=0
An interesting question may arise now, as to what solution
is actually realized in a real power system. Is it necessarily the Thus dk [0] = 1/ck [0] and the system (7) is readily solved to
continuation of the white solution, as proposed by the method? obtain the coefficients ck [1]. The same procedure can now be
To answer this question satisfactorily we need to remind repeated order after order:
ourselves that the load flow equations are just describing
X
Yik ck [n] = Si∗ d∗i [n − 1] (8)
the steady-state of the physical system. To be certain about k
the actual state of the system we would need to integrate
and by the convolution formulas:
the differential equations of the full dynamical model, using
Pn−1
some initial conditions. This would yield a unique steady ci [n − m]di [m]
state because there is no ambiguity in the solution of those di [n] = − m=0 (9)
ci [0]
equations. This is almost impossible to achieve in practice, but
so that at every step a linear system is solved. Therefore the
it is certain that the steady state will be analytic in the problem
power series can be computed by a univocal and well-defined
parameters, given some minimal assumptions of smoothness
procedure which essentially consists in solving linear systems.
for the differential equations. In the absence of this modeling,
In terms of computational work, it should be remarked that
solving the algebraic equations for the steady state actually
the matrix remains constant, and therefore its factorization
requires to make a proposal as to what solution is the correct
only needs to be done once. Efficient implementations of the
one for the physical systyem. In this sense, and invoking
method should of course make use of modern sparse linear
the smoothness of the underlying dynamical model, it seems
algebra routines [26], which include reordering algorithms for
most reasonable to propose the solution that is the analytical
the minimization of the fill-up in the matrix factors. Sec-
continuation of the white solution at the reference limit s = 0.
tion III-E analyzes the performance characteristics in further
Note however that the HELM procedure can also be applied
detail.
to the calculation of black branches, and it possibly provides
the best framework in which this can be done in a systematic
manner. C. Analytic continuation
The power series for the white P solution can now be calcu- Since the procedure shown above can be carried out to
∞
lated as follows. Assume Vi (s) = n=0P ci [n]sn are the power arbitrary orders, a natural question arises as to how many
∞
series to be calculated, and 1/Vi (s) = n=0 di [n]sn are the orders it is needed to calculate, in order to obtain the final
corresponding power series for the functions appearing on the solution at the objective state s = 1. It should be first
right hand side of (3). Making use of the reflection condition noted that the solution will not in general be obtained by
V i (s) = Vi∗ (s∗ ), the following equation is obtained where a direct summation of the power series at s = 1, as the

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 5

radius of convergence is typically much smaller than 1. The 1) Stahl’s extremal domain theorem [29]–[31]: this result
powerful procedure of analytic continuation is used instead asserts that for any analytic function there exists a
[27]. In passing, it should be stressed that the process of unique set of cuts with the property of having minimal
analytic continuation does not have anything in common with logarithmic capacity [32], and such that the function
the concepts of numerical continuation (homotopy methods has single-valued analytic continuation in the domain
[19]) used in continuation load flow methods. In practice the consisting of the complex plane excluding the cuts (i.e.,
analytical continuation is carried out by means of rational the maximal domain). This provides a natural criteria for
approximants, among which Padé approximation is the method the choice of cuts, and ensures that such choice exists
of choice for reasons to be revealed shortly. As to the question and is unique.
of the number of terms needed in the power series in order 2) Stahl’s Padé convergence theorem [33], [34]: for any
to attain a given level of precision, practice has shown that analytic function whose singularities are finite (in fact,
typically anywhere from 10 to 40 terms suffice to reach 5-digit it suffices for the set of singularities to have zero
precision in large networks, and about 60 terms will exhaust logarithmic capacity), any close-to-diagonal sequence of
the limits of the computer arithmetic in double precision. Padé approximants converge in capacity to said function
The mechanics of the method have now been completely in the extremal domain. The poles of the diagonal and
described, but another important issue needs to be addressed: paradiagonal Padé approximants accumulate on the set
is the analytical continuation procedure “complete”? That is, of cuts with minimal logarithmic capacity.
does it always reach the solution when it exists? Conversely, Padé approximants [31] are rational approximants to power
will the procedure unambiguously signal non-existence when series, and they have been used extensively as a technique
the solution does not exist?. To answer these questions it for analytic continuation because their convergence has been
is needed to invoke some powerful results from Complex known to be much better than that of power series. In
Analysis. particular, the diagonal and paradiagonal Padé approximants
As it is well-known from the theory of Algebraic Curves, coincide with the continued fraction approximation to the
all solutions Vi (s) and V i (s) are holomorphic functions that power series, which are also known to have good convergence
can be analytically continued along any path in s, as long properties in general. Stahl’s results, after the seminal works of
as this path does not contain points of the exceptional set of Nuttall [35], reveal that Padé approximants are really a means
the curve. This result can be alternatively stated by saying for maximal analytic continuation. Therefore these two results
that Algebraic Curves are analytical functions whose only confer the method very strong additional guarantees: if the
singularities are branch points, since the exceptional points Padé approximants converge at s = 1, the result is guaranteed
of the curve (the points in s where the zeros of (5) have to be the analytic continuation of the white branch at s = 1;
multiplicity greater than one, or equivalently, the points on the conversely, if the Padé approximants do not converge at s = 1
curve where ∂P/∂V1 = 0) are in fact branch points. These (i.e. the point s = 1 lies on the set of cuts with minimal
are the points where two or more branches of the curve V1 (s) logarithmic capacity) then it is guaranteed that there is no
coalesce. The algebraic curve as a whole (all its branches) solution (that is, the system is beyond voltage collapse).
is a complete global analytic function [28], which means
that there exist paths of analytic continuation connecting any
points between any branches. In other words, knowledge of the D. The method in brief
power series at any (non-exceptional) point can be exploited
to calculate the solution anywhere else, on any branch, by To recapitulate, the Holomorphic Embedding Load Flow
means of a suitable path of analytic continuation. However, boils down to these steps:
note that analytical continuation paths enclosing branch points 1) Choose a suitable complex embedding by means of
yield the curve on a branch different from the starting one. This a complex parameter s. The embedding needs to be
is the subject of Monodromy Theory and is a foundational holomorphic (uses V ∗ (s∗ ), not V ∗ (s)). At s = 0, this
part of the theory of Riemann surfaces, which is the natural embedding should be such that the system becomes
setting in which to study multivalued complex functions (and linear and trivially simple to solve (the no-load, no-
in particular, algebraic curves). For the purposes of the load generation case). This unambiguously selects the ref-
flow method, the aim is to perform analytic continuation of erence solution at s = 0.
the reference solution at s = 0 along paths that ensure single- 2) Calculate the power series of V (s) corresponding to the
valuedness, in other words, remaining always within the white reference solution, by means of a sequence of linear
branch. This is accomplished by the well-known procedure systems that yield the coefficients progressively, order
of selecting branch cuts on the complex plane. Branch cuts after order. The matrix in those systems remains always
consist of lines connecting all branch points in such a way constant, so it needs to be factorized just once; and
that no analytic continuation paths can encircle isolated branch the right hand sides can always be calculated from the
points, in which case Monodromy Theory ensures that any results of the previous system.
chosen path for analytical continuation is single-valued. 3) Compute the solution at s = 1 as the analytical continua-
However, since the specific geometry of cuts is arbitrary, tion of the power series obtained in step 2, by using Padé
some criteria need to be defined in order to choose them. At Approximants. These are guaranteed to yield maximal
this point, two key results are invoked: analytical continuation, therefore the solution is obtained

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 6

when it exists, or a divergence is obtained when it does convenient to introduce adimensional variables by making
ZS ∗
not exist. U ≡ VV0 and σ ≡ |V 0|
2 , so that the equation becomes

Steps 2 and 3 are typically performed in an interleaved fashion, σ

that is, after a new coefficient of the series is obtained, a U =1+ (11)
U∗
new Padé approximant can be computed. The procedure stops
which exhibits the essential algebraic structure of the load
when the desired accuracy is obtained, or when oscillation
flow problem in its purest form. The exact solution is readily
or divergence is detected. Section IV below will illustrate the
obtained by simple algebra manipulations:
whole procedure step by step. r
1 1
E. Performance U= ± + σR − σI2 + jσI (12)
2 4
The performance characteristics of the method are fairly subject to the condition ∆ ≡ 14 + σR − σI2 ≥ 0, where σR , σI
easy to analyze in the asymptotic limit, when the size N of are the real and imaginary parts of σ, respectively.
the network is large. Most of the work consists in factorizing As summarized in Section III-D above, the HELM al-
a matrix system (8), and using the result to repeatedly solve gorithm proceeds as follows. First, a suitable holomorphic
different right-hand sides. Modern direct sparse solvers [26] embedding is proposed:
are able to factorize a matrix A into its LDLT Cholesky
sσ
decomposition with a cost of O(|L|), that is, linear in the U (s) = 1 + ∗ ∗ (13)
U (s )
number of non-zero entries obtained in the factorized matrix
L. Therefore it is of utmost importance to use a fill-reducing The second step consists
P in nsolving for the coefficients of the
reordering of the matrix. Solving each system adds O(N ) power series U (s) = cnP s . For this we need the coefficients
operations due to forward and back substitution. Computing of the function 1/U (s) = dn sn . The system can be solved
the right hand sides of each linear system involves computing order after order, yielding the solution:
a simple convolution formula, reaching at most the order of
c0 = 1
the coefficient we are solving in the series. Since in practice
there is a maximum of about 60 coefficients to solve before cn+1 = σd∗n (n = 0, . . . , ∞) (14)
one encounters the limits of double-precision accuracy, the where the coefficients dn can always be obtained from those
cost of calculation of these convolution formulas can consid- previously
ered to be fixed for each bus and thus contributes another Pcalculated,
n−1
through the simple convolution formula
dn = − k=0 cn−k dk . Here is the result for the first coeffi-
O(N ). The same can be said about the computation of Padé cients up to n = 4:
Approximants. The cost of computing approximants increases
approximately as the square of the number of coefficients c1 = σ c3 = σ 2 σ ∗ + σσ ∗2
being used, but since this number tops off at about a maximum c2 = −σσ ∗ c4 = −σ 3 σ ∗ − 3σ 2 σ ∗2 − σσ ∗3
of 60, the costs are still O(N ) asymptotically. Moreover, The third step consists in the construction of the Padé Approx-
both the computation of the right-hand sides and the Padé imants using the power series coefficients. This is a standard
approximants can be trivially parallellized. Therefore for large procedure for which several efficient methods are available,
enough networks the costs are dominated by the factorization such as Wynn’s  algorithm or the QD algorithm [31]. In
of the impedance matrix. These costs are then similar to this case it is instructive to use the general method, which
the FDNR method, since the matrix factorization is only involves solving a linear systems but provides the coefficients
performed once. of the Padé approximants, instead of just its value at s = 1.
To give some practical figures, in a real-world large trans- Usually it is the diagonal and near-diagonal sequences of Padé
mission network of about 3,000 electrical nodes the HELM approximants that give the best precision. Using the traditional
algorithm solves a load flow in about 10 to 20 ms, running [L/M ] notation, the first n coefficients of the series allow us
on Intel Xeon 5500 CPUs, for a required precision of about to obtain approximants up to L + M = n:
5 significant digits. The maximum orders for the Padé ap-
proximants are typically about 20, going up to a maximum of 1 + (σ + σ ∗ )s
[1/1] =
about 60 in scenarios where the network approaches voltage 1 + σ∗ s
collapse. 1 + (2σ + σ ∗ )s + σ 2 s2
[2/1] =
1 + (σ + σ ∗ )s
IV. A WORKED EXAMPLE : THE TWO - BUS SYSTEM 1 + (2σ + 2σ ∗ )s + (σ 2 + σσ ∗ + σ ∗2 )s2
[2/2] =
The two-bus model offers an excelent playground to exem- 1 + (σ + 2σ ∗ )s + σ ∗2 s2
plify the HELM method and test the results against the exact
It is now straightforward to evaluate these and compare them
known solution. Using the same sign convention as in (1), the
against the exact solution. Table I shows the precision obtained
equation for the system is:
with some numerically computed Padé approximants, for three
V − V0 S∗ example scenarios. The third column corresponds to a case on
= ∗ (10)
Z V the verge of collapse (|V | = 0.58).
where V0 represents the swing. In this case the linear algebra It is worthwhile mentioning that the two-bus model explic-
will be trivial since the system is one-dimensional. It is itly exposes further insights into the problem. Looking at (13)

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 7

TABLE I
E XAMPLES OF RELATIVE PRECISION OF PAD É A PPROXIMANTS VS . ORDER applications. The HELM method was actually born out of the
need for a fully reliable load flow in the context of some
Padé order σ = −0.07 − j0.08 −0.14 − j0.15 −0.2 − j0.22 AI-based applications that depend critically on the ability
(|V | = 0.92) (|V | = 0.81) (|V | = 0.58) to perform exploratory load flow studies, with absolutely
[2/2] 1.795e-03 2.30e-02 2.74e-01 no margin for failure. The most prominent examples are
[5/5] 7.02e-09 2.10e-05 5.85e-02 two decision-support tools, a Limits Violation Solver and a
[10/10] 0 1.89e-10 1.16e-02
[15/15] 0 1.55e-15 2.86e-03 Restoration Plan Builder. These tools are fully model-based
[20/20] 0 0 7.38e-04 thanks to a technique well-known in the AI community: guided
exploration in the state-space of the electrical system, using
the A∗ algorithm. The state-space consists of all possible
and its complex comjugate, it is realized that the equations are electrical (steady) states that the network can achieve, and the
telling us explicitely what the solution is, in continued fraction available SCADA actions provide transitions between them.
form. Just subsitute the denominators iteratively: The algorithm needs sophisticated heuristics to guide the
search efficiently, but the load flow method needs to be 100%
σs reliable, as it is used at each and every step of the exploration.
U (s) = 1 + (15)
σ∗ s Our experience has shown that these kind of tools would be
1+
σs impossible to build on top of iterative load flow methods. Of
1+ course, other real-time tools such as Contingency Analysis or
σ∗ s
1+ PV/QV Curves also benefit from increased reliability.
1 + ···
Another promising yet untapped potential of the method
This is another representation of the function, but continued lies in the new insights it brings into the analysis of the load
fractions have in general a much larger convergence radius flow problem. The treatment in terms algebraic curves could
and faster convergence compared to power series. In fact, it prove quite powerful. For instance, it provides a coherent
is well known that the convergents (i.e. truncations) of this framework for the characterization and computation of all
continued fraction coincide with the diagonal and paradiagonal the multiple solutions to the original problem (white, black,
Padé approximants [31]. This may be verified against the ghost solutions). Given the vast amounts of results in the
approximants computed above in (15). field of algebraic curves in Complex Analysis, it is reasonable
The general N -bus case does not easily lend itself into to think that this is just scratching the surface of what is
the elegant continued fraction approach, since there are linear potentially possible. The theory of approximants (rational or
systems involved. However, given the equivalence between other) is another source for insights and practical results.
Padé approximants and continued fractions, the analysis of As it has been shown, the zeros and poles of the rational
the two-bus model is quite relevant, as the essential algebraic approximants tend to accumulate on the (minimal) branch cuts
structure of the problem is already there. of the functions V (s). Therefore their values, or even their
patterns of appearance as the approximant order increases, may
V. F INAL R EMARKS be used as new indicators, such as the proximity to voltage
A. Controls colapse. One may think of all this as a sort of a new language
Thus far only the pure PQ case has been considered. Real for the analysis of an old problem.
power systems have all sorts of automated controls imposing
additional constraints on the solution of (1), the most pervasive VI. C ONCLUSIONS
being the PV controls for generator buses. These and other
This paper has presented a novel load flow method that rad-
controls such as transformer ULTCs, FACTS, o HVDC links
ically breaks away from the established iterative methods. Its
are easily accommodated under the HELM methodology.
most salient features are that it is non-iterative, deterministic,
The explicit procedure for each type of constraint, although
and non-ambiguous: it guarantees, backed by mathematical
straightforward, is lengthy and will be the subject of a follow
proof, obtaining the right solution to the multivalued load
on paper [36]. Here it is pointed out that the only requirement
flow problem, and otherwise signals unambiguously the non-
is that the constraints can be expressed as algebraic equalities,
existence of solution when the system is beyond voltage
which is true in all the aforementioned cases. Since no
collapse. The method is based on a holomorphic embedding
approximations are needed, the solutions benefit from all the
procedure that extends the voltage variables into analytic
mathematical guarantees that have been shown in this paper.
functions in the complex plane. This provides a framework to
This is in contrast to iterative methods, where controls are
study and obtain the solutions using the full power of complex
introduced as adjustments at each iteration. Such adjustments
analysis. The method provides a procedure for constructing the
do not lend themselves to rigorous analysis and their effects
complex power series at a well-defined reference point, where
on the solutions have to be studied empirically [37], [38].
it is trivial to identify the correct branch of the multivalued
problem, and then uses analytical continuation by means of
B. Significance of HELM algebraic approximants to reach the objective. It can be proven
Arguably the single most important impact of the HELM that the continuation is maximal in logarithmic capacity,
algorithm is the enabling of reliable, real-time, intelligent thus propagating the chosen branch to the maximal possible

© 2012 IEEE
IEEE PES GENERAL MEETING, JULY 2012 8

domain on the complex plane. If the objective point is not [14] J. Thorp and S. Naqavi, “Load-flow fractals draw clues to erratic
in this domain, the initial correct branch does not have an behaviour,” IEEE Comput. Appl. Power, vol. 10, no. 1, pp. 59–62, jan
1997.
analytical continuation there and the problem has no solution. [15] Grupo AIA, “Problems with iterative load flow,” http://www.elequant.
Therefore, by construction, the process is totally deterministic com/products/agora/demo/iterativeloadflow, 2002–11.
thus ensuring that if there is a solution the method will find it [16] I. Dobson, “Observations on the geometry of saddle node bifurcation
and voltage collapse in electrical power systems,” IEEE Trans. Circuits
and, conversely, if there is no solution (voltage collapse) the Syst. I, vol. 39, no. 3, pp. 240–243, mar 1992.
method will unequivocally signal such condition as well. [17] C. Cañizares, “On bifurcations, voltage collapse and load modeling,”
Although advanced results from geometric function theory IEEE Trans. Power Syst., vol. 10, no. 1, pp. 512–522, feb 1995.
[18] V. Ajjarapu and C. Christy, “The continuation power flow: a tool for
are used to prove the method properties, the numerical im- steady state voltage stability analysis,” IEEE Trans. Power Syst., vol. 7,
plementation is straightforward and has performance charac- no. 1, pp. 416–423, feb 1992.
teristics that make it competitive and even superior to fast- [19] E. Allgower and K. Georg, Introduction to numerical continuation
methods, ser. Classics in applied mathematics. SIAM, 2003.
decoupled algorithms, thus making it a good general-purpose [20] A. Montes, “Algebraic solution of the load-flow problem for a 4-nodes
load flow method for power systems of any size. The method electrical network,” Mathematics and Computers in Simulation, vol. 45,
has been implemented in industrial-strength EMS applications no. 1-2, pp. 163–174, 1998.
[21] J. Ning, W. Gao, G. Radman, and J. Liu, “The application of the
now operating at several large transmission operators, and has groebner basis technique in power flow study,” in North American Power
been granted two US Patents [39]. Experience has proven the Symposium (NAPS), 2009, oct. 2009, pp. 1–7.
foremost practical impact of the method, namely the enabling [22] W. Xu, Y. Liu, J. Salmon, T. Le, and G. Chang, “Series load flow:
a novel noniterative load flow method,” Generation, Transmission and
of reliable real-time applications. Distribution, IEE Proceedings, vol. 145, no. 3, pp. 251–256, may 1998.
[23] A. Zambroni De Souza, C. Rosa, B. Lima Lopes, R. Leme, and
O. Carpinteiro, “Non-iterative load-flow method as a tool for voltage
ACKNOWLEDGMENT stability studies,” Generation, Transmission Distribution, IET, vol. 1,
no. 3, pp. 499–505, may 2007.
The author would like to thank J. A. Marqués and V. Gaitán [24] P. Sauer, “Explicit load flow series and functions,” IEEE Trans. Power
for their numerical implementations and extensive testing of App. Syst., vol. PAS-100, no. 8, pp. 3754–3763, aug. 1981.
the HELM algorithms since its inception, and J. L. Marı́n [25] B. Sturmfels, Solving systems of polynomial equations, ser. Regional
conference series in mathematics. AMS, 2002.
for his help in the preparation of this manuscript. The author [26] T. Davis, Direct methods for sparse linear systems, ser. Fundamentals
also thanks the reviewers of the manuscript, whose comments of algorithms. SIAM, 2006.
helped improve the clarity of the paper. [27] T. Gamelin, Complex analysis. Springer, 2001.
[28] L. Ahlfors, Complex analysis: an introduction to the theory of analytic
functions of one complex variable. McGraw-Hill, 1979.
R EFERENCES [29] H. Stahl, “Extremal domains associated with an analytic function I, II,”
Complex Variables, Theory and Application: An International Journal,
[1] A. Gómez-Expósito and F. L. Alvarado, “Load flow,” in Electric Energy vol. 4, no. 4, pp. 311–324, 325–338, 1985.
Systems: Analysis and Operation, A. Gómez-Expósito, A. J. Conejo, and [30] ——, “The structure of extremal domains associated with an analytic
C. Cañizares, Eds. CRC Press, 2009, pp. 95–126. function,” Complex Variables, Theory and Application: An International
[2] J. B. Ward and H. W. Hale, “Digital computer solution of power-flow Journal, vol. 4, no. 4, pp. 339–354, 1985.
problems,” Power Apparatus and Systems, Part III. Transactions of the [31] G. Baker and P. Graves-Morris, Padé approximants, ser. Encyclopedia
American Institute of Electrical Engineers, vol. 75, no. 3, pp. 398–404, of mathematics and its applications. Cambridge University Press, 1996.
jan. 1956. [32] R. Rumely, Capacity theory on algebraic curves, ser. Lecture notes in
[3] W. Tinney and C. Hart, “Power flow solution by newton’s method,” mathematics. Springer-Verlag, 1989.
IEEE Trans. Power App. Syst., vol. PAS-86, no. 11, pp. 1449–1460, [33] H. Stahl, “On the convergence of generalized Padé approximants,”
nov. 1967. Constructive Approximation, vol. 5, pp. 221–240, 1989.
[4] W. Tinney and J. Walker, “Direct solutions of sparse network equations [34] ——, “The convergence of Padé approximants to functions with branch
by optimally ordered triangular factorization,” Proc. IEEE, vol. 55, points,” Journal of Approximation Theory, vol. 91, no. 2, pp. 139–204,
no. 11, pp. 1801–1809, nov. 1967. 1997.
[5] S. Despotovic, B. Babic, and V. Mastilovic, “A rapid and reliable method [35] J. Nuttall, “On convergence of Padé approximants to functions with
for solving load flow problems,” IEEE Trans. Power App. Syst., vol. branch poits,” in Padé and Rational Approximation, E. B. Saff and R. S.
PAS-90, no. 1, pp. 123–130, jan. 1971. Varga, Eds. Academic Press, New York, 1977, pp. 101–109.
[6] B. Stott, “Decoupled newton load flow,” IEEE Trans. Power App. Syst., [36] A. Trias, “Holomorphic embedding load flow: Controls,” 2012, preprint.
vol. PAS-91, no. 5, pp. 1955–1959, sept. 1972. [37] S.-K. Chang and V. Brandwajn, “Solving the adjustment interactions in
[7] B. Stott and O. Alsac, “Fast decoupled load flow,” IEEE Trans. Power fast decoupled load flow,” IEEE Trans. Power Syst., vol. 6, no. 2, pp.
App. Syst., vol. PAS-93, no. 3, pp. 859–869, may 1974. 801–805, may 1991.
[8] R. van Amerongen, “A general-purpose version of the fast decoupled [38] J. Vlachogiannis, “Control adjustments in fast decoupled load flow,”
load flow,” IEEE Trans. Power Syst., vol. 4, no. 2, pp. 760–770, may Electric Power Systems Research, vol. 31, no. 3, pp. 185–194, 1994.
1989. [39] A. Trias, “System and method for monitoring and managing electrical
[9] G. Leonidopoulos, “Approximate linear decoupled solution as the initial power transmission and distribution networks,” US Patents 7 519 506
value of power system load flow,” Electric Power Systems Research, and 7 979 239, 2009–2011.
vol. 32, no. 3, pp. 161–163, 1995.
[10] R. Klump and T. Overbye, “Techniques for improving power flow
convergence,” in Power Engineering Society Summer Meeting, 2000.
IEEE, vol. 1, 2000, pp. 598–603.
[11] F. Wu, “Theoretical study of the convergence of the fast decoupled load
flow,” IEEE Trans. Power App. Syst., vol. 96, no. 1, pp. 268–275, jan.
1977. Antonio Trias received a PhD degree in Physics from the University of
[12] J. Hubbard, D. Schleicher, and S. Sutherland, “How to find all roots of Barcelona in 1974. From 1974 to 1976 he was a visiting postdoctoral research
complex polynomials by newton’s method,” Inventiones Mathematicae, fellow at the Lawrence Berkeley Laboratory of the University of California.
vol. 146, pp. 1–33, 2001. He is a founding member of Aplicaciones en Informática Avanzada S.A.,
[13] R. Klump and T. Overbye, “A new method for finding low-voltage power where he currently is R&D Vice President and CTO. He has worked in the
flow solutions,” in Power Engineering Society Summer Meeting, 2000. load flow problem and other algorithms of interest for power systems for over
IEEE, vol. 1, 2000, pp. 593–597. twenty years.

© 2012 IEEE