NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

Multiphase Optimal and NonSingular Power Flow by

Successive Linear

Approximations

Preprint

Andrey Bernstein 

National Renewable Energy Laboratory

Cong Wang and Jean-Yves Le Boudec 

Ecole Polytechnique F'ed'erale de Lausanne (EPFL)

Presented at the IEEE Power Systems Computation Conference 

(PSCC)

Dublin, Ireland

June 11-15, 2018

© 2018 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

Suggested Citation

Bernsetin, Andrey, Cong Wang, and Jean-Yves Le Boudec. 2018. 

“Multiphase Optimal and Non-Singular Power Flow by Successive Linear

Approximations: Preprint.” Golden, CO: National Renewable Energy

Laboratory. NREL/CP-5D00-70368 

https://www.nrel.gov/docs/fy18osti/70368.pdf

Conference Paper

NREL/CP-5D00-70368

August 2018
NOTICE

This work was authored by the National Renewable Energy Laboratory, operated by Alliance for Sustainable 

Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding 

provided by U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Solar Energy 

Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the 

U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, 

acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to 

publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

This report is available at no cost from the National Renewable 

Energy Laboratory (NREL) at www.nrel.gov/publications.

U.S. Department of Energy (DOE) reports produced after 1991 

and a growing number of pre-1991 documents are available 

free via www.OSTI.gov.

Cover Photos by Dennis Schroeder: (left to right) NREL 26173, NREL 18302, NREL 19758, NREL 29642, NREL 19795.

NREL prints on paper that contains recycled content.
Multiphase Optimal and Non-Singular Power Flow

by Successive Linear Approximations

Andrey Bernstein, Member, IEEE

National Renewable Energy Laboratory (NREL)

Golden (CO), USA

andrey.bernstein@nrel.gov

Cong Wang, Student Member, IEEE

Jean-Yves Le Boudec, Fellow, IEEE

Ecole Polytechnique F  ́ ed ́ erale de Lausanne (EPFL)  ́

Lausanne (VD), Switzerland

{cong.wang, jean-yves.leboudec}@epfl.ch

Abstract—We propose an optimal power flow method for

multiphase power systems. The method works for both radial

and meshed networks, and is able to accommodate different types

of load/source connections such as wye, delta, and combination

thereof. Moreover, the method generates a sequence of nonsingular points (i.e., associated with non-singular power-flow

Jacobian) in the electrical state space and a corresponding

sequence of power injections, which can be used in real-time

control settings. Under certain conditions, the sequences attain

the local minimum. The method has been tested using several

typical networks, including the IEEE 37-bus and 123-bus test

feeders.

Index Terms—Multiphase power networks; AC optimal power

flow; successive linear approximation; non-singularity; existence;

uniqueness; feasibility.

I. INTRODUCTION

The AC optimal power flow (OPF) problem is one of the

fundamental problems in power system operation and analysis.

Mainly due to the nonlinearity of the power flow equations, the

OPF problem is nonconvex and NP-hard in general (see, e.g.,

[1]). In transmission systems, OPF methods are typically based

on DC power-flow models and on the assumption that the

system is balanced. In distribution systems, however, these two

simplifying assumptions are no longer valid, since distribution

systems generally have high R/X ratio and are unbalanced with

a variety of different types of connections. Therefore a full

multiphase power-flow model has to be used, and AC OPF

methodology has to be developed and applied.

There is a wide range of literature that addressed the

AC OPF problem by using relaxation techniques, such as

semidefinite relaxation (SDR) [2], [3] and generalizations

thereof (e.g., [4]); due to the scope of this paper, we will not

present an exhaustive list here. Several papers addressed the

The author from NREL (1st) and the authors from EPFL (2nd and 3rd)

have equal contribution in the paper. The work of A. Bernstein, was supported

by the U.S. Department of Energy under Contract No. DE-AC36-08GO28308

with the National Renewable Energy Laboratory. The U.S. Government retains

and the publisher, by accepting the article for publication, acknowledges that

the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide

license to publish or reproduce the published form of this work, or allow

others to do so, for U.S. Government purposes. The work of C. Wang and

J.-Y. Le Boudec was supported by the Swiss National Foundation SNSF under

grant NRP70 - “Energy Turnaround”.

AC OPF problem for multiphase distribution networks (see [5]

for an extensive review). A majority of them (see, e.g., [6]–

[8] and pertinent references therein) utilized generic solvers

to identify a solution of this challenging nonconvex (and

NP-hard) optimization problem; nonetheless, no theoretical

convergence guarantees are offered. In [9], an SDR approach

was proposed, and a feasible (and globally optimal) solution

can be identified under certain conditions. Similarly, in [10], an

SDR approach was applied for multiphase radial distribution

networks. Finally, [11] developed a successive convex approximation methodology, which is proven to converge to a KKT

point of the original non-convex AC OPF. Overall, with the

exception of [7], the network models utilized in the existing

literature can support only wye-connected sources, and they

cannot be straightforwardly extended to delta connections.

In this paper, we propose to solve the AC OPF for multiphase unbalanced systems using a new solution methodology

based on successive linear approximations. The distinctive

characteristics of the proposed method are highlighted as

following: 1) We leverage the general multiphase power-flow

model recently proposed in [12]. This allows us to explicitly

consider different types of load/source connections, such as

wye, ungrounded delta, and a combination thereof; 2) We

propose the following iterative method. At each iteration:

(1) We convexify the original non-convex problem by replacing the exact AC power-flow constraints with an

