IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

1

Convergence of the Z-Bus Method for Three-Phase
Distribution Load-Flow with ZIP Loads

arXiv:1605.08511v3 [math.OC] 9 May 2017

Mohammadhafez Bazrafshan, Student Member, IEEE, and Nikolaos Gatsis, Member, IEEE

Abstract—This paper derives a set of sufficient conditions
guaranteeing that the load-flow problem in unbalanced threephase distribution networks with wye and delta ZIP loads has
a unique solution over a region that can be explicitly calculated
from the network parameters. It is also proved that the wellknown Z-Bus iterative method is a contraction over the defined
region, and hence converges to the unique solution.
Index Terms—Three-phase distribution load-flow, Z-Bus
method, ZIP load, contraction mapping, power flow

I. I NTRODUCTION

E

XISTENCE of three-, two-, and one-phase transmission
lines with high R/X ratios in distribution networks is
predominant. Therefore, the Z-Bus iterative method is preferred over classical Newton-Raphson methods for distribution
load-flow [1]. This paper deals with the Z-Bus method when
applied to practical distribution networks including three-phase
wye and delta loads with constant-power, constant-current, and
constant-impedance portions (ZIP loads), as well as potentially
untransposed lines. Conditions on the network components and
loads guaranteeing the existence of a unique load-flow solution
and the convergence of the Z-Bus method are derived.
The availability of conditions guaranteeing power flow
solution existence has several applications in power system
planning and operations. In planning studies, such as conventional and distributed generation planning, conditions for
solution existence serve as assurance that desirable operating conditions are feasible under future load demands. For
example, such conditions provide network-specific guides on
whether a certain lateral will have valid load-flow solutions
over the course of years, as the loads increase [2], [3].
In operations, conditions for solution existence allow power
engineers to determine whether the set of power injections and
flows respect security constraints of power system equipment
[4]. Furthermore, the availability of solvability regions in
terms of nodal voltages renders improved initial estimates
for numerical algorithms and accelerates the online analysis
required for real-time monitoring and control applications [5].
Specifically in regards to distribution networks, their
weakly-meshed structure limits the number of feasible voltage
profiles. This characteristic implies that it is harder to meet
different demand patterns in distribution networks [6]. In this
case, explicit sufficient conditions for solution existence that
Manuscript received May 26, 2016; revised November 22, 2016 and March
1, 2017; accepted April 26, 2017. The authors are with the Department of
Electrical and Computer Engineering, The University of Texas at San Antonio. E-mails: {mohammadhafez.bazrafshan, nikolaos.gatsis}@utsa.edu. This
material is based upon work supported by the National Science Foundation
under Grant No. CCF-1421583.

take into account various load types—that is, wye and delta
ZIP loads—prove to be very useful.
The literature on existence of solutions to the power flow
equations is extensive; see e.g., [3] for a review on general
power networks, and recent works in [5], [7], [8] geared toward
distribution networks. Here, the prior art most closely related
to the present work is reviewed.
The single-phase equivalent of a power distribution network
is considered in [5]. Power flow equations are formulated so
that injection currents can be expressed as a fixed point of a
quadratic map parametrized by constant-power loads and the
network admittance matrix. Sufficient conditions on the loads
and the network admittance matrix are derived to guarantee
existence of a power flow solution.
The work in [7] highlights that the injection currents can be
expressed by a family of fixed-point quadratic maps similarly
to [5] but parametrized by a diagonal matrix. An algorithm
is proposed for finding the diagonal matrix that expands the
regions given in [5] where a power flow solution exists.
Convergence of the Z-Bus method interpreted as a fixedpoint iteration is considered in [8] and [9] for singlephase distribution networks. Specifically, explicit conditions
on constant-power loads and the network admittance matrix
are presented in [8] ensuring convergence of the Z-Bus method
and existence of a unique solution for distribution load-flow
under balanced assumptions. Interestingly, it is numerically
verified in [9] that when distributed generation units are modeled as constant-power nodes, the voltage iterates produced by
the Z-Bus method are indeed contracting.
Single-phase power system modeling is valid under balanced three-phase operation, which is a reasonably accurate assumption for transmission networks. In distribution networks,
there is a prevalence of untransposed three-, two- and onephase transmission lines with unbalanced loading. There is
a gap in the literature, since theoretical results on solution
existence and convergence of load-flow methods overlook the
underlying unbalanced nature of distribution systems.
This paper considers general unbalanced distribution networks that include transmission lines with one, two, or three
available phases as well as wye and delta ZIP loads, in contrast
to [5], [7], [8], which only consider single-phase networks
with constant-power loads. Sufficient conditions for the Z-Bus
method to be a contraction over a region explicitly defined in
terms of the loads and the bus admittance matrix are derived.
A consequence of the contraction is that power-flow equations
are guaranteed to have a unique solution in the defined region.
In order to prove that the Z-Bus method is a contraction
mapping, a candidate region of contraction is defined and
the self-mapping and contraction properties are proved. These

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

steps were followed in [8] for single-phase networks with
constant-power loads. However, the extension to general unbalanced three-phase networks with wye and delta ZIP loads
and missing phases adds a level of complexity that necessitates
the development of novel methods to prove the convergence
of the Z-Bus method.
The major contribution of this paper is to derive sufficient
convergence conditions that handle the complexities arising
from various combinations of load types and phasings in
practical three-phase distribution networks. The conditions
are expressed directly in terms of the network loads and
admittance matrix, and guarantee the existence of a unique
solution to the power flow equations. It becomes immediately
apparent from these conditions that different load types contribute differently to the network solvability. Moreover, the
generality of the network model renders results that are readily
applicable to practical distribution networks without simplifying assumptions. To this end, the conditions are validated on
several IEEE distribution test feeders.
In order to prove the aforementioned conditions, this paper
develops a set of tools that advance the power flow analysis
via contraction methods. These tools are particularly tailored
to handle wye and delta connections, ZIP loads, and missing
phases, and are expected to be useful to other researchers
working in this area. In addition, this paper constructs interesting and non-trivial examples of networks where the Z-Bus
method exhibits oscillatory behavior, and the conditions are
not satisfied. One example amounts to a two-node network
where the non-slack bus has constant-current and impedance
loads, but a constant-power injection (e.g., from a renewable
energy source). The second example features a nontrivial
unbalanced network with phase couplings in its lines.
Paper Organization: The network and the ZIP load models
with wye and delta configurations are introduced in Section
II. The Z-Bus method is briefly reviewed in Section III.
Section IV presents the main theorem, which contains the
sufficient conditions for the Z-bus method to be a contraction.
The theorem is verified numerically in Section V. Section VII
provides pointers to future work, and Appendices A and B
include a detailed proof of the main result.
Notation: For a vector u, uj denotes the j-th element, and
diag(u) is a square matrix with elements of u on the main
diagonal. The notation v = {vn }n∈N is used to denote a
vector constructed from vertically stacking vectors vn for n ∈
N . For a matrix Z, Z•k denotes the k-th column, Zj• denotes
the j-th row, Zjk denotes the element in the j-th row and k-th
column; and [Z•k ]k∈NY denotes a matrix constructed from the
columns of Z in the set NY . The infinity norm for a vector u
and a P
matrix Z is defined as kuk∞ = maxj |uj | and kZk∞ =
maxj k |Zjk |, respectively.
II. N ETWORK

AND LOAD MODEL

Motivated by the various three-phase connections in IEEE
distribution test feeders, a detailed network model is presented
in this section.

2

A. Three-phase network modeling preliminaries
A power distribution network is represented by a graph
(N , E), where N := {1, 2, . . . , N } ∪ {S} is the set of nodes
and E ⊂ N × N := {(n, m)|n, m ∈ N } is the set of
edges. Node S is considered the slack bus, i.e., the point of
interconnection to the transmission network. Furthermore, let
N be partitioned as N = NY ∪ N∆ ∪ {S} where NY and N∆
collect nodes with wye and delta connections, respectively.
The set of available phases at node n is denoted by Ωn ⊆
{a, b, c}. For wye nodes, i.e., n ∈ NY , Ωn may have one,
two, or three phases. For delta nodes, i.e., n ∈ N∆ , Ωn has at
least two phases. For future use, define r(φ) as the right shift
of phase φ as follows:
r(a) = b, r(b) = c, r(c) = a.

(1)

