Global Optimal Power Flow over Large-Scale Power Transmission NetworksI
Y. Shia , H. D. Tuana , P. Apkarianb , A. V. Savkinc
a School

of Electrical and Data Engineering, University of Technology, Sydney, NSW 2007, Australia
Control System Department, 2, avenua Edouard Belin, 31055 Toulouse, France
c School of Electrical Engineering and Telecommunications, the University of New South Wales, NSW 2052, Australia
b CERT-ONERA,

arXiv:1802.04456v1 [] 13 Feb 2018

Abstract
Optimal power flow (OPF) over power transmission networks poses challenging large-scale nonlinear optimization problems, which involve a large number of quadratic equality and indefinite quadratic inequality constraints. These computationally intractable constraints are often expressed by linear constraints plus matrix additional rank-one constraints on
the outer products of the voltage vectors. The existing convex relaxation technique, which drops the difficult rank-one
constraints for tractable computation, cannot yield even a feasible point. We address these computationally difficult
problems by an iterative procedure, which generates a sequence of improved points that converge to a rank-one solution.
Each iteration calls a semi-definite program. Intensive simulations for the OPF problems over networks with a few
thousands of buses are provided to demonstrate the efficiency of our approach. The suboptimal values of the OPF
problems found by our computational procedure turn out to be the global optimal value with computational tolerance
less than 0.01%.
Keywords: Optimal power flow (OPF) problem; large-scale transmission networks; rank-one matrix constraint;
nonsmooth optimization; semi-definite programming (SDP).

1. Introduction
Smart grids are operated by the advanced distribution management system (DMS), which is responsible for
supervisory control and data acquisition in reactive dispatch, voltage regulation, contingency analysis, capability
maximization and other smart operations. The optimal
power flow (OPF) problem, which determines a steady
state operating point that minimizes the cost of electric
power generation or the transmission loss is the backbone of DMS (see e.g. [1, 2, 3, 4] and references therein).
Mathematically, the OPF problem is highly nonlinear and
nonconvex due to numerous quadratic equality and indefinite quadratic inequality constraints for bus interconnections, hardware operating capacity and the balance between power demand and supply. These nonlinear constraints are mathematically troublesome so the state-ofthe-art nonlinear optimization solvers may converge to just
stationary points (see [5] and references therein), which
are not necessarily feasible. To handle these nonlinear
constraints, it is common to reformulate them as linear
constraints on the outer product W = V V H of the voltage vector V = (V1 , V2 , ..., Vn )T ∈ Cn . As a result, the
OPF problem is recast by a semi-definite program (SDP)
I Corresponding author H.D. Tuan. Tel. +61.02.9405 9435. Fax
+61.02.9405 9039
Email addresses: ye.shi@student.uts.edu.au (Y. Shi),
tuan.hoang@uts.edu.au (H. D. Tuan), Pierre.Apkarian@onera.fr
(P. Apkarian), a.savkin@unsw.edu.au (A. V. Savkin)

Preprint submitted to Journal of System & Control Letters

plus the additional rank-one constraint on outer product
matrix W [6, 7]. The rank-one constraint is then dropped
for semi-definite relaxation (SDR). However, the optimal
solution of SDR is of rank-more-than-one in general and
cannot help retrieval of a feasible point or stationary point
of the OPF problem [7, 8, 9]. In [10], we have extended the
technique of [11] for solving the beamforming optimization problems in signal processing to optimize the outer
product matrix W , which works very well and is practical for moderate-scale power distribution networks up to
n = 300 buses. There is another approach (see e.g. [12]
and references therein), which is based on hierarchies of
moment-based relaxation for nonconvex quadratic problems to tackle large networks with simple nonconvex constraints.
Power transmission networks in modern smart grids are often devised with a few thousand buses [13, 14, 15]. Under a
such large number n of buses it is impossible to use the single matrix W ∈ Cn×n , which involves n(n+1)/2 ≈ O(107 )
complex variables. On the other hand, the number of the
flow lines for bus connection is relatively moderate so only
a small portion of the crossed nonlinear terms Vk Vm∗ appears in the nonlinear constraints. The common approach
is to use the outer products of overlapped groups of the
voltage variables to cover them [16, 17, 9]. All rank-one
constraints on these outer products are then dropped for
SDR. Obviously, the optimal solution of this SDR usually is not of rank-one and thus does not have any physical meaning. There is no technique to retrieve a feasiFebruary 14, 2018

• I is the complex current vector, I = Y V = [I1 , I2 , · · ·
In ]T ∈ Cn , where Ik is the complex current injected
to bus k ∈ N .

ble rank-one point from the rank-more-than-one solution
of SDR.1 Multiple matrix rank constrained optimization
has received a great attention due to its potential application in robust control synthesis [19, 20] but to our best
knowledge there is no effective computation so far. The
contribution of this paper is two-fold:

• Ikm isX
the complex current
Xin the power line (k, m) ∈
L,
Ikm = Ik =
ykm Vm .
m∈N (k)

• An effective decomposition for large-scale OPF problems, which involves essentially reduced numbers of
the rank-one constraints on matrices of moderate
size for expressing the network nonlinear constraints;

m∈N (k)

• Skm = Pkm +jQkm is the complex power transferred
from bus k to bus m, where Pkm and Qkm represent
the real and reactive transferred power.
• SGk = PGk + jQGk is the complex power injected by
bus k ∈ G, where PGk and QGk represent the real
and reactive generated power.

• A new iterative procedure for rank-one constrained
optimization, which is practical for computational
solutions of large-scale indefinite quadratic programming. Simulations for the large-scale OPF problems
show that it is capable of finding the global optimal
solution with the computational tolerance less than
0.01%.

For each bus k, it is obvious that
SGk − SLk

The paper is structured as follows. Section 2 is devoted
to the OPF problem formulation and its difficulties. Its
computational solution is developed in Section 3. Section
4 provides simulation to show the efficiency of our method.
The conclusions are drawn in Section 5.
Notation. j denotes the imaginary unit; M  0 means
that M is a Hermitian symmetric positive semi-definite
matrix; rank(M ) is the rank of the matrix M ; <(·) and =(·)
denote the real and imaginary parts of a complex quantity;
a ≤ b for two complex numbers a and b is componentwise
understood, i.e. <(a) ≤ <(b) and =(a) ≤ =(b); h., .i is the
dot product of matrices, while diag{Ai } denotes the matrix
with diagonal blocks Ai and zero off-diagonal blocks; the
cardinality of a set L is denoted by |L|.

= (PGk − PLk )X
+ j(QGk − QLk )
∗
= Vk Ik∗ = Vk
Vm∗ ykm
.
m∈N (k)

Therefore, the real generated power PGk and reactive generated power QGk at bus k are the following nonconvex
quadratic functions of the bus voltage vector
Xvariable V :=
∗
T
n
(V1 , V2 , ..., Vn ) ∈ C : PGk = PLk + <(
Vk Vm∗ ykm
)
m∈N (k)

and QGk = QLk + =(

X

∗
Vk Vm∗ ykm
).

m∈N (k)

Bus 1

2. Optimal power flow problem and challenges

G

y13

Consider an AC electricity transmission network with a
set of n buses N := {1, 2, · · · , n}. The buses are connected
through a set of flow lines L ⊆ N × N , i.e. bus m is
connected to bus k if and only if (m, k) ∈ L. Accordingly,
N (k) := {m ∈ N : (m, k) ∈ L}. The power demanded at
bus k ∈ N is SLk = PLk + jQLk , where PLk and QLk are
the real and reactive power. A subset G ⊆ N of buses is
supposed to be connected to generators. Any bus k ∈ N \G
is thus not connected to generators.
Other physical parameters are following [21, 22, 23]:

y 23

Bus 3
y 35

Bus 5

• Y = [ykm ](k,m)∈N ×N ∈ Cn×n is the admittance matrix [23]. Each ykm is the mutual admittance between bus k and bus m, so ykm = ymk ∀ (k, m) ∈ L.

y12

y45

G

Bus 2

y24
Bus 4

Figure 1: WB5-Five Bus Network [24]

For illustrative purpose, Figure 1 provides a diagram
of WB5-Five Bus Network [24] with N = {1, 2, 3, 4, 5}.
It is seen from the bus connection in this Figure that
N (1) = {2, 3} because buses 2 and 3 are connected to
bus 1, N (2) = {1, 3, 4}, N (3) = {1, 5}, N (4) = {2, 5} and
N (5) = {3, 4}. Also G = {1, 5} because buses 1 and 5 are
connected to generators.

• V is the complex voltage vector, V = [V1 , V2 , · · · , Vn ]T
∈ Cn , where Vk is the complex voltage injected to
bus k ∈ N .
1 There is an algorithm of finding a rank-one solution [18, Alg. 1],
which however is applicable to simple nonconvex constraints and is
not guaranteed to convergence

2

problem in W = [Wkm ]k,m=1,...,n ∈ Cn×n :

The objective of OPF is to minimize either the total
cost of the generated power defined by [22]
X
f (PG ) =
(ck2 PG2 k + ck1 PGk + ck0 ),
(1)

min
X

W ∈Cn×n

−PLk − jQLk =

k∈G

W  0,

(6a)

∗
Wkm ykm
, k ∈ N \ G,

(6b)

∗
Wkm ykm
) ≤ PGmax
,
k

(6c)

F (W )

s.t.

m∈N (k)

with given ck2 > 0, ck1 and ck0 and real active generated
power PG , or the total transmission losses defined by
X
X
Floss (PG ) =
(PGk − PLk ) −
PL k .
(2)

PGmin
≤ PLk + <(
k

m∈N (k)

k ∈ G,

k∈N \G

k∈G

Qmin
Gk

X

∗
Vk Vm∗ ykm
)) + ck0 ].