appropriate linear approximation, leveraging the model

in [12]. The linearization is performed around the powerflow solution obtained in the previous iteration.

(2) In this convex optimization problem, we explicitly impose

constraints on the power injections to guarantee the

existence of the exact AC power-flow solution that is

non-singular (namely, having a non-singular power-flow

Jacobian) and unique in an analytically specified domain.

(3) After obtaining the optimal power injections in the convexified problem, we compute its exact AC power-flow

solution using the algorithm in [12].

a) If the exact AC power-flow solution is feasible in the

sense that it truly satisfies the security constraints, then

we set this solution as the next linearization point and

move on to the next iteration.

b) Otherwise, we shrink the constraints imposed by The1 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
orem 2 in [12], recompute the solution and repeat this

step until the conditions in (a) are satisfied.

We prove that the algorithm generates a non-singular sequence of power-flow solutions that, under certain conditions,

approaches a local minimum in the original multiphase AC

OPF problem. Moreover, all the members in this sequence

belong to the same non-singular path-connected component in

the voltage space, which means that we can move from one

member to another without passing the steady-state stability

limit.

Overall, the main advantages of the proposed algorithm are:

(i) It generates a sequence of non-singular points in the state

space and a corresponding sequence of power injections;

the latter can be used in real-time control settings to

dispatch setpoints.

(ii) It is applicable to the general multiphase model and

accommodates different types of connections, including

wye, delta, and a combination thereof.

(iii) It does not require radiality of the network topology.

It is worth mentioning that feature (i) above is important in

practice as it ensures that the entire trajectory from the current

operating point to the optimal one satisfies the operational

constraints and does not pass through the stability limits.

Recently, there were several works on successive linear

or convex approximations for solving OPF problems (see,

e.g., [11], [13], [14]). In particular, in [13], [14], a sequential

linearization with consecutive power-flow verification (of the

similar type proposed in our paper) was used to solve a

certain planning problem for transmission networks; however,

the approach was tailored to balanced networks (hence, no

wye/delta connections were considered), and no convergence

results were provided.

The paper is structured as follows. In Section II, we present

the multiphase distribution network model. In Section III, we

formulate the AC OPF problem and its convex relaxation using

linearized power-flow model. In Section IV, we present our

algorithm and main results. In Section V, we evaluate the

performance of the algorithm on several typical test feeders.

Finally, we present some concluding remarks in Section VI.

Notation. Upper-case (resp. lower-case) boldface letters are

used for matrices (resp. column vectors); (·)

T

for transposition;

and | · | for the absolute value of a number or the componentwise absolute value of a vector or a matrix. For a complex

number c ∈ C, <{c} and ={c} denote its real and imaginary

part, respectively; and c denotes the conjugate of c. For an N ×

1 vector x ∈ C

N , kxk∞ := max(|x1|...|xn|), kxk1 := PN

i=1 |xi|,

and diag(x) returns an N × N matrix with the entries of x in

its diagonal. For an M ×N matrix A ∈ C

M×N , the `∞-induced

norm is defined as ||A||∞ = maxi=1,...,M PN

j=1 |(A)ij |.

II. SYSTEM MODEL

In this section, we briefly overview the general multiphase

power-flow model recently proposed and analyzed in [12]. A

simpler version, without delta loads and with only fully threephase scenarios, can be found in the two-page letter [15].

A. Power-Flow Model

Consider a generic three-phase distribution network with

one slack bus and N three-phase P Q buses. Each bus in the

network can have different types of load/source connections,

such as wye, ungrounded delta, and a combination thereof.

This model is compactly described by the following set of

non-linear equations:

s

Y = diag(v)



YL0v0 + YLLv



− diag H

T

i∆



v, (1a)

s

∆ = diag (Hv)i∆. (1b)

In (1), s

Y

, s

∆ ∈ C3N are the complex power injection vectors

for wye and delta connections, respectively; v ∈ C

3N collects

the complex voltages at every bus and phase of the P Q buses;

v0 ∈ C

3

collects the complex voltages at the slack bus; i

∆ ∈

C

3N collects the phase-to-phase currents at every bus; Y00 ∈

C

3×3

, YL0 ∈ C

3N×3

, Y0L ∈ C

3×3N , and YLL ∈ C3N×3N are

the submatrices of the three-phase admittance matrix

Y := 

Y00 Y0L

YL0 YLL

∈ C

3(N+1)×3(N+1)

, (2)

which can be formed from the topology of the network, the

π-model of the transmission lines, and other passive network

devices, as shown in, e.g., [16]; and H is a 3N × 3N blockdiagonal matrix defined by

H :=







Γ

.

.

.

Γ





 , Γ :=





1 −1 0

0 1 −1

−1 0 1



 . (3)

Note that (1) defines an explicit mapping from the state

vector (v

T

,(i

∆)T

)

T ∈ C6N to the vector of power injections

s := ((s

Y

)

T

,(s

∆)T

)

T ∈ C6N . For further development, we

sometimes find it useful to represent the power-flow equations in real (rather than complex) space. To this end, let

x := <{s

Y

}

T

, ={s

Y

}

T

, <{s

∆}T

, ={s

∆}T

T

collect the active

and reactive power injections of wye and delta sources. Also,

let y := <{v}

T

, ={v}

T

, <{i

∆}T

, ={i

∆}T

T

denote the vector

of the state variables. Then, the power-flow equations can be

written as

x = h(y), (4)

where h : R

12N → R12N is the mapping defined explicitly

by (1). When need to revert the real-valued expressions back

to the complex-valued counterparts, we define the operator

comp() such that comps(x) and compv(y) represent the s and

v that correspond to x and y, respectively.

B. Conditions for Existence, Uniqueness, and Non-Singularity

In [12], explicit conditions for existence, uniqueness, and nonsingularity of the solution to the multiphase power-flow equations

were derived. We next summarize these conditions, as they will be

used explicitly in our proposed AC OPF algorithm in Section IV;

see [12] for a thorough review of other similar conditions recently

proposed in the literature.

Let w ∈ C

3N denote the zero-load voltage, which is given

explicitly by w := −Y−1

LLYL0v0. Also, let W := diag(w), and

L := |H| be the component-wise absolute value of the matrix H

defined in (3). For s := ((s

Y

)

T

,(s

∆)T

)

T ∈ C6N define the following

norm ξ(·) on C

6N (see Lemma 1 in [12]):

ξ

Y

(s) :=










W−1Y

−1

LLW−1

diag(s

Y

)










∞

, (5a)

ξ

∆(s) :=










W−1Y

−1

LLH

T

diag(L|w|)

−1

diag(s

∆)










∞

, (5b)

ξ(s) := ξ

Y

(s) + ξ

∆(s), (5c)

2 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
where |w| is the component-wise absolute value of the vector w, and

kAk∞ is the induced `∞-norm of a complex matrix A.