If |Ωn | = 3, then we have that r(φ) ∈ Ωn for all phases
φ ∈ Ωn . If |Ωn | = 2, then there exists only one phase φ ∈ Ωn
such that r(φ) ∈ Ωn . For example, if Ωn = {a, c}, only phase
c has the property that r(c) = a ∈ Ωn .
For every node n ∈ N , the complex line to neutral voltage
vector is denoted by vn := {vnφ }φ∈Ωn ∈ C|Ωn | . The slack
2π
2π
bus voltage is fixed at vS = {1, e−j 3 , ej 3 }. The line to
neutral voltages are collected in a vector v = {vn }n∈N \{S} .
′
Moreover, we use a vector eφφ
∈ C|Ωn | with entries of 1, -1,
n
′
φφ′ T
and possibly 0, such that (en ) vn = vnφ − vnφ gives the line
′
to line voltages between phases φ and φ . It is easy to see that
′
T
φ′ φ T
(eφφ
n ) vn = −(en ) vn .
further an index set J := {1, . . . , J} where J =
PDefine
N
|Ω
n | is the length of vector v, and each j ∈ J is
n=1
a linear index corresponding to a particular pair (n, φ) with
n ∈ N \ {S} and φ ∈ Ωn . In this case, we write the
operation j = Lin(n, φ). The operation Node[Lin(n, φ)] =
Node[j] = n relates the index j ∈ J to the corresponding
node n. Define the set Jn := {j | Node[j] = n} as the set of
all linear indices j that correspond to node n.SFinally, set J
is partitioned
as J = JY ∪ J∆ , where JY = n∈NY Jn and
S
J∆ = n∈N∆ Jn .
B. Three-phase load models
For a load at node n, the complex vector of net current
injections is denoted by in = {iφn }φ∈Ωn and we define i =
{in }n∈N \{S} . Due to existence of loads, the nodal net current
injection i is a function of nodal voltages v. This dependence is
denoted by in (vn ). According to the ZIP load model, in (vn )
is
iPQn =
 φcomposed of currents from constant-power loads
iPQn φ∈Ωn , constant-current loads iIn = iφIn φ∈Ωn , and

constant-impedance loads iZn = iφZn φ∈Ωn . For n ∈ N \{S}
and φ ∈ Ωn we have that
iφn (vn ) = iφPQn (vn ) + iφIn (vn ) + iφZn (vn )

(2)

where functions iφPQn (vn ), iφIn (vn ), and iφZn (vn ) are defined
in (3) for wye connections and in (4) for delta connections.
Specifically, for wye connections (n ∈ NY and φ ∈ Ωn ),
the ZIP load components are
iφPQn (vn ) = −(sφLn /vnφ )∗

(3a)

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

iφZn (vn ) =

vnφ

iφ ,
φ Ln
|vn |
φ
−yL
vφ ,
n n

iφIn (vn ) = −

3

(3b)
(3c)

where sφLn is the nominal power in the constant-power portion
of the ZIP model, iφLn is the nominal current for the constantφ
current portion of the ZIP model, and yL
is the nominal
n
admittance in the constant-impedance portion of the ZIP model
at node n ∈ NY and phase φ ∈ Ωn .
For delta connections (n ∈ N∆ and φ ∈ Ωn ), the ZIP load
components are
!∗
φφ′
X
sL
φ
n
(4a)
iPQn (vn ) = −
φφ′ T
(e
)
v
n
n
′
φ ∈Ωn \{φ}
′

iφIn (vn )
iφZn (vn )
′

′

=−
=−

X

φ′ ∈Ωn \{φ}

X

′
iφφ
Ln

′

T
(eφφ
n ) vn

′
T
|(eφφ
n ) vn |
′

φφ
T
yL
(eφφ
n ) vn ,
n

(4b)
(4c)

φ′ ∈Ωn \{φ}

where the dependence of injected currents on nodal voltages
is shown explicitly.
Using (2), i(v) can be decomposed into three parts:
i(v) = iPQ (v) + iI (v) + iZ (v),

(7)

where iPQ (v) = {iPQn }n∈N \{S} , iI (v) = {iIn }n∈N \{S} ,
iZ (v) = {iZn }n∈N \{S} . Moreover, (3c) and (4c) reveal that
iZ (v) can be written as a linear function of v as follows:
iZ (v) = −YL v,
where YL ∈ C

J×J

(8)

has entries given by

φ
YL (j, j) = yL
, if j = Lin(n, φ), n ∈ NY ,
n
X
φφ′
, if , j = Lin(n, φ), n ∈ N∆ ,
yL
YL (j, j) =
n
φ′ ∈Ωn \{φ}
′

φφ
YL (j, k) = −yL
, j = Lin(n, φ), k = Lin(n, φ′ ), n ∈ N∆ ,
n

and the remaining entries of YL are all zero.
Substituting (7) and (8) in (6) and after some manipulation
the following fixed-point equation for v is obtained:

′

φφ
φφ
where sφφ
Ln , iLn , and yLn are respectively the nominal power,
current, and admittance in the ZIP model for nodes n ∈ N∆
and over phases
φ, φ′ ′ ∈ Ωn′ . For′ n ∈ N∆ and
φ, φ′′ ∈ Ωn , we
φφ′
φ φ φφ
φφ
φφ′
φφ
have that sLn = sLn , iLn = iLn , and yLn = yL
.
n
Notice that the constant-current loads in (3b) and (4b)
rightfully adhere to the proper definition given in [10], where
it is emphasized that the power varies directly with voltage
magnitude. To see this, one can calculate the apparent power
consumption for the load model in (3b) which gives sφn =
vnφ [(vnφ /|vnφ |)iφLn ]∗ = |vnφ |(iφLn )∗ , i.e., the apparent power is
only a function of the magnitude of the voltage |vnφ |. It is
easy to see that in this case, the load current magnitude is

|iLn | and the load power factor is

Re[iφ
Ln ]
|iφ
Ln |

, both of which are

constant values, concluding that these models are in line with
the definition in [11, pp. 315] as well. It should be noted that
(3b) and (4b) are different than the constant-current-phasor
model employed in [12].
For future use, if k = Lin(n, φ) where n ∈ NY and φ ∈
Ωn , we define skL = sφLn and ikL = iφLn . For k = Lin(n, φ)
where n ∈ N∆ and φ ∈ Ωn , if r(φ) ∈ Ωn , define skL =
φr(φ)
φr(φ)
φr(φ)
. If r(φ) ∈
/ Ωn , then
sLn , ikL = iLn , and ek = en
′
′
k
Ωn = {φ, φ }, and define sL = 0, ikL = 0, and ek = eφφ
n .

v = Z [iPQ (v) + iI (v)] + w,
−1

where Z = (Y + YL )

(6)

(11)

Equation (10) lends itself to the well-known Z-Bus method
for three-phase networks, outlined in e.g., [1]:
v[t + 1] = Z [iPQ (v[t]) + iI (v[t])] + w,

(12)

where v[t] is the value of v at iteration t. Notice that (12) is an
application of Picard’s iteration for solving nonlinear systems
of equations [14, pp. 182].
IV. S UFFICIENT CONDITIONS FOR CONVERGENCE OF THE
Z-B US METHOD TO A UNIQUE SOLUTION
In this section, sufficient conditions are presented such that
the iteration of the Z-Bus method (12) is contracting, i.e., the
iteration converges geometrically to a unique solution of (10).
To present these conditions, it is beneficial to make a change
of variables in (12). Specifically, let u = Λ−1 v[t] where Λ ∈
CJ×J is an invertible diagonal matrix with entries λj . Matrix
Λ will serve as a design parameter which can be leveraged
to potentially expand the convergence region of the Z-Bus
method. Parameterizing v by u and multiplying both sides of
(12) by Λ−1 yield the following iteration for u:
u[t + 1] = T(u[t]),

where iS is the complex current injection of the slack bus.
Matrices Y, YNS , YSN , and YSS are formed by concatenating
the admittance matrices of transmission lines, transformers and
voltage regulators, given e.g., in [13]. It follows from (5) that
i(v) = Yv + YNS vS ,

, and

w = −ZYNS vS .

III. T HE Z-B US M ETHOD
For a three-phase network, the multidimensional Ohm’s law
is given by
 
  
v
Y
YNS
i
,
(5)
=
YSN YSS vS
iS

(10)

J

(13)

J

where T(u) : C → C is the following mapping

T(u) = Λ−1 Z [iPQ (Λu) + iI (Λu)] + Λ−1 w,

(14)

and the Z-Bus method of (12) can be equivalently written as
v[t + 1] = ΛT(Λ−1 v[t]).

(15)

Since the iterations in (13) and (15) have a one-to-one correspondence, convergence of (13) to a fixed point ufp = T(ufp )
is equivalent to the convergence of (15) to a respective solution
vfp = Λufp of (10).

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