∗
Wkm ykm
)

≤ Qmax
Gk ,

(6d)

k ∈ G,
(Vkmin )2
Wkk +

≤ Wkk ≤ (Vkmax )2 , k ∈ N ,
∗
max
|Wkm ykm
| ≤ Skm
, (k, m) ∈ L,
max 2
Wmm − Wkm − Wmk ≤ (Vkm
) ,

(6e)
(6f)
(6g)

(k, m) ∈ L,

m∈N (k)

+ck1 (PLk + <(

≤ QLk + =(

X
m∈N (k)

In this paper, we consider only the objective function (1).
It is quite clear that our approach can be adjusted to solve
the objective (2) in a straightforward manner.
Write the objective (1) as the following function of the bus
voltages V :
X
X
∗
f (V ) =
[ck2 (PLk + <(
Vk Vm∗ ykm
))2
k∈G

X

(3)

max
=(Wkm ) ≤ <(Wkm ) tan θkm
, (k, m) ∈ L,

(6h)

rank(W ) = 1,

(6i)

m∈N (k)

X

where F (W ) =

Accordingly, the OPF problem is formulated as follows

[ck2 (PLk + <(

k∈G

minn f (V )

V ∈C

X

−PLk − jQLk =

∗
Vk Vm∗ ykm
,k

s.t. (4a)

+ ck1 (PLk + <(

≤ PLk + <(

X
X

∗
Vk Vm∗ ykm
) ≤ PGmax
, k ∈ G (4c)
k
∗
Vk Vm∗ ykm
) ≤ Qmax
Gk , k ∈ G (4d)

m∈N (k)

Vkmin ≤ |Vk | ≤ Vkmax , k ∈ N , (4e)
∗
max
|Skm | = |Vk Vm∗ ykm
| ≤ Skm
, ∀(k, m) ∈ L (4f)
max
|Vk − Vm | ≤ Vkm
, (k, m) ∈ L, (4g)
max
|arg(Vk ) − arg(Vm )| ≤ θkm
, (k, m) ∈ L, (4h)

where (4b) is the equation of the balance between the demand and supply power at bus k ∈ N \ G, (4c)-(4d) are
the power generation bounds with the lower bound PGmin
,
k
min
max
max
QGk and upper bound PGk , QGk of the real power reactive power generations, (4e) are the voltage amplitude
bounds, (4f)-(4h) are capacity limitations with the line
currents between the connected buses constrained by (4f),
while (4g)-(4h) guarantee the voltage balance in terms of
their magnitude and phases [23].
One can see that (4) is a highly nonlinear optimization subject to nonlinear inequality and equality constraints and
thus poses a real challenge for its computation. A common
approach is to use the slack variables
Wkm = Vk Vm∗ , k = 1, ..., n; m = 1, ..., n

m∈N (k)
∗
Wkm ykm
)) + ck0 ], which is convex

quadratic in Wkm , while all constraints (6b)-(6h) are linear. All the problem nonconvexity is now concentrated in
the single rank-one constraint (8b), which is then dropped
for SDR. If the solution of this SDR is of rank-one, then
it obviously leads to the global solution of the nonconvex
optimization problem (6). Otherwise even a feasible point
of (6) is hardly retrieved from the SDR solution. Moreover formulation (6) cannot be practically used for n more
than a few thousands in large-scale networks as it involves
n(n + 1)/2 variables–a prohibitively large number. The
next section will address and resolve all these issues.

m∈N (k)

Qmin
Gk ≤ QLk + =(

∗
Wkm ykm
))2

m∈N (k)

∈ N \ G, (4b)

m∈N (k)

PGmin
k

X

X

3. Nonsmooth optimization based solution
The first issue is to decompose large-size matrix W
in (6) into matrices of smaller size to make the problem
tractable even with limited computational power. This is
also prompted by the fact that there is only a small portion of the crossed terms Vk Vm∗ appearing in the nonlinear
constraints (4b)-(4h) so the large-size matrix variable W ∈
Cn×n contains many redundant terms Vk Vm∗ . The main result of [16, 17, 9] is to decompose the set N := {1, 2, ..., n}
of buses into I overlapped subsets Ni = {i1 , ..., iNi } of
buses, called bags, such that i` ∈ N (i`+1 ), ` = 1, ..., Ni −
1 and iNi ∈ N (i1 ), for each i = 1, 2, ..., I, i.e. the buses
in the same bag are serially connected. The set of bags
can be reset to make bags of relatively same size. Define
the Hermitian symmetric matrix variables

(5)

and then recast the problem (4) in V ∈ Cn to the following

W i = [Wik im ]k,m=1,..,Ni ∈ CNi ×Ni , i = 1, 2, ..., I.
3

(7)

By replacing Wkm = Vk Vm∗ in (6) we have the following
equivalent reformulation for (6)
min

W =diag{W i }

F (W )

s.t.

(6b) − (6h),

Indeed, (8b) implies Trace(W i ) − λmax (W i ) ≥ 0 ∀ i, so (9)
means that Trace(W i ) = λmax (W i ), i.e. W i has only one
nonzero eigenvalue
so it is of rank-one. The nonnegative
PI
i
quantity
(Trace(W
) − λmax (W i )) can therefore be
i=1
used to measure the degree of satisfaction of the rankone constraints (9). Without squaring, the penalization
PI
i
i
i=1 (Trace(W ) − λmax (W )) is exact, meaning that (9)
can be satisfied by a minimizer of the problem

(8a)

W i  0, i = 1, ..., I,

(8b)

rank(W i ) = 1, i = 1, ..., I.

(8c)

Reference [16, 17] just dropped all rank-one constraints
in (8c) for SDR without any justification. Reference [9]
also dropped all rank-constraints in (8c) but then used a
penalized SDR for locating low-rank semi-definite matrices
W i in (7). Based on these low-rank matrices, [9, Sec. IV]
also proposed to find rank-one matrices, which however
are not necessarily feasible to (8).
PI
The variable number in (8) is i=1 Ni (Ni + 1)/2. To keep
this number reasonably moderate, it is desired that both
I and Ni are sufficiently moderate. However, one can see
that the above described decomposition [16, 17, 9] leads
to a large number I of bags as well as few large size Ni
that result in many rank-one constraints in (8c), which are
much less probably satisfied by solving SDR.
Our first step toward to computation of (4) is to develop a new decomposition with many fewer bags involved.
Recalling that N (k) is the set of the buses that are connected to bus k, the cardinality |N (k)| is small in largescale networks. We resort N = {1, 2..., n} as N = {N1 , ...,
Nn } such that the cardinality |N (Nk )| is in decreased order:
|N (N1 )| ≥ |N (N2 )| ≥ ... ≥ |N (Nn )|.

min

W =diag{W i }

−λmax (W ))

(Trace(W i )

s.t.

(6b) − (6h), (8b),

(10)

with a finite value of µ > 0 (see e.g. [25, Chapter 16]).
This is generally considered as a sufficiently nice property
to make such exact penalization attractive.
For any W i,(κ) feasible for the convex constraints (6b)(6h), (8b), function λmax (W i ) is nonsmooth and is lower
bounded by
i,(κ) H
i,(κ)
λmax (W i ) = max wH W i w ≥ (wmax
) W i wmax
,
||w||=1

(11)

i,(κ)

where wmax is the normalized eigenvector corresponding
to the eigenvalue λmax (W i,(κ) ), i.e.
i,(κ) H
i,(κ)
λmax (W i,(κ) ) = (wmax
) W i,(κ) wmax
.

(12)
i,(κ)

Accordingly, µλmax (W i )−µλmax (W i,(κ) ) ≥ µ((wmax )H W i
i,(κ)
i,(κ)
i,(κ)
i,(κ)
i,(κ)
.wmax − (wmax )H W i,(κ) wmax ) = µhwmax (wmax )H , W i −
i,(κ)
i,(κ)
W i,(κ) i, so µwmax (wmax )H is a subgradient of the function
i,(κ)
i,(κ)
µλmax (W i ) at W i,(κ) . Then µdiag{wmax (wmax )H } is a
PI
subgradient of the function µ i=1 λmax (W i ) at diag{W i,(κ) }.

N2 = {i ∈ N (N2 ) : {i, N2 } 6⊂ N1 }.
Note that the crossed term Vi VN∗2 is already treated in the
previous bag N1 whenever {i, N2 } ⊂ N1 so we exclude
such bus i in defining bag N2 .
Similarly, for ` ≥ 3 the `-th bag is defined as

The following SDP provides an upper bound for the
nonconvex optimization problem (10)
min

0

N` = {i ∈ N (N` ) : {i, N` } 6⊂ N`0 ∀1 ≤ ` ≤ ` − 1}

W =diag{W i }

F (κ) (W ) := F (W ) + µ

i,(κ) H
i,(κ)
−(wmax
) W i wmax
)

Vi VN∗`

to exclude those buses i, whose crossed term
already
is treated in a previous bag.
As each |Ni | is obviously small, such decomposition is very
efficient, leading to a substantial reduction of involved bags
in comparison with that used in [16, 17, 9].
Our next step is to tackle the numerous difficult rankone constraints in (8c), not dropping them for SDR as in
all the previous works.
Firstly we express I rank-one constraints in (8c) by the
following single spectral constraint
(Trace(W i ) − λmax (W i )) = 0,

I
X
i=1

i

Accordingly, the first bag of buses is defined as N1 =
N (N1 ). The second bag is defined as

I
X

Fµ (W ) := F (W ) + µ

I
X

(Trace(W i )

i=1

s.t.

(6b) − (6h), (8b) (13)

because F (κ) (diag{W i }) ≥ Fµ (diag{W i }) ∀ W i  0 according to (11). Suppose that W (κ+1) = diag{W i,(κ+1) } is
the optimal solution of SDP (13). Since W (κ) = diag{W i,(κ) }
is also feasible to (13) with Fµ (W (κ) ) = F (κ) (W (κ) ), it is
true that
Fµ (W (κ+1) ) ≤ F (κ) (W (κ+1) ) ≤ F (κ) (W (κ) ) = Fµ (W (κ) ),
so W (κ+1) is a better feasible point of (10) than W (κ) . Initialized by any feasible point W (0) = diag{W i,(0) } of SDP
constraint (8a)-(8b), the sequence {W (κ) } = {diag{W i,(κ) }}
with W (κ+1) = diag{W i,(κ+1) } iteratively generated as the
optimal solution of SDP (13) is a sequence of improved feasible points of the nonconvex optimization problem (10).

(9)

i=1

where λmax (W i ) stands for the maximal eigenvalue of W (i) .
4

Since W (κ) are uniformly bounded, the sequence {W (κ) }
has a limit point W̄ = diag{W̄ i }, which is the optimal
solution of the optimization problem
min

W =diag{W i }

F (W ) + µ

i
i
−(w̄max
)H W i w̄max
)

I
X

iteration may turn to be more than one with
Trace(W i,(κ+1) ) − λmax (W i,(κ+1) ) >
Trace(W i,(κ) ) − λmax (W i,(κ) ).

(Trace(W i )

Consequently, it is very difficult to achieve rank-one for all
W i,(κ) as desired. It is also impossible to add a ”weight” to
each term under the sum in the objective in (13) to handle
the individual convergence of Trace(W i ) − λmax (W i ).
We now develop a systematic way to resolve this issue as
follows. For κ = 0, 1, ..., and W (κ) = diag{W (i,(κ)) } define

i=1

s.t.

(6b) − (6h), (8b),

(14)

i
where w̄max
is the normalized eigenvector corresponding
to the eigenvalue λmax (W̄ i ) of W̄ i . Particularly,

F (W ) + µ

I
X

i
i
(Trace(W i ) − (w̄max
)H W i w̄max
)

L(κ) = {i ∈ {1, ..., I} : rank(W i,(κ) ) = 1}

≥

and generate W (κ+1) = diag{W i,(κ+1) } as the optimal solution of the following SDP instead of SDP (13)

i=1

F (W ) + µ

I
X

i
i
(Trace(W̄ i ) − (w̄max
)H W̄ i w̄max
),

i=1

min

or equivalently, under the definition g(W ) = F (W ) +
I
X
µ
Trace(W i ),

W =diag{W i }

Note that Trace(W i ) ≥ wH W w for all ||w|| = 1 and it
is obvious that rank(W i ) = 1 if and only if Trace(W i ) −
H
W i wmax = 0 for some normalized wmax . Therefore,
wmax
the constraint (18b) for some tolerance tol is introduced
to warrant the rank-one of all W i,(κ+1) , i ∈ L(κ) . As a
result L(κ) ⊂ L(κ+1) and L(κ) → {1, ..., I} is expected to
have all W i,(κ) of rank-one. Unlike (13), the iterations (18)
leads to achieving rank-one of all W i while the objective
function Fµ is still decreased.
In summary, we propose the following Large-Scale Nonsmooth Optimization Algorithm (Large-scale NOA) for the
multiple rank-one constrained optimization problem (10).

i
i
g(W ) − g(W̄ ) − hµdiag{w̄max
(w̄max
)H }, W − W̄ i ≥ 0

for all feasible points W = diag{W i } in (6b)-(6h), (8b).
As a result, W̄ is the optimal solution of the convex optimization problem
min

i
i
g(W ) − hµdiag{w̄max
(w̄max
)H }, W − W̄ i

s.t.

(6b) − (6h), (8b),

so it must satisfy the optimality condition
i
i
h∇g(W̄ ) − µdiag{w̄max
(w̄max
)H }, W − W̄ i ≥ 0

for all feasible points diag{W i } in (6b)-(6h), (8b). The
latter is also the first order necessary optimality condition
i
i
(w̄max
)H } is a subgradient of
for (10) because µdiag{w̄max
PI
i
the function µ i=1 λmax (W ) at W̄ . As our simulations
will show, W̄ is indeed the global optimal solution of (10)
and (8).
However, unlike [10] with only a single rank-one constrained matrix, although quantity
I
X
(Trace(W i,(κ) ) − λmax (W i,(κ) ))

Initialization. Solve SDP
min

W =diag{W i }

F (W )

s.t.

(6b) − (6h), (8b)

(19)

to generate W (0) := diag{W i,(0) }. If rank(W i,(0) ) ≡ 1
stop: W (0) is the global solution of the nonconvex optimization problem (8). Otherwise set κ = 0 and define
L(κ) by (17).
κ-th iteration. For κ = 0, 1, .., solve (18) to generate
W (κ+1) := diag{W i,(κ+1) }. Reset κ = κ + 1 and define
L(κ) by (17). Stop whenever L(κ) = {1, ..., I}. Otherwise
go to the next iteration.

(15)

i=1

in (10) is iteratively decreased, not all individual quantities
Trace(W i,(κ) ) − λmax (W i,(κ) )

I
X
F (W ) + µ [Trace(W i )

i=1
i,(κ) H
i i,(κ)
−(wmax ) W wmax ] s.t. (6b) − (6h), (8b), (18a)
i,(κ) H
i,(κ)
Trace(W i ) − (wmax
) W i wmax
≤ tol , i ∈ L(κ) .(18b)

i=1

W =diag{W i }

(17)

4. Simulation results

(16)

The computation facilities for our implementation are
Processor Intel(R) Core i5-3470 CPU @3.20GHz, Matlab
version R2013b and CVX with SDPT3. We set the tolerances  = tol = 10−5 and the penalty
µ = 106
PI parameter (i,(0))
which makes the penalty term µ i=1 (Trace(W
)−
λmax (W i,(0) )) at similar magnitude with the objective F (W (0) ).
The data source for all examples is Matpower version 5.1

are iteratively decreased so the rank of each matrix W i,(κ)
is no longer iteratively reduced to one as expected. Worse,
W i,(κ) is rank-one but the rank of W i,(κ+1) in the next

5

[23]. All examples were considered in [16, 17, 9] by SDR
only. We recall that T is the number of matrix variables
W i in the OPF problem (8) and L(κ) is defined by (17)
is the set of rank-one matrices W i found after κ-th iteration. The capability of our large-scale NOA in locating
the global optimal solution of the OPF problems is demonstrated by showing that the global optimality tolerance
(GOT) of its found solution defined as

4.4. Polish-2746wop system
There are n = 2746 buses, 514 generators and 3514
transmission lines, which lead to 2232 nonlinear constraints
in (4b).
Initialization. A lower bound 1.2039 × 106 of (8) is obtained by solving SDP (19). |L(0) | = 1538 and there are 8
matrices W i,(0) of rank-more-than-one. Their largest size
(smallest size, resp.) is 6 × 6 (3 × 3, resp.).
Stage 1. |L(2) | = I = 1546 is achieved. The found value
of the objective is 1.2040 × 106 with GOT 8.3063e − 05.

the found value- lower bound
lower bound

4.5. Polish-2746wp system

is almost zero.
The numerical examples are presented as follows.

There are n = 2746 buses, 520 generators and 3514
transmission lines, which lead to 2226 nonlinear constraints
in (4b).
Initialization. A lower bound 1.626590 × 106 of (8) is obtained by solving SDP (19). |L(0) | = 1545 and there are 2
matrices W i,(0) of rank-more-than-one. Their size is 4 × 4.
Stage 1. |L(1) | = I = 1547 is achieved. The found value
of the objective is 1.626591 × 106 with GOT 6.1478e − 07.

4.1. Polish-2383wp system
There are n = 2383 buses, 327 generators and 2896
transmission lines, leading to 2056 nonlinear constraints
in (4b).
Initialization. A lower bound 1.8490×106 of (8) is found by
solving SDP (19). |L(0) | = 1210 and there are 32 matrices
W i,(0) of rank-more-than-one. Their largest size (smallest
size, resp.) is 10 × 10 (2 × 2, resp.).
Stage 1. |L(10) | = 1234 is achieved. There are 8 matrices
W i,(10) of rank-more-than-one. Their largest size (smallest
size, resp.) is 10 × 10 and (3 × 3, resp.).
Stage 2. |L(19) | = 1237 is achieved. There are 5 matrices
W i,(20) of rank-more-than-one. Their largest size (smallest
size, resp.) is 9 × 9 (3 × 3, resp.).
Stage 3. |L(25) | = I = 1242 is achieved. The found value
of the objective is 1.8408 × 106 with GOT 4.3267e − 04.

4.6. Polish-3012wp system
There are n = 3012 buses, 502 generators and 3572
transmission lines, which lead to 2510 nonlinear constraints
in (4b).
Initialization. A lower bound 2.5717 × 106 of (8) is obtained by solving SDP (19). |L(0) | = 1682 and there are 7
matrices W i,(0) of rank-more-than-one. Their largest size
(smallest size, resp.) is 7 × 7 (2 × 2, resp.).
Stage 1. |L(4) | = I = 1689 is achieved. The found value
of the objective is 2.5727 × 106 with GOT 3.8885e − 04.

4.2. Polish-2736sp system
4.7. Polish-3120sp system

There are n = 2736 buses, 420 generators and 3504
transmission lines, which lead to 2316 nonlinear constraints
in (4b).
Initialization. A lower bound 1.3041 × 106 of (8) is obtained by solving SDP (19). |L(0) | = 1534 and there are 4
matrices W i,(0) of rank-more-than-one. Their largest size
(smallest size, resp.) is 6 × 6 (4 × 4, resp.).
Stage 1. |L(9) | = I = 1538 is achieved. The found value
of the objective is 1.3042 × 106 with GOT 7.6681e − 05.

There are n = 3120 buses, 505 generators and 3693
transmission lines, which lead to 2615 nonlinear constraints
in (4b).
Initialization. A lower bound 2.1314 × 106 of (8) is obtained by solving SDP (19). |L(0) | = 1749 and there are 8
matrices W i,(0) of rank-more-than-one. Their largest size
(smallest size, resp.) is 8 × 8 (2 × 2, resp.).
Stage 1. |L(9) | = I = 1757 is achieved. The found value
of the objective is 2.1391 × 106 with GOT 0.0036.

4.3. Polish-2737sop system
4.8. Numerical summary

There are n = 2737 buses, 399 generators and 3506
transmission lines, which lead to 2338 nonlinear constraints
in (4b).
Initialization. A lower bound 7.7571 × 105 of (8) is obtained by solving SDP (19). |L(0) | = 1532 and there are 6
matrices W i,(0) of rank-more-than-one. Their largest size
(smallest size, resp.) is 6 × 6 (3 × 3, resp.).
Stage 1. |L(2) | = I = 1538 is achieved. The found value
of the objective is 7.7572 × 105 with GOT 1.2891e − 05.

One can observe that GOT of the solutions computed
by the large-scale NOA is very small, proving its capability
to provide the global solution of (8). Table 1 and Table
2 summarize the main points in our simulation. The second and third columns of Table 1 are the number I of
bags in (7) by our decomposition and by that in [9], while
the fourth and fifth columns give the maximum size Ni
in (7). One can see that both I and the maximum Ni
by our decomposition are substantially smaller than their
counterparts by [9]. This leads to far smaller numbers of
variables in (8), which are provided in the seventh and
6

Table 1: Comparison of bags number I, largest bag size Mi and number of variables
System
Polish-2383wp
Polish-2736sp
Polish-2737sop
Polish-2746wop
Polish-2746wp
Polish-3012wp
Polish-3120sp

I
1242
1538
1538
1546
1547
1689
1757

I by [9]
2383
2736
2737
2746
2746
3012
3120

Max. Ni
10
10
10
10
10
10
10

Max. Ni by [9]
23
23
23
23
24
24
24

Var. # in (6)
2,840,536
3,744,216
3,746,953
3,771,631
3,771,631
4,537,578
4,868,760

Var. # in (8)
23,199
27,298
27,034
29,024
28,257
30,996
32,637

Var. # in (8) by [9]
89,893
104,388
103,720
108,950
107,148
116,799
121,869

Found value
1.8408 × 106
1.3042 × 106
7.7572 × 105
1.2040 × 106
1.6266 × 106
2.5727 × 106
2.1391 × 106

Found by [9]
1.8742 × 106
1.3082 × 106
7.7766 × 105
1.2085 × 106
1.6324 × 106
2.6089 × 106
2.1608 × 106

Table 2: Performance comparison
System
Polish-2383wp
Polish-2736sp
Polish-2737sop
Polish-2746wop
Polish-2746wp
Polish-3012wp
Polish-3120sp

Found value
1.8408 × 106
1.3042 × 106
7.7572 × 105
1.2040 × 106
1.6266 × 106
2.5727 × 106
2.1391 × 106

Found by [9]
1.8742 × 106
1.3082 × 106
7.7766 × 105
1.2085 × 106
1.6324 × 106
2.6089 × 106
2.1608 × 106

eighth columns. The number n(n + 1)/2 of complex variables in (6) is also provided in the sixth column to contrast
to the number of complex variables in (8) in the seventh
column. Furthermore, the second column of Table 2 provides the best values of (4) found by our large-scale NOA,
which are far smaller than ones in the third and fourth
columns found by [9] and Matpower6.0 [23] (using an interior point method), respectively. In short, our computation approach to the OPF problem (8) outperforms other
existing approaches in terms of computational efficiency
and performance.

Found by [23]
1.8685 × 106
1.3078 × 106
7.7763 × 105
1.2083 × 105
1.6317 × 106
2.5917 × 106
2.1427 × 106

[4] K. S. Pandya, S. K. Joshi, A survey of optimal power flow methods, J. of Theoretical and Applied Information Technology 4 (5)
(2008) 450–458.
[5] W. Bukhsh, A. Grothey, K. McKinnon, P. Trodden, Local solutions of the optimal power flow problem, IEEE Trans. Power
Systems 28 (4) (2013) 4780–4788.
[6] X. Bai, H. Wei, K. Fujisawa, Y. Wang, Semidefinite programming for optimal power flow problems, Elect. Power Energy
Syst. 30 (6-7) (2008) 383–392.
[7] J. Lavaei, S. H. Low, Zero duality gap in optimal power flow
problem, IEEE Trans. Power Systems 27 (1) (2012) 92–107.
[8] R. Madani, S. Sojoudi, J. Lavaei, Convex relaxation for optimal power flow problem: Mesh networks, IEEE Trans. Power
Systems 30 (2015) 199–211.
[9] R. Madani, M. Ashraphijuo, J. Lavaei, Promises of conic relaxation for contigency-constrained optimal power flow problem,
IEEE Trans. Power Systems 30.
[10] Y. Shi, H. D. Tuan, S. Su, H. H. M. Tam, Nonsmooth optimization for optimal power flow over transmission networks, in:
Proc. 2015 IEEE Global Conf. on Signal and Information Processing (GlobalSIP), 2015.
[11] H. A. Phan, H. D. Tuan, H. H. Kha, D. T. Ngo, Nonsmooth
optimization for efficient beamforming in cognitive radio multicast transmission, IEEE Trans. Signal Processing 60 (6) (2012)
2941–2951.
[12] C. Josz, D. K. Molzahn, Moment/sum-of-squares hierarchy for complex polynomial optimization,
Available:http://arxiv.org/abs/1508.02068.
[13] Bamberger, et al, Vision and Strategy for Europes Electricity
Networks of the Future, Office for Official Publications of the
European Communities, 2006.
[14] P. Hallberg, et al, Smart Grids and Networks of the Future EURELECTRIC Views, Ref: 2009-030-0440, 2009.
[15] M. Hashmi, S. Hnninen, K. Mki, Survey of smart grid concepts,
architectures, and technological demonstrations worldwide, in:
Proc. 2011 IEEE PES Conf. on Innovative Smart Grid Technologies (ISGT Latin America), 2011, pp. 1–7.
[16] D. K. Molzahn, J. T. Holzer, B. C. Lesieutre, C. L. DeMarco,
Implementation of a large-scale optimal power flow solver based
on semidefinite programming, IEEE Trans. Power Systems 28
(2013) 3987–3998.
[17] M. S. Andersen, A. Hansson, L. Vandenberghe, Reducedcomplexity semidefinite relaxations of optimal power flow problems, IEEE Trans. Power Systems 29 (4) (2014) 1855–1863.
[18] D. K. Molzahn, C. Josz, I. A. Hiskens, P. Panciatici, A
Laplacian-based approach for finding near globally optimal solutions to OPF problems, IEEE Trans. Power Systems 32 (1)
(2017) 305–315.

5. Conclusion
The OPF problems over power transmission network
are large-scale optimization problems, which involve a large
number of quadratic equality and indefinite quadratic inequality constraints and thus are difficult computationally.
We have developed a large-scale nonsmooth optimization
algorithm to compute their optimal solutions, which is
efficient and practical for large-scale power transmission
networks of a few thousands of buses. Applications of
the developed large-scale NOA to the OPF problems over
three-phase power transmission networks are currently under investigation.
Acknowledgements. We thank Cedric Josz, a coauthor of [12] and [18] for sending us these references.
References
[1] J. Carpentier, Contribution to the economic dispatch problem,
Bulletin Society Francaise Electriciens 3 (8) (1962) 431–447.
[2] M. Huneault, F. D. Galiana, A survey of the optimal power flow
literature, IEEE Trans. Power Systems 6 (2) (1991) 762–770.
[3] J. A. Momoh, M. E. El-Hawary, R. Adapa, A review of selected
optimal power flow literature to 1993. part I: Nonlinear and
quadratic programming approaches, Solar Energy 14 (1) (1999)
96–104.

7

[19] P. Apkarian, H. D. Tuan, A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems, in: Proc. of the
38-th Conf. on Decision and Control (CDC), Phoenix, Arizona,
1999, pp. 2238–2245.
[20] P. Apkarian, D. Noll, O. Prot, A trust region spectral bundle
method for nonconvex eigenvalue optimization, SIAM J. Optimization 19 (1) (2008) 281–306.
[21] S. Hadi, Power system analysis, PSA, 2010.
[22] H. W. Dommel, W. F. Tinney, Optimal power flow solutions,
IEEE Trans. Power Apparatus Systems PAS-87 (10) (1968)
1866–1876.
[23] R. D. Zimmerman, C. E. Murillo-Sanchez, R. J. Thomas, MATPOWER: steady-state operations, planning, and analysis tools
for power systems research and education, IEEE Trans. Power
Systems 26 (1) (2011) 12–19.
[24] W. Bukhsh, A. Grothey, K. McKinnon, P. Trodden,
Test case archive of optimal power flow (OPF) problems
with local optima, http://www.maths.ed.ac.uk/optenergy/
LocalOpt/introduction.html.
[25] J. F. Bonnans, J. C. Gilbert, C. Lemarechal, C. Sagastizabal, Numerical Optimization Theoretical and Practical Aspects
(second edition), Springer, 2006.

8