Next, for v ∈ C

3N , let

α(v) := min

j

|(v)j |

|(w)j |

(6a)

β(v) := min

j

|(Hv)j |

(L|w|)j

(6b)

γ(v) := min {α(v), β(v)} . (6c)

Finally, for any ρ ≥ 0, define

Dρ(v) := {v

0

: |(v

0

)j − (v)j | ≤ ρ|(w)j |, j = 1 . . . 3N}. (7)

We next wrap up the main results from [12] for completeness.

Theorem 1 (Theorems 2 and 3 in [12]). Let vb be a given solution

to the power-flow equations with power injection bs satisfying:

ξ(bs) < (γ(vb))2, (8)

where ξ(·) and γ(·) are given in (5) and (6), respectively. Consider

some other candidate power injections vector s, and assume that

ξ(s − bs) <

1

4



(γ(vb))2 − ξ(bs)

γ(vb)

2

. (9)

Let

ρ

‡

(vb,bs) := 1

2



(γ(vb))2 − ξ(bs)

γ(vb)



(10a)

ρ

†

(vb,bs, s) := ρ

‡

(vb,bs) −

q

(ρ

‡(vb,bs))2 − ξ(s − bs) (10b)

Then:

(i) The power-flow solution (vb,bs) is non-singular, in the sense that

the Jacobian matrix of the mapping h defined in (4) evaluated

at vb,bi

∆ := diag−1



Hvb



bs

∆



is invertible;

(ii) There exists a unique power-flow solution v in Dρ(vb) defined

in (7) with ρ = ρ

‡

(vb,bs);

(iii) The solution is located in Dρ(vb) with ρ = ρ

†

(vb,bs, s);

(iv) The solution v satisfies ξ(s) < (γ(v))2, hence is non-singular

as well;

(v) The solution v can be found by iteration (11) starting from

anywhere in Dρ(vb) with ρ = ρ

‡

(vb,bs):

v ← w + Y

−1

LL 

diag(v)

−1

s

Y + HT

diag (Hv)

−1

s

∆



. (11)

Once v is obtained, i

∆ can be uniquely recovered using (1).

Remark 1. For general multiphase networks, the vectors v, s

Y

, and

w collect their corresponding electrical quantities for the existent

phases; the vectors i

∆ and s∆ collect the electrical quantities for the

existent phase-to-phase connections; H contains rows that correspond

to the existent phase-to-phase connections. More precisely, H is an

N

∆ × Nphase matrix, where N∆ is the total number of phase-tophase connections and N

phase is the total number of phases. When

there are no phase-to-phase connections in the network, everything

still holds after removing the terms that involve H, L, i

∆ and s∆.

III. OPF PROBLEM FORMULATION

We consider the following prototypical multiphase AC OPF problem

(P0) min

x∈R

12N ,y∈C12N

f(x) (12a)

subject to : x ∈ X (12b)

x = h(y) (12c)

V

min ≤ |compv

(y)| ≤ V

max

, (12d)

where f(·) is a convex function, X ⊆ R

12N is a convex compact set,

h(·) is the multiphase power-flow mapping, and |compv(y)| ∈ R

3N

is the vector of voltage magnitudes that is decided by y. (Recall the

definition of comp() at the end of Section II-A.)

As it is well-known, (P0) is a non-convex problem due to the

non-linear equality constraints (12c) and the lower bound constraint

in (12d). Also, from the practical point of view, there might be several

solutions to the power-flow equations (12c), some of which might be

singular1.

To tackle these challenges, we propose a successive linear approximation of the constraints (12c) and (12d). To this end, let xb be a

feasible power injection vector with a voltage vb that satisfy condition

(8). Consider now a linear approximation of the voltage magnitudes

as a function of power injections:

|v| ≈ K(vb,xb)x + a(vb,xb) (13)

with the property that K(vb,xb)xb + a(vb,xb) = |vb|. This linear approximation can be efficiently computed, for example, by computing

the inverse of the load-flow Jacobian, or any other method; see,

e.g., [12], [18], [19]. In this paper, we assume that the matrix

K(vb,xb) ∈ R