4

In what follows, conditions such that the mapping T(u)
in (14) is a contraction over a region D ⊆ CJ are derived.
As a consequence, T(u) has a unique fixed point over D, and
iterations (13) converge geometrically to the unique fixed point
[15, Prop. 3.1.1]. The definition of contraction is as follows:
Definition. For a closed set D ⊆ CJ , a mapping T(u) is a
contraction on D if the following two properties hold:
1) Self-mapping property: u ∈ D ⇒ T(u) ∈ D.
2) Contraction property: u, ũ ∈ D ⇒ kT(u) − T(ũ)k ≤
αku − ũk,
where k.k is some norm and α ∈ [0, 1) is a constant.

kΛ−1 ZY diag(sY )diag(wY )−1 k∞
1 − Rλ̄/w
−1 ∆
kΛ Z diag(s∆ )diag(w∆ )−1 k∞
+
1 − 2Rλ̄/ρ

+ kΛ−1 ZY diag(iY )k∞ + kΛ−1 Z∆ diag(i∆ )k∞ ≤ R (C3)
kΛ−1 ZY diag(sY )ΛY [diag(wY )−1 ]2 k∞
(1 − Rλ̄/w)2
−1 ∆
2kΛ Z diag(s∆ )Λ∆ [diag(w∆ )−1 ]2 k∞
+
(1 − 2Rλ̄/ρ)2

2kΛ−1 ZY diag(iY )ΛY diag(wY )−1 k∞
1 − Rλ̄/w
4kΛ−1 Z∆ diag(i∆ )Λ∆ diag(w∆ )−1 k∞
+
≤ α < 1.
1 − 2Rλ̄/ρ
+

In regards to the previous definition, any number α ∈ [0, 1)
˜
u)k
uniformly for all
that upperbounds the ratio kT(u)−T(
ku−ũk
u, ũ ∈ D (u 6= ũ), is called contraction modulus.

In order to state the conditions for (14) to be a contraction
map, the following quantities are defined:

(C4)

sY = {skL }k∈JY ∈ C|JY | , iY = {ikL }k∈JY ∈ C|JY | , (16b)
wY = {wk }k∈JY ∈ C|JY | , ZY = [Z•k ]k∈JY CJ×|JY | , (16c)

If conditions (C1)–(C4) hold for some R > 0, then T(u) is
a contraction mapping and has as a consequence a unique fixed
point ufp in the ball defined by DR . In other words, the power
flow equations have a unique solution in DR . The value of R
can be interpreted as the infinity norm of the voltage difference
from the no-load power flow solution scaled by Λ. The latter
is obtained as u = Λ−1 w upon setting iPQ = iI = 0 in (14).

w∆ = {eTk wNode[k] }k∈J∆ ∈ C|J∆ | ,

Moreover, the contraction property of T(u) implies the
following relationships for successive iterates of the Z-Bus
method upon initialization with u[0] ∈ DR :

w = min |wk |, λ̄ = max |λk |, ρ = min {|eTk wNode[k] |}, (16a)
k∈J

k∈J

k∈J∆

ΛY = diag({λk }k∈JY ) ∈ C|JY |×|JY | ,
s∆ = {skL }k∈J∆ ∈ C|J∆ | , i∆ = {ikL }k∈J∆ ∈ C|J∆ | ,

(16d)
(16e)

J×|J∆ |

(16g)

∆

Z =

[Z∆
•lk ]k∈J∆

∈C

,

(16f)

where lk ∈ {1, . . . , |J∆ |} is the order of k in J∆ , and

Z•k − Z•k′ , if φ, r(φ) ∈ Ωn ,
Z∆
=
•lk
0J ,
if φ ∈ Ωn , r(φ) ∈
/ Ωn ,
′
with k = Lin(n, φ), k = Lin(n, r(φ)),
∆
Λ = diag({maxl∈JNode[k] |λl |}k∈J∆ ) ∈ C|J∆ |×|J∆ | . (16h)
In (16), the notation wn = {wk }k∈Jn relates the indices in
w to corresponding ones in v. Respectively for wye and delta
connections, quantities sY and s∆ collect all constant-power
loads, iY and i∆ collect all constant-current loads, ZY and
Z∆ collect corresponding columns of the impedance matrix
Z, and wY , w∆ collect the voltages induced by the slack bus.
Matrix ΛY contains the elements of the design matrix Λ that
correspond to wye nodes and matrix Λ∆ selects the maximum
of the per node entries of the design matrix Λ for delta nodes.
Theorem 1 establishes the convergence of the Z-Bus method
in three-phase distribution networks with ZIP loads, and is
proved in Appendix B using intermediate results in Appendix A.
Theorem 1. Define the closed ball DR := {u ∈ CJ | ku −
Λ−1 wk∞ ≤ R} where R > 0. Then T(u) in (14) is a
contraction mapping on DR with contraction modulus α if
the following four conditions hold:
1 − Rλ̄/w > 0,

(C1)

1 − 2Rλ̄/ρ > 0,

(C2)

ku[t + 1] − ufp k∞ ≤ αku[t] − ufp k∞ ,

ku[t + 1] − u[t]k∞ ≤ αku[t] − u[t − 1]k∞ .

(17a)
(17b)

Equation (17a) implies that the iterates in (13) converge to ufp
with a convergence rate that is upperbounded by α. Equation
(17b) can be used to numerically evaluate the convergence
ku[t+1]−u[t]k∞
. We will see in
rate by computing the ratio ku[t]−u[t−1]k
∞
Section V for the IEEE test feeders that the aforementioned
ratio turns out to be much smaller than the estimate provided
by the left-hand side of (C4).
Theorem 1 on T(u) implies certain convergence properties
for iteration (12), as stated next.
Theorem 2. If conditions (C1)–(C4) hold, then by initializing
v[0] ∈ D′R := {v|kΛ−1 (v−w)k∞ ≤ R}, the Z-Bus iterations
in (12) remain within D′R and converge to a unique solution.
Moreover, the convergence is R-linear with rate α [i.e., the
sequence generated by (12) is dominated by the geometric
sequence Bαt where B is a non-negative scalar].
Proof. The self-mapping property of T(u) and the one-toone correspondence between v and u guarantee that the
iterates v[t] remain within D′R . Regarding the convergence
rate of (12), note that the contraction property of T(u) implies
(cf. [15, Prop. 3.1.1])
ku[t] − ufp k∞ ≤ αt ku[0] − ufp k∞ .
Since we have u = Λ

−1

(18)

v, it follows from (18) that

kΛ−1 (v[t] − vfp )k∞ ≤ αt kΛ−1 (v[0] − vfp )k∞ . (19)

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

5

From the definition of the infinity norm it follows that
1
(vj [t] − vjfp )
λj


1
= max
|vj [t] − vjfp |
j∈J
λj
1
≥ max |vj [t] − vjfp |
λ̄ j∈J
1
(20)
⇒kΛ−1 (v[t] − vfp )k∞ ≥ kv[t] − vfp k∞ .
λ̄
Moreover, we have that
kΛ−1 (v[t] − vfp )k∞ = max
j∈J

1
(vj [0] − vjfp )
λj
1
(vj [0] − vjfp )
≤ max
j∈J mink∈J |λk |
1
≤
max |vj [0] − vjfp |
mink∈J |λk | j∈J
1
kv[0] − vfp k∞ . (21)
⇒kΛ−1 (v[0] − vfp )k∞ ≤
mink∈J |λk |
kΛ−1 (v[0] − vfp )k∞ = max
j∈J

Combining (20), (21), and (19) yields

kv[t] − vfp k∞ ≤ Bαt ,

(22)

λ̄

where B = minj∈J |λj | kv[0] − vfp k∞ ≥ 0. Since the quantity
kv[t] − vfp k∞ is upperbounded by a geometric sequence,
the rate of convergence is (informally) said to be geometric;
formally, the convergence is R-linear with rate α [16, pp.
619].
Conditions (C1)–(C4) can be expressed via inequalities
involving polynomials of degree at most four in the radius
of the contraction region. The coefficients of the polynomials
can be readily written in terms of the network parameters
and loads. Quartic polynomials and their roots have been
completely characterized, and can be routinely computed in
terms of the polynomial coefficients—see for example [17].
The convergence region can thus be easily computed. Significant analytical advantages can also be afforded by the
presented conditions, as the convergence region can potentially
be expressed explicitly in terms of the network loads and
parameters. An alternative set of nonlinear non-polynomial
conditions for convergence of the Z-Bus method is presented
in [18].
Conditions (C1)–(C4) can also be graphically solved, as
explained next. First, the quantities in (16) must be computed
which depend only on the feeder data. The set of R values
satisfying each condition can simply be plotted on the real
line. The intersection of the feasible regions in these plots
gives a range for R. In the next section, the valid range
of R satisfying conditions (C1)–(C4) is computed for IEEE
distribution test feeders, and the convergence of the Z-Bus
method is numerically illustrated.
V. N UMERICAL

VERIFICATION ON

IEEE FEEDERS

In this section, the range of R where conditions (C1)–(C4)
are satisfied is given for IEEE-13, IEEE-37, and IEEE-123
distribution test feeders [19]. Furthermore, the self-mapping

property as well as the convergence of the Z-Bus method are
numerically confirmed for iterations on voltages in (12). This
is because the Z-Bus method is routinely implemented on the
variable v as opposed to u.
The feeder data are available online [20]. We model
transformer nodal admittances according to [1]. However, to
avoid singularities in the bus admittance matrix, delta-delta
transformers are modeled using a modification suggested in
[21] which entails connecting a small shunt admittance to the
ground. For voltage regulators, the models are derived from
[11, Ch. 7].1
To numerically verify conditions (C1)–(C4), the design
parameter Λ is selected to be Λ = diag(w). This choice
allows for the center of DR to be the vector of all one’s. This
practically means that the no-load solution for scaled voltages
is u = Λ−1 w = 1.
For the IEEE-13 test feeder, conditions (C1)–(C4) yield the
following inequalities for R:
1 − 1.043R > 0; 1 − 1.203R > 0
0.041
0.127
+
+ 0.032 + +0.01 ≤ R
1 − 1.043R 1 − 1.203R
0.127
0.048
+
(1 − 1.043R)2 (1 − 1.203R)2
0.024
0.064
+
< 1.
+
1 − 1.043R 1 − 1.203R
For the IEEE-37 test feeder, the conditions are
1 − 1.037R > 0; 1 − 1.198R > 0
0.037
+ 0.018 ≤ R
1 − 1.198R
0.042
0.043
+
< 1.
(1 − 1.198R)2 1 − 1.198R

(23a)
(23b)

(23c)

(24a)
(24b)
(24c)

Finally, for the IEEE-123 feeder, the conditions are
1 − 1.078R > 0; 1 − 1.238R > 0
(25a)
0.001
0.129
+
+ 0.030 + 0.012 ≤ R
(25b)
1 − 1.078R 1 − 1.238R
0.001
0.129
+
2
(1 − 1.0777R)
(1 − 1.238R)2
0.060
0.027
+
+
< 1. (25c)
1 − 1.078R 1 − 1.238R
For each feeder, the subsets of the real line where conditions
(C1)–(C4) are satisfied are depicted in Fig. 1, by solving
(23), (24), and (25). The intersection of these regions for the
IEEE-13, IEEE-37 and IEEE-123 feeders is [Rmin , Rmax ] =
[0.29, 0.48], [Rmin , Rmax ] = [0.06, 0.64], and [Rmin , Rmax ] =
[0.22, 0.54], respectively. Both Rmin and Rmax reveal important information about the Z-Bus convergence. Specifically,
larger Rmax guarantees convergence when the initialization is
far from the unique solution. Smaller Rmin guarantees that the
unique solution is close to the center of the ball. The latter
implies that the deviation of the load-flow voltage solution
from the no-load solution is tightly bounded. Notice that these
regions are dependent on the design matrix Λ = diag(w).
1 The MATLAB scripts that compute the convergence regions for the IEEE
test feeders and run the Z-Bus method are provided online at the following
page: https://github.com/hafezbazrafshan/contraction-mapping.

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

(C1)

(C2)

(C3)

6

vn [0] = vS
vn [0] = wn
vn [0] = wn − Rmin wn
vn [0] = wn − Rmax wn
Rmin
Rmax

(C4)
0.5

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kΛ−1 (v[t] − w)k∞

0

R
(a)
(C1)

(C2)

(C3)

(C4)

0.4
0.3
0.2
0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9
0

R
(b)

0

(C1)

(C2)

(C3)

1

2

3

4

Iteration t

(C4)

Demonstration of the self-mapping property of the Z-Bus
iterations in (12) for IEEE-123 test feeder. As long as the Z-Bus
method is initialized with v[0] ∈ D′R , where R satisfies conditions
(C1)–(C4), successive iterations v[t] remain within D′R .
Fig. 3.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

R
(c)

Fig. 1. Regions on the real line where conditions (C1)–(C4) hold

for (a) IEEE-13, (b) IEEE-37, and (c) IEEE-123 test feeders. The
presence of a mark indicates that the condition is satisfied.

vn [0] = vS
vn [0] = wn
vn [0] = wn − Rmin wn
vn [0] = wn − Rmax wn

103

101

kv[t + 1] − v[t]k∞

Lower bound of contraction modulus α

1

IEEE-13
IEEE-37
IEEE-123

0.9
0.8

10−1

10−3

0.7
0.6

10−5
0.5
0.4

10−7

0.3

0

1

2

3

4

Iteration t

0.2
0.1

Fig. 4. Convergence of the Z-Bus iterations in (12) for IEEE-123 test

0
0

0.1

0.2

0.3

0.4

0.5

0.6

R

Fig. 2. Lower bound of contraction modulus α versus R in the range for
which conditions (C1)–(C4) are satisfied.

For example, setting Λ = IJ for the IEEE-13 feeder yields
Rmax = 0.60 and also shifts the center of the ball from 1
to w. A study on the impact of Λ on the contracting region
remains for future work.
Figure 2 depicts the lower bound of the contraction modulus
α evaluated from (C4) as a function of the feasible R for
the three IEEE feeders. It is revealed in Fig. 2 that a more
localized initialization (i.e., smaller R) yields smaller α, with
the smallest α occurring at Rmin . For the IEEE-13, IEEE37, and IEEE-123 feeders, the smallest α is respectively
0.50, 0.09, and 0.34. Upon initializing within DRmin , the
convergence rate is guaranteed to be at most as much as the
aforementioned values.
The self-mapping property of the Z-Bus iterations (12) is
demonstrated in Fig. 3 for the IEEE-123 test feeder. It is
shown that when the initialization v[0] is within D′R , where
R ∈ [Rmin , Rmax ], the sequence v[t] produced by the ZBus method in (15) remains within D′R . In Fig. 3, a case of
special interest is the graph corresponding to the initialization
with the flat voltage profile, that is, when vn [0] = vS , where

feeder. It is observed that the Z-Bus method converges to the unique
solution very fast—within 4 iterations.

π

π

vS = {1, e−j2 3 , ej2 3 } is the voltage at the slack bus. The
graph shows that the typical initialization with the flat voltage
profile is within the contracting region; and upon initializing
with vn [0] = vS ∈ D′R for all nodes, the iterates remain
within D′R .
The distance between two successive iterates of the Z-Bus
method for the IEEE-123 test feeder is plotted in Fig. 4.
It is shown that the distance kv[t + 1] − v[t]k∞ practically
converges to zero within 4 iterations, implying that the fixedpoint solution is obtained. For the sequences of Fig. 4, the
empirical convergence rate is also numerically calculated using
(17b) and is depicted in Fig. 5. All the numerical convergence
rates are smaller than 0.34, which is the smallest theoretically
obtained contraction modulus for the IEEE-123 network.
The merits of this work is that it characterizes certain
regions with unique power flow solution in three-phase distribution networks and that it guarantees convergence of the
Z-Bus method. Working with networks that satisfy solvability
conditions can provide a degree of assurance that the numerical
algorithm—the Z-Bus method in this case—will converge to
a solution.
Further research is required to obtain conditions that ac-

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

7

0.12

vn [0] = vS
vn [0] = wn
vn [0] = wn − Rmin wn
vn [0] = wn − Rmax wn

kΛ−1 (v[t+1]−v[t])k∞
kΛ−1 (v[t]−v[t−1])k∞

0.1

0.08

0.06

0.04
1

2

3

4

Iteration t

Fig. 5. Numerical convergence rate of the Z-Bus iterations in (12) for
the iterations in Fig. 4. The numerical convergence rate at iteration t
is calculated using (17b). Notice that these convergence rates are all
smaller than 0.34, which is the tightest upper bound given in Fig. 2
for IEEE-123.

curately characterize the largest region with a unique power
flow solution. As was previously alluded to, it is possible to
leverage the design parameter Λ for this purpose. By setting
Λ = diag(w), the conditions are satisfied for the practical
IEEE test networks; developing comprehensive methods for
the selection of parameter Λ are left for future work.
It should be noted at this point that this work derives
sufficient conditions for convergence. Generally speaking, sufficient conditions are conservative by nature. That is, the existence of loads and network parameters that do not satisfy the
conditions yet yield convergence is not precluded. It is clear
from conditions (C3) and (C4) that by increasing the constantpower and constant-current loads, there will eventually be no
R > 0 that ensures contraction.
As an example, in the IEEE-123 test feeder, after 50% of
load increase, no R > 0 satisfies condition (C3), however the
Z-Bus converges to a solution within 5 iterations. Therefore,
the ultimate goal would be to provide a set of conditions that
divide the voltage space to two regions, namely, one where
there is convergence and another one of non-convergence.
However, this is still an open research question. The next
section provides two network examples where the Z-Bus
method fails to converge in cases where conditions (C1)–(C4)
are not satisfied yet a solution may exist.

VI. N UMERICAL

EXAMPLES OF NON - CONVERGENCE

In this section, we provide two example networks where the
Z-Bus method does not converge. The first example (Section
VI-A) depicts a simple two-node network for which an analytical solution is found; however, convergence conditions are not
satisfied and the Z-Bus method oscillates. The second example
(Section VI-B) features a non-trivial unbalanced network with
constant-power loads where the convergence occurs only when
the loads are reduced by a certain percentage. The convergence
conditions also are satisfied for a close percentage, which
reveals that they are relatively tight in this example.

A. Balanced two-node network with power injection
Consider a two-node network where the slack bus S is connected to a secondary bus through an ideal phase-decoupled
three-phase line with real-valued series impedance yt . The
secondary bus contains real-valued ZIP components sL , iL ,
and yL on each phase, where iL , yL > 0 but sL < 0, that is,
the constant-power portion of the ZIP component is considered
to be a power injection. In this case, we have that sL = sL 1,
iL = iL 1, YL = yL I where 1 ∈ R3 is a vector of all one’s and
I ∈ R3×3 is the 3 × 3 identity matrix. Moreover, the network
admittance matrix is Y = −YNS = yt I. Denote the voltage
vector at the non-slack bus with v = {v a , v b , v c } and let the
2π
2π
slack bus voltage be fixed at vS = {1, e−j 3 , ej 3 }.
Due to the fact that the network is balanced and the values
of loads are real, it turns out that a solution of the form
2π
2π
v a = v b ej 3 = v c e−j 3 = v exists where v is a real number.
To demonstrate this, we write the Ohm’s law for phase a as
follows:
v
sL
− iL
= (yt + yL )v − yt .
(26)
−
v
|v|
Assuming that we are only interested in non-negative voltages v, that is, when |v| = v, the following quadratic equation
can be obtained for v:
(yt + yL )v 2 + (iL − yt )v + sL = 0.
The solution to (27) is explicitly given as
√
yt − i L
∆
v=
±
,
2(yt + yL ) 2(yt + yL )

(27)

(28)

with ∆ = (yt − iL )2 − 4(yt + yL )sL . Keep in mind that (28)
is not necessarily the only admissible voltage solution for this
network since we have only limited the search to v ∈ R≥0 .
In order to simplify the ensuing computations, let us select
yt = yL = 21 pu, and set the current injection iL = 21 pu.
Substituting these values into (28), the real non-negative
voltage solution v parameterized by sL is thus computed
√
v = −sL .
(29)
It is clear that for values sL √
< 0, the network admits at least
one voltage solution vsol = −sL vS .
At this point, it is our intention to show that even though
voltage solutions exist in this example network, the Z-Bus
method may not converge and the convergence conditions
(C1)–(C4) are not satisfied.
1
We first compute Z = (Y +YL )−1 = yt +y
I = I and w =
L
1
2 vS . By selecting Λ = I, convergence conditions (C1)–(C4)
for this network yield the following inequalities parameterized
by sL :
1 − 2R > 0,
1
2|sL |
+ ≤ R,
1 − 2R 2
4|sL |
2
+
< 1,
(1 − 2R)2
1 − 2R

(30a)
(30b)
(30c)

where in (30), due to the absence of delta-connected nodes,
(C2) is not present.

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

2

8

v[0] = vs
v[0] = 2vsol
v[0] = 0.5vsol
v[0] = vsol + 0.1vS

1.8

101

kv[t + 1] − v[t]k∞

1.6

102

kv[t]k∞

1.4
1.2
1
0.8
0.6

100
10−1
10−2
10−3
10−4

0.4

10−5

0.2

10−6

0
0

2

4

6

8

0

10

Iterations of the Z-Bus method performed on the two-node
test network while setting sL = −0.5 pu. Since Z-Bus convergence
conditions in (30) are not feasible for any sL , the convergence of the
Z-Bus method is not guaranteed, even when there is knowledge that
solution (29) exists.

Inequality (30b) can be transformed to the quadratic inequality 2R2 − 2R + 12 + 2|sL | ≤ 0, which does not have
a feasible solution in R for any sL < 0. However, existence
of power solutions for sL < 0 is easily perceived by (29).
Figure 6 plots the infinity norm of the iterates of the Z-Bus
method, that is, kv[t]k∞ , when sL = −0.5 pu. Figure 6 is
representative of the fact that the Z-Bus method does not
converge to the solution (29).
B. Unbalanced three-node network

2

4

6

8

10

12

14

16

18

20

Iteration t

Iteration t

Fig. 6.

θ=1
θ = 0.5
θ = 0.3
θ = 0.2
θ = 0.12
θ = 0.11

Iterations of the Z-Bus method performed on the threenode test network for several values of parameter θ. When the
value of θ is such that Z-Bus convergence conditions (33) are not
feasible, the convergence of the Z-Bus is not guaranteed. The black
line corresponding to θ = 0.12 depicts a very long oscillation.
As θ approaches θsol = 0.1060, the Z-Bus method converges.
sol
When
, a√ unique solution exists within DR where R ∈
√ θ ≤ θ
1− 1−9.4360θ 1+ 1−9.4360θ
[
,
].
2
2
Fig. 7.

the following inequalities in R parameterized by θ:
1−R>0⇒R<1
(33a)
2.359θ
≤ R ⇒ R2 − R + 2.359θ ≤ 0
(33b)
1−R
2.359θ
< 1 ⇒ R2 − 2R + 1 − 2.359θ > 0. (33c)
(1 − R)2

It is not hard to see that (33) has a valid solution for R only
sol
√ Z-Bus is
Consider a three-node network with the set of nodes if θ ≤ θ = 0.1060. The range of √R for convergent
1− 1−9.4360θ 1+ 1−9.4360θ
,
].
N = {1, 2, S} and the set of edges {(1, S), (1, 2)}. The series then given as [Rmin , Rmax ] = [
2
2
The Z-Bus method is run on this example network for several
admittances for edges (1, S) and (1, 2) are as follows:
representative values of θ using the initialization v[0] = w.


0.077 − j5.33 0.01 − j0.09 0.02 − j0.08
The corresponding Z-Bus iterations are shown in Fig. 7. It is
0.087 − j8 0.03 − j0.07 , (31a) revealed in Fig. 7, that for almost all percentages θ higher than
Y1S =  0.01 − j0.09
0.02 − j0.08 0.03 − j0.07 0.07 − j1.5
θsol the Z-Bus method fails to converge. However, for θ ≤ θsol ,


0.056 − j8.66
0
0.02 − j0.07
the Z-Bus iterations are guaranteed
to converge
√
√ to the unique
1+ 1−9.4360θ
0
0.02 − j4.8 0.03 − j0.05 , (31b) solution in DR where R ∈ [ 1− 1−9.4360θ
Y12 = 
,
].
2
2
0.02 − j0.07 0.03 − j0.05 0.03 − j3.8

where all values are per unit. The load buses {1, 2}
only contain
values of
 constant-power loads with per unit
T
and sL2 =
sL1 = θ 0.7 + j1.5 0.8 + j1.5 0.8 + j2.5

T
θ 0.6 + j2.5 0.6 + j0.5 0.3 + j0.5 where θ ∈ (0, 1] is
a parameter that is going to be varied to obtain convergence
threshold for the Z-Bus method.
Our intention is to find θ such that the convergence of the
Z-Bus method is guaranteed and numerically investigate how
conservative the theoretical bound is. First, network parameters
are computed:


Y1S + Y12 −Y12
Y=
(32a)
−Y12
Y12


 
−Y1S
v
YNS =
(32b)
, w= S .
O
vS
Next, to find θ such that Z-Bus convergence is guaranteed,
by setting Λ = I, the convergence conditions (C1)–(C4) yield

VII. C ONCLUDING R EMARKS

A set of sufficient conditions was derived that guarantee
the convergence of the Z-Bus method to a unique solution of
the power flow problem in three-phase distribution networks
with ZIP loads. It was numerically demonstrated that the IEEE
test feeders with wye and delta ZIP loads satisfy the derived
conditions. Characterizing the largest regions with unique
power flow solution is an open research question. Investigating
the effect of the design parameter Λ to this end is the subject
of future work.
A PPENDIX A
U SEFUL I NTERMEDIATE R ESULTS
This Appendix proves intermediate results on the product
Z[iPQ +iI ] that shows up in (12), and an inequality on complex
numbers, which will be later used.

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

A. Rearranging

P

k∈J∆

9

Zjk iPQ (v)k

B. Rearranging

∆

Lemma 1. Consider the matrix Z defined in (16g). For n ∈
N∆ , φ ∈ Ωn , k = Lin(n, φ) ∈ J∆ , j ∈ J , and indices
lk ∈ {1, . . . , |J∆ |}, the following equality holds:
 k ∗
X
X
sL
Z∆
Zjk iPQ (v)k = −
.
(34)
jlk
Tv
e
k n
k∈J
k∈J
∆

∆

Proof. Using (4a) for constant-PQ loads with delta connections, we can write
X
Zjk iPQ (v)k =
k∈J∆

−

X

n∈N∆




′∗

X

ZjLin(n,φ)

φ∈Ωn

sφφ
Ln

X

′

T ∗
(eφφ
n ) vn

φ′ ∈Ωn /{φ}



.

(35)

Consider the term in parenthesis in (35). Since n ∈ N∆ , then
|Ωn | ≥ 2, i.e., node n has at least two available phases. If
|Ωn | = 2, then without loss of generality we can assume
Ωn = {φ, φ′ } where φ′ = r(φ) and r(φ′ ) ∈
/ Ωn . The term
in parenthesis in (35) can be written as follows:
′∗

ZjLin(n,φ)
=
=

sφφ
Ln
′

T ∗
(eφφ
n ) vn

′

+ ZjLin(n,φ′ )

sφLnφ

∗

(36)

where k ′ = Lin(n, r(φ)), and equality (a) results from the
fact that Z∆
/ Ωn [cf. (16g)].
jlk is zero when φ ∈ Ωn and r(φ) ∈
If node n has three phases then Ωn = {a, b, c} and r(a) = b,
r(b) = c, and r(c) = a. The term in parenthesis in (35) equals:
∗

∗

sac
sab
ZjLin(n,a) abLnT ∗ + ZjLin(n,a) acLnT ∗
(en ) vn
(en ) vn
∗
∗
sbc
sba
+ ZjLin(n,b) baLnT ∗ + ZjLin(n,b) bcLTn ∗
(en ) vn
(en ) vn
∗
∗
sca
scb
+ ZjLin(n,c) caLnT ∗ + ZjLin(n,c) cbLTn ∗
(en ) vn
(en ) vn
∗
ab∗
sL n
sab
Ln
= ZjLin(n,a) ab T ∗ − ZjLin(n,b) ba T ∗
(en ) vn
(en ) vn
∗
bc∗
sbc
sL n
Ln
+ ZjLin(n,b) bc T ∗ − ZjLin(n,c) bc T ∗
(en ) vn
(en ) vn
∗
∗
sca
sca
Ln
Ln
+ ZjLin(n,c) ca T ∗ − ZjLin(n,a) ca T ∗
(en ) vn
(en ) vn
∗
X 
 sφr(φ)
Ln
=
ZjLin(n,φ) − ZjLin(n,r(φ))
φr(φ) T ∗
(en
) vn
φ∈Ωn
∗
∗
k
k
X
X
s
sL
[Zjk − Zjk′ ] TL ∗ =
=
Z∆
jlk T ∗ .
ek vn
ek vn
k∈J
k∈J
n

k∈J∆

Zjk iI (v)k

Lemma 2. Consider the matrix Z∆ defined in (16g). For n ∈
N∆ , φ ∈ Ωn , k = Lin(n, φ) ∈ J∆ , j ∈ J , and indices
lk ∈ {1, . . . , |J∆ |}, the following equality holds:
X
X
ikL eTk vn
Z∆
.
(38)
Zjk iI (v)k =
jlk
|eTk vn |
k∈J
k∈J
∆

n

′∗

Proof. Follow Appendix A-A with the terms
∗
sk
L
∗
eT
v
k n

replaced by

′
φφ′ T
iφφ
Ln (en ) vn
′
T
|(eφφ
n ) vn |

and

T
ik
L ek vn
|eT
k vn |

sφφ
Ln
′

T ∗
(eφφ
n ) vn

and

, respectively.

C. A useful inequality for complex numbers
Lemma 3. For two complex numbers x, y 6= 0 we have that
y
|x − y|
x
−
.
≤2
|x| |y|
|x|

(39)

Proof. The following hold:

x
y
y
y
y
x
≤
+
−
−
−
|x| |y|
|x| |x|
|x| |y|
|y| − |x|
|x − y|
|x − y|
=
+ |y|
≤2
. (40)
|x|
|x||y|
|x|

Notice that (40) holds since ||y| − |x|| ≤ |x − y|.

′

(eφn φ )T vn∗

φr(φ)∗
φr(φ)∗
sL n
sL n
ZjLin(n,φ) φr(φ)
− ZjLin(n,r(φ)) φr(φ)
(en
)T vn∗
(en
)T vn∗
∗
∗
∗
sk
skL (a) X ∆ skL
Z
[Zjk − Zjk′ ] TL ∗ = Z∆
=
jlk T ∗ ,
jlk T ∗
ek vn
ek vn
ek vn
k∈Jn

P

A PPENDIX B
P ROOF OF T HEOREM 1
In this section, the proof of Theorem 1 is given. First,
to guarantee that for all u ∈ DR the corresponding lineto-neutral and line-to-line voltage magnitudes are positive,
conditions (C1) and (C2) are derived in Appendices B-A and
B-B respectively. Condition (C3) ensures the self-mapping
property of T(u) and is derived in Appendix B-C. Appendix
B-D derives condition (C4), which guarantees the contraction
property of T(u).
A. Bounds on line to neutral voltages
This section proves two results. First, a lower bound on all
line to neutral voltages is derived in (42) which will be later
used in proving the self-mapping and contraction properties.
Second, condition (C1) is proved.
Due to condition u ∈ DR , the following holds for k ∈ J :
wk
wk
| ≤ R ⇒ |uk −
| ≤ R,
max |uk −
k∈J
λk
λk
⇒ |λk uk − wk | ≤ |λk |R ⇒ |vk − wk | ≤ R|λk |,
⇒ |wk | − |vk | ≤ R|λk | ⇒ |wk | − R|λk | ≤ |vk |
⇒ |wk |(1 − R|λk |/|wk |) ≤ |vk |.

(37)

n

For |Ωn | = 2 or |Ωn | = 3, substituting the appropriate term
from (36) or (37) for the parenthesis in (35) yields (34).

(41)

(42)

Equation (42) provides a lower bound for the magnitude of
line to neutral voltages in terms of R and w.
The proof of (C1) is based on (42). In (42), in the term in
parenthesis, we replace |wk | by its minimum w [defined in
(16a)], and replace |λk | by its maximum λ̄ [also defined in
(16a)] to obtain a lower bound for the left-hand side of (42)
which is common for all k ∈ J :
|wk |(1 − R|λk |/|wk |) ≥ |wk |(1 − Rλ̄/w).

(43)

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

10

To ensure that for all u ∈ DR , the magnitude of the line to
neutral voltages is positive, the right-hand side of (43) must
be greater than zero, yielding condition (C1).

The right-hand side of (49) is the sum of four terms. An
upper bound for each term is provided next.
For the first term in (49), the following holds:
skL ∗
)
vk
k∈JY
k∈JY
" ∗  ∗  ∗ #
X
skL
skL
skL
=
Zjk
−
+
vk
wk
wk