12N×3N and vector a(vb,xb) ∈ R3N are continuous in

(vb, xb). These assumptions hold for the approximations proposed in,

e.g., [12]. Let



max

(vb,xb)

:=

1

4



(γ(vb))2 − ξ(bs)

γ(vb)

2

, (14)

and for any  ∈ (0, max

(vb,xb)

), consider the following local optimization

problem:

(P1(vb, xb, )) min

x∈R

12N

f(x) (15a)

subject to : x ∈ X (15b)

V

min ≤ K(vb,xb)x + a(vb,xb) ≤ Vmax (15c)

ξ (comps(x − xb)) ≤ , (15d)

where ξ(·) is defined in (5)2. Observe that (P1(vb, xb, )) is a convex

optimization problem, and constraint (15d) guarantees existence,

uniqueness, and non-singularity of the exact power-flow solution

as prescribed by Theorem 1. An iterative algorithm based on

(P1(vb, xb, )) is proposed in the next section.

Remark 2. Other constraints can be added to (P0) and (P1). For

example, if we plan to consider the ampacity of branch currents, we

could approximate the branch current by affine functions of x that

are similar to (13). Indeed, this can be efficiently conducted, since the

branch currents are linear functions of the nodal voltages in generic

two-port device modeling. However, for the ease of exposition, we

prefer to analyze here a simpler problem formulation.

IV. ALGORITHM AND MAIN RESULT

In this section, we propose a successive convex approximation

algorithm that is based on solving a sequence of (P1(vb, xb, )) until

convergence. In order to initialize the algorithm, a strictly feasible

solution (v

(0)

, x

(0)) to the power-flow problem is required, in the

sense that

x

(0) ∈ X , V min < |v(0)| < V max

. (16)

In addition, it is required that this solution satisfies (8). Practically,

such initial values can be obtained from state estimation procedures.

However, in the cases where state estimation is not available, one can

either simply choose v

(0) = w and x(0) = 0 as w is almost always

strictly feasible, or utilize techniques as in [11] to find a feasible

point.

The algorithm is given in Algorithm 1.

1See [17] for a practical example where singular state can be feasible.

2Note that (5) defines a norm on C6N which is equivalent to the norm on

R12N .

3 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
Algorithm 1 Multiphase AC OPF

Input: Strictly feasible power-flow solution (v

(0)

, x

(0)), satisfying (16) and (8).

Input: A parameter β ∈ (0, 1).

Output: (v

∗

, x

∗

), e, {v

(k)

, x

(k)

, 

(k)

seq}

1: k ← 0, flagA← 1

2: while flagA do

3: Compute the linear model (K(v(k),x(k)), a(v(k),x(k)))

using (v

(k)

, x

(k)

).

4: n ← 0, flagB← 1

5: while flagB do

6: Set 

(n)

:= β

n+1max

(v(k),x(k))

, where 

max

(v(k),x(k))

is

given in (14).

7: {SOLVE THE LOCAL OPTIMIZATION PROBLEM}

8: Solve (P1(v

(k)

, x

(k)

, (n))), obtain the optimal solution xe.

9: {SOLVE THE POWER-FLOW EQUATIONS FOR

THE EXACT SOLUTION}

10: Solve the power-flow equations for xe using the

method in [12], let ve denote the solution,

11: {CHECK FEASIBILITY}

12: if ve is strictly feasible then

13: v

(k+1) ← ve, x(k+1) ← xe, 

(k)

seq ← 

(n)

, flagB← 0

14: else

15: n ← n + 1

16: end if

17: end while

18: if



f(x

(k+1)) − f(x(k)

)



 < ErrBound then

19: flagA← 0, e ←





x

(k+1) − x(k)







2

20: end if

21: k ← k + 1

22: end while

23: v

∗ ← v(k)

, x

∗ ← x(k)

24: return (v

∗

, x

∗

), e, {v

(k)

, x

(k)

, 

(k)

seq}

Remark 3. In fact, each obtained (v

(k)

, x

(k)

) has an equivalent

(x

(k)

, y

(k)

). The latter will be used in the statement of the main

theorem.

Remark 4. Some of the practical scenarios may require that the

objective function explicitly involves the electrical state y. If that is

the case, then we want to minimize f(x, y). To extend the proposed

method to these scenarios, we have to make following modifications:

• In (P1(vb, xb, )), an extra constraint y = K0

(vb,xb)x + a

0

(vb,xb)

should be added, where K0

(vb,xb)

and a

0

(vb,xb)

can be obtained in

the way given by [12];

• In line 12 of the algorithm, the condition should be replaced by

“If ve is strictly feasible and f(xe, ye) ≤ f(x

(k)

, y

(k)

)”.

For brevity, in this paper, we focus on the cases where f(·) does not

explicitly involve y.

Below is our convergence result regarding Algorithm 1. The proof

is deferred to the Appendix.

Theorem 2. The following statements hold for Algorithm 1:

(i) The sequence of solutions (x

(k)

, y

(k)

) is

• feasible for problem (P0);

• non-singular, in the sense that the Jacobian matrix of h(·)

evaluated at y

(k)

is invertible;

• connected to each other via non-singular paths.

(ii) f(x

(k)

) ≤ f(x

(k−1)) for all k.

(iii) When ErrBound > 0, the algorithm converges (stops) in a

finite number of steps. Upon convergence, if e = 0, then

(x

∗

, y

∗

) is a local minimum of (P0).

(iv) In addition, when ErrBound is set to 0, we obtain an infinite

sequences such that:

(a) If