B. Bounds on line to line voltages
Similar to Appendix B-A, this section proves two results.
First, we find a lower bound for line to line voltages in delta
connections. Next we derive condition (C2).
For n ∈ N∆ , any line to line voltage can be written as
eTk vn , where k = Lin(n, φ) and φ ∈ Ωn . Using the triangle
inequality, we find the following:
|eTk vn | = |eTk (vn − wn ) + eTk wn |
≥

|eTk wn |

−

|eTk (vn

X

≤
≤

(44)

Due to Hölder’s inequality, for n = Node[k] we have that

=

|eTk (vn − wn )| ≤ kek k1 kvn − wn k∞ ,
(a)

l∈Jn

⇒ |eTk (vn − wn )| ≤ 2 max{|λl ||ul −
l∈Jn

wl
|}
λl

(b)

⇒ |eTk (vn − wn )| ≤ 2R max{|λl |},
l∈Jn

(45)

where (a) holds because kek k1 = 2 for k ∈ J , and (b) holds
due to the assumption u ∈ DR . Using (45) in (44) we find:
|eTk vn |

≥

|eTk wn |

− 2R max{|λl |}

(46)

l∈Jn

≥ |eTk wn |(1 − 2R max{|λl |}/|eTk wn |).
l∈Jn

l∈Jn

(48)

To guarantee that for all u ∈ DR , the magnitude of line-to-line
voltages is positive, (48) must be greater than zero, yielding
condition (C2).
C. Self-mapping property
To obtain a sufficient condition for the self-mapping property, it is assumed that u ∈ DR , and an upper bound for
the term kT(u) − Λ−1 wk∞ is sought. The j-th entry of
kT(u) − Λ−1 wk∞ can be upper-bounded as follows:

≤

1
|λj |

+

X

k∈JY

wj
1
wj
wj
|=
Zj• [iPQ (v) + iI (v)] +
−
λj
λj
λj
λj
X
X
Zjk iPQ (v)k
Zjk iPQ (v)k +

|Zjk |

k∈JY

X

X

k∈J∆

Zjk iI (v)

!

.

k∈JY

|skL |R|λk |

(|wk | − R|λk |) |wk |

|Zjk ||skL |

(b)

|wk | − R|λk |

≤
≤

X

|Z∆
jlk |

X

|Z∆
jlk |

k∈J∆

k∈J∆

≤

X

k∈JY

X

k∈JY

|Zjk |

|skL |
|wk |

|Zjk ||skL |
,
|wk |(1 − Rλ̄/w)

(50)

X
|skL ||eTk (wn − vn )|
skL
|Z∆
+
jlk |
T
T
T
|ek vn ||ek wn |
e k wn
k∈J
∆

|skL |2R max{|λl |}
l∈Jn

(|eTk wn | − 2R max{|λl |})|eTk wn |
l∈Jn

skL
|Z∆
|
+
jlk
eTk wn
k∈J∆

k 
X |Z∆
2R maxl∈Jn {|λl |}
jlk ||sL |
≤
+
1
|eTk wn | |eTk wn | − 2R maxl∈Jn {|λl |}
k∈J
∆

≤
(b)

≤

k
|Z∆
jlk ||sL |

X

|eTk wn |(1 − 2R maxl∈Jn {|λl |}/|eTk wn |)

X

|eTk wn |(1 − 2Rλ̄/ρ)

k∈J∆

k
|Z∆
jlk ||sL |

k∈J∆

,

(51)

where n = Node[k], index lk is defined in (16g), inequality
(a) is due to (45) and (46), and inequality (b) is due to (48).
The upper bound for the third term in (49), is given by
X

iI (v) =

k∈JY

X

k∈JY

Zjk

X
vk k
iL ≤
|Zjk ||ikL |.
|vk |

(52)

k∈JY

Lemma 2 is invoked to provide an upper bound for the
fourth term in (49):
X

Zjk iI (v)k =

X

Zjk iI (v)k ≤

k∈J∆

(49)

+

X

k∈J∆

k∈JY

Zjk iI (v) +

X

k∈JY

X
|sk |
|skL ||wk − vk |
|Zjk | L
+
|vk ||wk |
|wk |

∆

(a)

|eTk wn |(1 − 2R max{|λl |}/|eTk wn |)

|T(u)j −

|Zjk |

where (a) is due to (41), and (b) is due to (42) and (43).
For the second term in (49), we use Lemma 1, and we write
 k ∗
X
X
sL
∆
Zjk iPQ (v)k =
Zjlk
Tv
e
k n
k∈J∆
k∈J∆
"



 k ∗ #
∗
∗
X
skL
sL
skL
∆
−
+
Zjlk
=
T
T
T
ek vn
e k wn
ek vn
k∈J

(47)

Equation (47) provides a lower bound for the magnitude of
line to line voltages in terms of R and w.
The proof of (C2) is based on (47). In (47), in the term
in parenthesis, we replace the term |eTk wn | by its minimum
ρ [given in (16a)], and replace |λk | by its maximum λ̄ [also
given in (16a)] to obtain a lower bound for the line to line
voltages which is common for all k ∈ J∆ and n = Node[k]:
≥ |eTk wn |(1 − 2Rλ̄/ρ).

X

k∈JY

⇒ |eTk (vn − wn )| ≤ 2 max{|vl − wl |}

Zjk (

k∈JY

(a)

− wn )|

X

Zjk iPQ (v)k =

⇒

k∈J∆

X

Z∆
jlk

k∈J∆

X

k∈J∆

eTk vn k
i
|eTk vn | L

k
|Z∆
jlk ||iL |.

(53)

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

11

Using equations (50), (51), (52), and (53) in (49) yields
"
X
|Zjk ||skL |
1
wj
|≤
|T(u)j −
λj
|λj |
|wk |(1 − Rλ̄/w)
k∈J

≤
≤

X

+

k∈J∆

X

+

k∈JY

k
|Z∆
jlk ||sL |

|eTk wNode[k] |(1 − 2Rλ̄/ρ)
|Zjk ||ikL | +

X

k∈J∆

#

k
|Z∆
jlk ||iL | .

≤
(54)

The self-mapping property holds if the right-hand side of (54)
is upper-bounded by R for all j ∈ J , which is ensured by
(C3).
D. Contraction property
To obtain a sufficient condition for the contraction property,
we seek to find an upper-bound for the term kT(u)−T(ũ)k∞
proportional to ku − ũk∞ . Setting v = Λu and ṽ = Λũ, the
j-th entry of T(u) − T(ũ) is upper-bounded as follows:
|T(u)j − T(ũ)j |
1
1
=
Zj• [iPQ (v) + iI (v)] −
Zj• [iPQ (ṽ) + iI (ṽ)]
λj
λj
X
1
Zjk [iPQ (v)k − iPQ (ṽ)k ]
≤
|λj |
X

Zjk [iPQ (v)k − iPQ (ṽ)k ]

X

Zjk [iI (v)k − iI (ṽ)k ]

X

Zjk [iI (v)k − iI (ṽ)k ]

k∈J∆

+

k∈JY

+

k∈J∆

!

k∈JY

(c)

≤

(b)

≤

X

k∈JY

k∈JY

|λk uk − λk ũk |
(|wk | − R|λk |)2
|wk

− ũk∞ ,

(56)

where inequality (a) is due to (41) and (b) is due to (43).
Using Lemma 1 for the second term in (55) yields:
X

Zjk [iPQ (v)k − iPQ (ṽ)k ]

X

Z∆
jlk

k∈J∆

=

k∈J∆

X

|eTk wn |2 (1 − 2Rλ̄/ρ)2

k
2|Z∆
jlk ||sL |

k∈J∆

X

Zjk [iI (v)k − iI (ṽ)k ]

X

Zjk

max{|λl |}kun − ũn k∞
l∈Jn

max{|λl |}ku − ũk∞ , (57)
l∈Jn



vk k
ṽk k
i −
i
|vk | L |ṽk | L



X 2|Zjk ||ik |
(b) X
2|Zjk ||ikL |
L
|vk − ṽk | ≤
|vk − ṽk |
|vk |
|wk | − R|λk |
k∈JY

k∈JY

X

k∈JY

(55)

|Zjk ||skL |

|Zjk ||skL |
|λk |ku
|2 (1 − Rλ̄/w)2

|eTk wn |2 (1 − 2Rλ̄/ρ)2

− 2R maxl∈Jn {|λl |}/|eTk wn |)2