min := lim inf

k→∞



(k)

seq > 0, (17)

we have that any limit point of the sequence of solutions

{x

(k)

, y

(k)

} produced by Algorithm 1 is a non-singular

local minimum of (P0).

(b) If

max := lim sup

k→∞



(k)

seq > 0, (18)

then there exists a limit point of the sequence {x

(k)

, y

(k)

}

that is a non-singular local minimum of (P0).

Remark 5. We note that Theorem 2 does not establish convergence

of the sequence of solutions to a local minimum in general. Rather,

under condition (17), it is guaranteed that any limit point is a local

minimum, while under a milder condition (18), we know that there

exists such a limit point. In case the sequence indeed converges,

both part (iv) (a) and (b) state that it converges to a local minimum.

Conditions under which convergence in solutions is guaranteed is left

for future work.

Remark 6. Conditions (17) and (18) are a-posteriori conditions

that cannot be verified before the algorithm’s run. As we show in

the numerical examples in Section V, these conditions are typically

satisfied in practice.

V. PERFORMANCE EVALUATION

In this section, we first illustrate the performance of the algorithm

using a simple example, which involves a single-phase two-bus

network3and a quadratic objective function. Then, we adopt the

IEEE 37-bus and 123-bus test feeders to evaluate performance of

the proposed method. In particular, (i) the regulators in the two

feeders are set to their nominal values; (ii) the objective function

is linear with randomly generated coefficients. In addition, note that

(15d) is essentially a collection of linear constraints imposed on

|s − bs|. Therefore, we introduce auxiliary variables and formulate

(P1(vb, xb, )) as a second-order conic programming (SOCP) problem,

which can be efficiently solved by the interior-point method.

Finally, we assume that (i) the slack bus voltages are known and

set to nominal values; (ii) the voltage bounds are 0.9 and 1.1 p.u.;

(iii) ErrBound is 0.0001 in all numerical analyses unless otherwise

specified.

A. Illustrative Example

Here, we consider a single-phase two-bus network, where the

P Q bus is directly attached to the slack bus via a line with series

impedance 0.3+ 0.7 p.u. Suppose that (i) there is no shunt element;

(ii) the P Q-bus power injection belongs to the triangular region

shown on Figure 1(b), which defines X ; (iii) the objective function

is kx − (3, 0)Tk

2

2. Note that, with this objective function, we try to

encourage the local (active) power generation.

In this network, it is easy to check that v

(0) = 1 and x(0) =

(0, 0)Tform a valid initial point. With this initial point and β = 0.9,

the algorithm converges in 3 steps (i.e., 4 instances of (P1(vb, xb, ))

are solved). We plot the trajectory of voltage v

(k)

in Figure 1(a) and

3

In fact, it is equivalent to a balanced three-phase network.

4 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
0.8 0.9 1 1.1

Real

-0.1

0

0.1

0.2

0.3

0.4

0.5

Imag

Voltage at the PQ Bus

0.9 p.u.

1.1 p.u.

(a)

Power Injection at the PQ Bus

6

6

6.5 6.5

7

7

7.5 7.5

8

8

8.5 8.5

9

9

9.5 9.5

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

Real

-0.2

-0.1

0

0.1

0.2

Imag

(b)

Figure 1: Trajectories in (a) voltage space and (b) power

injection space. For each trajectory of (a), the lowest point

represents v

(0); for each trajectory of (b), the leftmost point

stands for x

(0). Note that in voltage space, a point is nonsingular if its real part is not 0.5 p.u.

its corresponding trajectory x

(k)

in Figure 1(b) via thick solid lines.

Clearly, the obtained v

(k)

are all strictly feasible and x

∗ = (0.5, 0)T

is the unique minimum in X . Moreover, the obtained points are nonsingular and connected to each other via nonsingular paths. Further

consider that this network can be solved analytically, therefore we are

able to generate other valid initial values for v

(0) and x(0). Similar to

the case where v

(0) = 1 and x(0) = (0, 0)T

, we randomly generate

some v

(0)

, x

(0) and plot their corresponding trajectories in Figure 1

via dashed lines of different thickness. As can be seen, all of them

are feasible and nonsingular in voltage space and converge to the

same optimal value.

It is also worth noticing that, when set ErrBound to 0, we obtain

numerically lim infk→∞ 

(k)

seq = lim supk→∞ 

(k)

seq = 0.1211. This

means that the conditions in the fourth item of Theorem 2 hold.

B. IEEE 37-Bus Feeder

In the IEEE 37-bus feeder, all the loads are delta-connected. We

choose v

(0) = w, x(0) = 0, which is feasible for this network.

Now, for a positive real number κ, construct set X to be: X =

{x : κxbenchmark

j ≤ xj ≤ 0, j ∈ {1, ..., 12N}}, where xj is the jth entry of x and x

benchmark

j

is the benchmark value of xj (coming

from the feeder data) that is negative due to consumption.

First, we analyze how the required number of local problem

(P1(vb, xb, )) (i.e., the number of inner loops) depends on κ. In

other words, we want to test whether a larger set X results in

a much greater number of local problem (P1(vb, xb, )) to solve.

Indeed, this is important for the algorithm performance as the total

complexity includes the complexity of (P1(vb, xb, )) and the number

of (P1(vb, xb, )) instances solved. For a good algorithm, the number

of required (P1(vb, xb, )) should not change dramatically with the size

of X .

We fix β at 0.8 and gradually increase κ from 1 to 6. The

numerical results show that the algorithm finishes after solving (i)

1 local problem (P1(vb, xb, )) when κ = 1; (ii) 5 to 7 local problem

(P1(vb, xb, )) when κ = 6. Thus, the algorithm scales well with the

size of X .

Next, we fix the value of κ at 6, and check the impact of β

on the required number of local problem (P1(vb, xb, )) to solve. By

letting β be from {0.2, 0.4, 0.6, 0.8}, we observe that the number of

local problem (P1(vb, xb, )) decreases as β increases. In fact, this is

partially attributed to the accuracy of the linear approximation (13).

In particular, the fact that the linear model is accurate means that

the true values of the state variables are accurately represented by

their approximations. Therefore, if the approximations are feasible,

then the true values are also feasible in most of the cases. By such

observation, any small β becomes conservative and only leads to

more instances of (P1(vb, xb, )) to solve.

C. IEEE 123-Bus Feeder

In the IEEE 123-bus feeder, we add an additional delta-connected

power injection (−0.03 −0.01, −0.03 −0.01, −0.03 −0.01) p.u.

to bus 1, and hence create a bus with mixed wye-delta load. Denote

the modified benchmark power vector as x

benchmark and let (i) v(0)

be the solution to x

benchmark that is guaranteed in Dρ‡(w,0)

(w) by

Theorem 1; (ii) x

(0) = xbenchmark. Then, after repetition of the

analysis in the IEEE 37-bus feeder, we obtain similar results. In

particular, for κ = 6, the number of required (P1(vb, xb, )) is between

8 and 11. This proves that the algorithm is able to work on networks

that have sizes comparable to the IEEE 123-bus feeder.

However, we also notice a limitation of the performance. For

large networks, the number of constraints coming from (15d) is

very large. In order to reduce the complexity, we have to manually

remove some of the constraints that are redundant (i.e., the ones

that are implied by some other constraints). For radial networks with

negligible shunt elements, this can be done through the way given in

[20]. Nonetheless, how to reduce the complexity in general networks

still requires further study.

VI. CONCLUSION

In the paper, we propose an optimal power flow method for

multiphase power systems. The method does not require any specific

network topology and is able to accommodate different types of

load/source connections. Moreover, the method generates a sequence

of non-singular points in the state space and a corresponding sequence

of power injections, which can be used in real-time control settings.

Under certain conditions, the sequences attain the local minimum.

Despite the advantages of this method, there are several limitations

as we mentioned in the paper. The improvement is left as a future

work.

APPENDIX

AUXILIARY SETS AND LEMMAS

A. Auxiliary Sets

For any  > 0, define the open set

Ye(vb, xb, ) := nx ∈ X : V

min < K(vb,xb)x + a(vb,xb) < V max

,

ξ(comps(x − xb)) < o(19)

5 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
and set

Y(xb, ) := nx ∈ X : ∃y ∈ C

12N s.t. x = h(y),

V

min < |compv

(y)| < V max, ξ(comps(x − xb)) < o.

(20)

Moreover, let B(xb, ) := {x ∈ X : ξ(comps(x − xb)) < }. These

auxiliary sets form the core of following lemmas.

B. Lemmas

Lemma 1. Assume that xb ∈ X and the power flow solution vb

to xb satisfy (i) V

min < K(vb,xb)xb + a(vb,xb) < V max and (ii)

ξ(comps(xb)) < (γ(vb))2. Then, there exists 0(vb, xb) > 0 such that

for all  ∈ (0, 0(vb, xb)], we have that

B(xb, ) = Ye(vb, xb, ) = Y(xb, ). (21)

Proof. First, consider that:

1) K(vb,xb)x + a(vb,xb) is an affine function from R

12N to R3N ;

2) V

min < K(vb,xb)xb + a(vb,xb) < V max;

3) ξ(·) is a norm and defines a metric that is Lipschitz equivalent

to the Euclidean distance.

Therefore, there exists 1(vb, xb) > 0 such that

ξ(comps(x − xb)) < 1(vb, xb)

⇒V

min < K(vb,xb)x + a(vb,xb) < V max (22)

Namely, we have

Ye(vb, xb, ) = B(xb, ), ∀ ∈ (0, 1(vb, xb)]. (23)

Next, consider:

1) |vb| = K(vb,xb)xb + a(vb,xb) ⇒ V

min < |vb| < V max;

2) ξ(comps(xb)) < (γ(vb))2;

3) Based on item 2), if ξ(comps(x − xb)) < max

(vb,xb)

, then there is

a solution in Dρ† (vb);

4) The value of ρ

†

(vb, comps(xb), comps(x)) is controlled by

ξ(comps(x − xb)) as

ρ

†

(vb, comps(xb), comps(x))

=ρ

‡

(vb, comps(xb)) −

p

(ρ

‡(vb, comps

(xb)))2 − ξ(comps(x − xb))

=ρ

‡

(vb, comps(xb)) 1 −

s

1 −

ξ(comps(x − xb))

(ρ

‡(vb, comps

(xb)))2

!

≤

ξ(comps(x − xb))

ρ

‡(vb, comps

(xb)) . (24)

Therefore, there exists 2(vb, xb) ∈ (0, max

(vb,xb)

] such that

ξ(comps(x − xb)) < 2(vb, xb)

⇒∃y s.t. x = h(y), V min < |compv(y)| < V max (25)

Namely, we have

Y(xb, ) = B(xb, ), ∀ ∈ (0, 2(vb, xb)]. (26)

Let 0(vb, xb) = min{1(vb, xb), 2(vb, xb)}, we complete the proof.