k
2|Z∆
jlk ||sL |

k∈J∆

k∈JY

≤

 k

X
X |Zjk ||sk |
sL ∗
skL ∗
L
Zjk ( ) − ( )
=
≤
|vk − ṽk |
vk
ṽk
|vk ||ṽk |
≤

X

k
2|Z∆
jlk ||sL |kvn − ṽn k∞

k∈JY

k∈JY

(a)

|eTk wn |2 (1

where n = Node[k], (a) is due to Hölder’s inequality and
(46), and (b) is due to (48).
For the third term in (55) we have that:

(a)

The expression in (55) is a sum of four terms. In what
follows, we find an upper bound for each term.
For the first term we have that
X
Zjk [iPQ (v)k − iPQ (ṽ)k ]

X

X

2|Zjk ||ikL |

|wk |(1 − Rλ̄/w)

|λk |ku − ũk∞ ,

X

Zjk [iI (v)k − iI (ṽ)k ]

X

Z∆
jlk

k∈J∆

=

k∈J∆
(a)

≤
(b)

≤
(c)

≤
≤



eTk vn k
eT ṽn k
iL − Tk
i
T
|ek vn |
|ek ṽn | L

k
X 2|Z∆
jlk ||iL |

k∈J∆

|eTk vn |

skL
eTk vn

∗

−



skL
eTk ṽn

∗ #



|eTk vn − eTk ṽn |

k
X 2|Z∆
jlk ||iL |kek k1 kvn − ṽn k∞

k∈J∆

|eTk wn | − 2R maxl∈Jn {|λl |}

X

k
4|Z∆
jlk ||iL |
max{|λl |}kun
T
|ek wn |(1 − 2Rλ̄/ρ) l∈Jn

X

|eTk wn |(1 − 2Rλ̄/ρ)

k∈J∆

k∈J∆

k
4|Z∆
jlk ||iL |

− ũn k∞

max{|λl |}ku − ũk∞ , (59)
l∈Jn

where n = Node[k], (a) is due to Lemma 3, (b) is due to (46)
and Hölder’s inequality, and (c) is due to (48).
Using (56), (57), (58), and (59) in (55) yields:
"
X
|Zjk ||skL ||λk |
1
|T(u)j − T(ũ)j | ≤
|λj |
|wk |2 (1 − Rλ̄/w)2
k∈J
Y

"

(58)

where inequality (a) is due to Lemma 3, (b) is due to (41), and
(c) is due to (43) and the fact that kuk − ũk k ≤ ku − ũk∞ .
For the fourth term, Lemma 2 is used to obtain:

.

k∈JY

(|eTk wn | − 2R maxl∈Jn {|λl |})2

|eTk (vn − ṽn )|

k
|Z∆
jlk ||sL |kek k1 kvn − ṽn k∞

k∈J∆

=

k∈JY

+

≤

X

k∈J∆

(b)

≤

|(eTk vn )||(eTk ṽn )|

k∈J∆

(a)

Y

k
|Z∆
jlk ||sL |

X

+

k
X 2|Z∆
jlk ||sL | maxl∈JNode[k] {|λl |}

k∈J∆

|eTk wNode[k] |2 (1 − 2Rλ̄/ρ)2

IEEE TRANSACTIONS ON POWER SYSTEMS (ACCEPTED)

+

X

k∈JY

+

2|Zjk ||ikL ||λk |
|wk |(1 − Rλ̄/w)

k
X 4|Z∆
jlk ||iL | maxl∈JNode[k] {|λl |}

k∈J∆

12

|eTk wNode[k] |(1 − 2Rλ̄/ρ)

#

Mohammadhafez Bazrafshan received the B.S. degree from Iran University
of Science and Technology in 2012, and the M.S. degree in 2014 from the
University of Texas at San Antonio both in electrical engineering, where he
is currently a PhD student with research interests in smart power grids.

ku − ũk∞

(60)

For T(u) to be a contraction, the coefficient multiplying ku −
ũk∞ must be less than 1 for all j ∈ J , which is ensured by
(C4).

R EFERENCES
[1] T. H. Chen, M. S. Chen, K. J. Hwang, P. Kotas, and E. A. Chebli,
“Distribution system power flow analysis-a rigid approach,” IEEE Trans.
Power Del., vol. 6, no. 3, pp. 1146–1152, July 1991.
[2] B. C. Lesieutre, P. W. Sauer, and M. A. Pai, “Existence of solutions
for the network/load equations in power systems,” IEEE Trans. Circuits
Syst. I: Fundam. Theory Appl., vol. 46, no. 8, pp. 1003–1011, Aug 1999.
[3] D. K. Molzahn, B. C. Lesieutre, and C. L. DeMarco, “A sufficient
condition for power flow insolvability with applications to voltage
stability margins,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2592–
2601, Aug. 2013.
[4] F. Wu and S. Kumagai, “Steady-state security regions of power systems,”
IEEE Trans. Circuits Syst., vol. 29, no. 11, pp. 703–711, Nov 1982.
[5] S. Bolognani and S. Zampieri, “On the existence and linear approximation of the power flow solution in power distribution networks,” IEEE
Trans. Power Syst., vol. 31, no. 1, pp. 163–172, Jan. 2016.
[6] H. D. Chiang and M. E. Baran, “On the existence and uniqueness of
load flow solution for radial distribution power networks,” IEEE Trans.
Circuits Syst., vol. 37, no. 3, pp. 410–416, Mar 1990.
[7] S. Yu, H. D. Nguyen, and K. S. Turitsyn, “Simple certificate of
solvability of power flow equations for distribution systems,” in Proc.
IEEE PES General Meeting, July 2015, pp. 1–5.
[8] C. Wang, A. Bernstein, J. Y. L. Boudec, and M. Paolone, “Explicit conditions on existence and uniqueness of load-flow solutions in distribution
networks,” IEEE Trans. Smart Grid, 2016.
[9] T. Q. Zhao, H. D. Chiang, and K. Koyanagi, “Convergence analysis of
implicit Z-bus power flow method for general distribution networks with
distributed generators,” IET Generation, Transmission & Distribution,
vol. 10, no. 2, pp. 412–420, Feb. 2016.
[10] “Load representation for dynamic performance analysis [of power systems],” IEEE Trans. Power Syst., vol. 8, no. 2, pp. 472–482, May 1993.
[11] W. H. Kersting, Distribution System Modeling and Analysis, 3rd ed.
CRC Press, 2002.
[12] A. Garces, “A linear three-phase load flow for power distribution
systems,” IEEE Trans. Power Syst., vol. 31, no. 1, pp. 827–828, Jan.
2016.
[13] T. H. Chen, M. S. Chen, T. Inoue, P. Kotas, and E. A. Chebli, “Threephase cogenerator and transformer models for distribution system analysis,” IEEE Trans. Power Del., vol. 6, no. 4, pp. 1671–1681, Oct. 1991.
[14] J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear
Equations in Several Variables. Academic Press, Inc., 1970.
[15] D. Bertsekas and J. N. Tsitsiklis, Parallel and Distributed Computation:
Numerical Methods. Massachussetts, USA: Athena Scientific, 1989.
[16] J. Nocedal and S. J. Wright, Numerical Optimization. New York, USA:
Springer, 2006.
[17] R. W. D. Nickalls, “The quartic equation: Invariants and euler’s solution
revealed,” The Mathematical Gazette, vol. 93, no. 526, pp. 66–75,
2009. [Online]. Available: http://www.jstor.org/stable/40378672
[18] M. Bazrafshan and N. Gatsis, “On the solution of the three-phase loadflow in distribution networks,” in 50th Asilomar Conference on Signals,
Systems, and Computers, Nov. 2016.
[19] W. H. Kersting, “Radial distribution test feeders,” in Proc. IEEE PES
Winter Meeting, vol. 2, 2001, pp. 908–912 vol.2.
[20] “Distribution Test Feeders,” https : //ewh.ieee.org/soc/pes/dsacom/
testfeeders/, accessed: 2016-10-12.
[21] M. J. Gorman and J. J. Grainger, “Transformer modelling for distribution
system studies part II: Addition of models to YBUS and ZBUS,” IEEE
Trans. Power Del., vol. 7, no. 2, pp. 575–580, April 1992.

Nikolaos Gatsis received the Diploma degree in Electrical and Computer
Engineering from the University of Patras, Greece, in 2005 with honors.
He completed his graduate studies at the University of Minnesota, where
he received the M.Sc. degree in Electrical Engineering in 2010, and the
Ph.D. degree in Electrical Engineering with minor in Mathematics in 2012.
He is currently an Assistant Professor with the Department of Electrical and
Computer Engineering at the University of Texas at San Antonio. His research
interests lie in the areas of smart power grids, renewable energy management,
communication networks, and cyberphysical systems, with an emphasis on
optimal resource management. Prof. Gatsis co-organized symposia in the
area of Smart Grids in IEEE GlobalSIP 2015 and IEEE GlobalSIP 2016.
He also served as a Technical Program Committee member for symposia
in IEEE SmartGridComm 2013, IEEE SmartGridComm 2014, and IEEE
SmartGridComm 2015.