Using Lemma 1, we are able to prove that at each iteration, we

either output a feasible solution with a smaller objective value, or

converge to a local minimum.

Lemma 2. Assume that the conditions of Lemma 1 hold. Let xe =

arg minx∈cl(Ye(vb,xb,)) f(x) and suppose that xe ∈ cl(Y(xb, )), for

some  ∈ (0, max

(vb,xb)

). Then xe = arg minx∈cl(Ye(vb,xb,))∩cl(Y(xb,)) f(x).

Moreover, if xe = xb, then xb is a local minimum of f(·) over

cl(Y(xb, ∞)).

Proof. First, observe that under the conditions of the lemma, we have

that

f(xe) ≤ f(x), ∀x ∈ cl(Ye(vb, xb, )) ∩ cl(Y(xb, )) ⊆ cl(Ye(vb, xb, ))

xe ∈ cl(Ye(vb, xb, )) ∩ cl(Y(xb, )).

Therefore, the first part of the lemma trivially follows.

Second, notice that by Lemma 1, we have that for any  > 0, there

exists δ > 0 such that

B(xb, δ) ⊆ Ye(vb, xb, )∩Y(xb, ) ⊆ cl(Ye(vb, xb, ))∩cl(Y(xb, )). (27)

Indeed, this is true for δ =  if  ≤ 0(vb, xb), and δ = 0(vb, xb)

otherwise. Therefore, by the definition of a local minimum, if

arg minx∈cl(Ye(vb,xb,))∩cl(Y(xb,)) f(x) = xb, then xb is a local minimum.

PROOF OF THEOREM 2

C. Properties of v

(k)

, x

(k)

and 

(k)

seq

It is easy to verify that the sequences v

(k)

, x

(k)

and 

(k)

seq have

following properties:

• x

(k+1) ∈ Y(x(k)

, 

(k)

seq);

• x

(k+1) minimizes f() over cl(Ye(v(k)

, x

(k)

, 

(k)

seq));

• 

(k)

seq < max

(v(k),x(k))

;

• ∀k, x

(k)

and v

(k)

satisfies the conditions of Lemma 1.

These properties will be used in the subsequent proof of Theorem 2.

D. Proof of Theorem 2

First note that by Theorem 1, Lemma 1 and Lemma 2, we have

that

• the inner loop on n in Algorithm 1 (lines 5-16) completes within

finite number of steps with a strictly feasible solution for the

constraints in (P0);

• f(x

(k)

) is an infinite monotonically non-increasing sequence if

without the lines 18-20 in the algorithm;

• the Jacobian matrix of h(·) evaluated at y

(k)

is invertible;

• ∀k, x

(k)

and x

(k+1) can be connected by a continuous path in

B(x

(k)

, max

(v(k),x(k))

), which must have a non-singular (hence

continuous) pre-image in the state space as guaranteed by

Theorem 1.

This proves items (i) and (ii) of the theorem.

Now, consider that f(·) is continuous and X is compact, we have

that the sequence f(x

(k)

) is bounded. By Monotone Convergence

Theorem, f(x

(k)

) converges. In other words, given ErrBound > 0,

∃M ∈ N such that





f(x

(k+1)) − f(x(k)

)





 < ErrBound, ∀k ≥ M.

Furthermore, by Lemma 2, x

∗

is a local minimum of f(·) over the

feasible set (12b)–(12d) whenever e = 0. This proves item (iii).

To prove item (iv) (a), note that the sequence {x

(k)

} stays in a

compact set, hence there exists a convergent subsequence. Let {σ(k)}

denote the indexes of a subsequence of {x

(k)

} with a limit point x

∗

,

namely x

σ(k) → x∗

as k → ∞. Moreover, observe that, in the

framework of the proposed Algorithm 1, the feasibility constraints

can be wrapped up as φ(xb, x) ≤ 0. In details, the continuous function

φ is defined as

φ(xb, x) := 

K(g(xb),xb)x + a(g(xb),xb) − V

max

V

min − K(g(xb),xb)x − a(g(xb),xb)



,

where g(·) is a continuous function defined by means of Theorem

1. Now we proceed by contradiction. Suppose that x

∗

is not a local

minimum of (P0). Therefore, it is not a minimizer of (P1(v

∗

, x

∗

, ))

for any  > 0. In particular, there exists xe such that φ(x

∗

, xe) < 0,

ξ(comps(x

∗ − xe)) < min, and

f(xe) < f(x

∗

). (28)

6 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
Now, from the fact that

(i) x

σ(k) → x∗

;

(ii) φ(·, x) is continuous for all x; and

(iii) ξ(compx(·)) is continuous;

(iv) hypothesis (17) holds;

there exists k0 such that for all k ≥ k0, we have that φ(x

σ(k)

, xe) ≤ 0

and ξ(comps(x

σ(k) − xe)) ≤ min ≤ 

σ(k)

seq . Therefore, xe ∈

cl(Ye(g(x

σ(k)

), x

σ(k)

, 

σ(k)

seq )). Since x

σ(k)+1 is an optimum in

cl(Ye(g(x

σ(k)

), x

σ(k)

, 

σ(k)

seq )) by its definition, it follows that

f(xe) ≥ f(x

σ(k)+1). (29)

However, since {f(x

(k)

)} is a monotonically non-increasing sequence, we have that f(x

σ(k)+1) ≥ f(x∗

). Combining this last

inequality with (28) and (29), we obtain that

f(x

σ(k)+1) ≥ f(x∗

) > f(xe) ≥ f(x

σ(k)+1),

a contradiction. Therefore, x

∗

is a local minimum of (P0).

For item (iv) (b), let {σ(k)} be a sequence such that 

σ(k)

seq → max

as k → ∞. Let 0 ∈ (0, max). There exists k0 such that for all

k ≥ k0, we have that 

σ(k)

seq > 0. Now let {σ

0

(k)} be a subsequence

of {σ(k)} such that x

σ

0

(k) → x∗

for some x

∗

. Then, the arguments

in proof for item (iv) (a) apply to x

∗

, by replacing {σ(k)} with

{σ

0

(k)} and min with 0. This completes the proof of theorem.

REFERENCES

[1] K. Lehmann, A. Grastien, and P. V. Hentenryck, “Ac-feasibility on tree

networks is np-hard,” IEEE Transactions on Power Systems, vol. 31,

no. 1, pp. 798–801, Jan 2016.

[2] J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow

problem,” IEEE Trans. Power Syst., vol. 1, no. 1, pp. 92–107, Feb.

2012.

[3] X. Bai, H. Wei, K. Fujisawa, and Y. Wang, “Semidefinite programming

for optimal power flow problems,” Int. J. of Electrical Power & Energy

Systems, vol. 30, no. 6–7, pp. 383–392, 2008.

[4] D. K. Molzahn and I. A. Hiskens, “Moment-based relaxation of the

optimal power flow problem,” in 2014 Power Systems Computation

Conference, Aug 2014, pp. 1–7.

[5] H. Abdi, S. D. Beigvand, and M. L. Scala, “A review of optimal power

flow studies applied to smart grids and microgrids,” Renewable and

Sustainable Energy Reviews, vol. 71, pp. 742–766, 2017.

[6] S. Paudyal, C. A. Canizares, and K. Bhattacharya, “Optimal operation

of distribution feeders in smart grids,” IEEE Trans. on Ind. Electron.,

vol. 58, no. 10, pp. 4495–4503, Oct. 2011.

[7] L. R. de Araujo, D. R. R. Penido, and F. de Alcntara Vieira, “A

multiphase optimal power flow algorithm for unbalanced distribution

systems,” International Journal of Electrical Power and Energy Systems,

vol. 53, pp. 632 – 642, 2013.

[8] S. Paudyaly, C. A. Canizares, and K. Bhattacharya, “Three-phase distribution OPF in smart grids: Optimality versus computational burden,”

in 2nd IEEE PES Intl. Conf. and Exhibition on Innovative Smart Grid

Technologies, Manchester, UK, Dec. 2011.

[9] E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal power

flow for smart microgrids,” IEEE Trans. Smart Grid, vol. 4, no. 3, pp.

1464–1475, Sep. 2013.

[10] L. Gan and S. H. Low, “Convex relaxations and linear approximation for

optimal power flow in multiphase radial networks,” in Power Systrem

Computation Conference, Wroclaw, Poland, August 2014.

[11] A. S. Zamzam, N. D. Sidiropoulos, and E. Dall’Anese, “Beyond Relaxation and Newton-Raphson: Solving AC OPF for Multi-phase Systems

with Renewables,” IEEE Transactions on Smart Grid, vol. PP, no. 99,

pp. 1–1, 2017.

[12] A. Bernstein, C. Wang, E. Dall’Anese, J.-Y. Le Boudec, and

C. Zhao, “Load-flow in multiphase distribution networks: Existence, uniqueness, and linear models,” 2017, [Online] Available at:

http://arxiv.org/abs/1702.03310.

[13] V. Frolov, P. G. Thakurta, S. Backhaus, J. Bialek, and M. Chertkov,

“Scalable optimal placement and sizing of facts devices using ac power

flow model,” 2016.

[14] V. Frolov and M. Chertkov, “Methodology for multi-stage, operationsand uncertainty-aware placement and sizing of facts devices in a large

power transmission system,” 2017.

[15] C. Wang, A. Bernstein, J.-Y. Le Boudec, and M. Paolone, “Existence and

uniqueness of load-flow solutions in three-phase distribution networks,”

IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 3319–3320,

2017.

[16] W. H. Kersting, Distribution System Modeling and Analysis. 2nd ed.,

Boca Raton, FL: CRC Press, 2007.

[17] R. Prada and L.J.Souza, “Voltage stability and thermal limit: Constraints

on the maximum loading of electrical energy distribution feeders,” IEE

Proceedings - Generation, Transmission and Distribution, vol. 145,

no. 5, pp. 573 – 577, 1998.

[18] K. Christakou, J.-Y. Le Boudec, M. Paolone, and D.-C. Tomozei,

“Efficient Computation of Sensitivity Coefficients of Node Voltages and

Line Currents in Unbalanced Radial Electrical Distribution Networks,”

IEEE Transactions on Smart Grid, vol. 4, no. 2, pp. 741–750, 2013.

[19] S. Bolognani and F. Dorfler, “Fast power system analysis via implicit

linearization of the power flow manifold,” in 2015 53rd Annual Allerton

Conference on Communication, Control, and Computing (Allerton), Sept

2015, pp. 402–409.

[20] C. Wang, A. Bernstein, J.-Y. Le Boudec, and M. Paolone, “Explicit

Conditions on Existence and Uniqueness of Load-Flow Solutions in

Distribution Networks,” IEEE Transactions on Smart Grid, vol. 9, no. 2,

pp. 953–962, 2018.

7 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.